seedr: Hydro and Thermal Time Seed Germination Models in R

Eduardo Fernández-Pascual

2020-10-25

Overview

seedr is an R package that provides functions to fit hydro and thermal time germination models. These models characterize seed lots by two sets of parameters: (i) the physiological thresholds (water, temperature) between which the seed lot can germinate, and (ii) the physiological-time units that the seed lot needs to accumulate before it can germinate. seedr fits the hydro time model of Bradford (Gummerson 1986, Bradford 1990, Bewley et al. 2013) and the thermal time model of Garcia-Huidobro (Garcia-Huidobro et al. 1982, Gummerson 1986, Bewley et al. 2013). It allows to fit models to grouped datasets, i.e. datasets containing multiple species, seedlots or experiments.

Installation

Install the last version of seedr from GitHub by running the following code in your R session. Note that you will need to install the package devtools if you don’t have it already:

install.packages("devtools")
devtools::install_github("efernandezpascual/seedr")

Using seedr

# once installed, load seedr
library(seedr)

To operate seedr you need just one function: physiotime(), a wrapper function that formats your data and fits a hydro/thermal time model. In this vignette, first we will see the basic operations behind physiotime(), and then we will see how to use physiotime().

physiodata, a way of storing germination data

The first step of a seedr analysis is processed by the physiodata() function, which will create a physiodata object storing the germination data in a format adequate for subsequent steps. Importantly, the physiodata can be used for basic visualization of your germination data.

Let’s say you have imported your data into R and saved it into a data.frame object called centaury, which looks like this:

head(centaury, 3)
#>                species population temperature dish times germinated germinable
#> 1 Centaurium somedanum   La Malva          20   A1     1          0         22
#> 2 Centaurium somedanum   La Malva          20   A1     2          0         22
#> 3 Centaurium somedanum   La Malva          20   A1     3          0         22

In your data, each column must be a variable and each row a record. Important: for each germination scoring time, you need to record the number of seeds that germinated at that time, not the cumulative germination count since the start of the experiment.

Now, you can transform your data.frame to a physiodata object like this:

cent <- physiodata(d = centaury, # the name of your data object
                   t = "times", # the column with the scoring times
                   x = "temperature", # the column with the experimental treatment,
                   g = "germinated", # the column with the number of germinated seeds,
                   pg = "germinable") # the column with the total number of viable seeds

You can use summary() on your physiodata object to get a table in which you will see, for each treatment, the final germination proportions (with 95% binomial confidence interval) and the median germination rate (i.e., the inverse of the time taken for 50% germination, 1/t50):

summary(cent)
#>     treatment time germination.mean germination.lower germination.upper
#>  1:       0.0   28        0.0000000         0.0000000        0.07707356
#>  2:       3.0   28        0.0000000         0.0000000        0.07707356
#>  3:       4.4   28        0.0000000         0.0000000        0.08029560
#>  4:       5.8   28        0.0000000         0.0000000        0.07410013
#>  5:       7.3   28        0.2391304         0.1391218        0.37935139
#>  6:       8.7   28        0.7083333         0.5682074        0.81758421
#>  7:      10.1   28        0.8888889         0.7650090        0.95159528
#>  8:      11.5   28        0.9583333         0.8602434        0.98849811
#>  9:      12.9   28        0.9347826         0.8249726        0.97757214
#> 10:      14.3   28        0.9795918         0.8930648        0.99638833
#> 11:      15.8   28        0.9591837         0.8628697        0.98873437
#> 12:      17.2   28        0.9800000         0.8950456        0.99646074
#> 13:      18.6   28        1.0000000         0.9273022        1.00000000
#> 14:      20.0   28        0.9333333         0.8214340        0.97706797
#> 15:      25.0   28        0.2380952         0.1348096        0.38527544
#> 16:      30.0   28        0.0000000         0.0000000        0.07555760
#>            r50
#>  1:         NA
#>  2:         NA
#>  3:         NA
#>  4:         NA
#>  5:         NA
#>  6: 0.04823151
#>  7: 0.06282723
#>  8: 0.07511737
#>  9: 0.08095238
#> 10: 0.10600707
#> 11: 0.12435233
#> 12: 0.13559322
#> 13: 0.14555256
#> 14: 0.14448669
#> 15:         NA
#> 16:         NA

Using plot() on your physiodata object will produce a plot of the cumulative germination curves:

plot(cent)

And using barplot() on your physiodata object will produce a plot of the final germination proportions and the median germination rate:

barplot(cent)

The example we have seen so far was about a single-group dataset, this is, a dataset about one experiment with one seedlot of one species. If your data has several groups (for example, several species, populations, or treatments of stratification), you must group your data with the groups argument of physiodata(). Also, if the column names in your dataset are the same as the default ones we defined, you don’t need to indicate them in the call to physiodata():

grass <- physiodata(grasses, 
                    x = "psi", # in this dataset the treatment is water potential
                    groups = "species") # this dataset has two different species
plot(grass)

Bradford’s hydrotime model

As long as you are working with a single-group dataset, you can fit Bradford’s hydrotime model (Gummerson 1986, Bradford 1990, Bewley et al. 2013) to your physiodata object using the function bradford():

anisantha <- physiodata(subset(grasses, species == "Anisantha rubens"), # select only one species
                        x = "psi") 
b <- bradford(anisantha$proportions) # bradford() uses the $proportions element within the physiodata object
b
#> Bradford's hydrotime model 
#>  Water potential levels in experiment: -1.6 -1 -0.8 -0.6 -0.4 -0.2 0 
#>  Theta - Hydrotime constant: 50.54 
#>  Psib50 - Base water potential (median): -1.59 
#>  Sigma of the base water potential: 0.4 
#>  R2: 0.94 
#> 

The result is a bradford object, which can be plotted using plot():

plot(b)

Garcia-Huidobro’s thermal time model

Likewise, you can fit Garcia-Huidobro’s thermal time model (Garcia-Huidobro et al. 1982, Gummerson 1986, Bewley et al. 2013) to your single-group physiodata object using the function huidobro():

malva <- physiodata(subset(centaury, population == "La Malva"), # select only one population
                    x = "temperature") 
h <- huidobro(malva$proportions) # huidobro() uses the $proportions element within the physiodata object
h
#> Garcia-Huidobro's thermal time model 
#>  Suboptimal temperature levels in experiment: 8.7 10.1 11.5 12.9 14.3 15.8 17.2 18.6 
#>  Supraoptimal temperature levels in experiment: 20 25 
#>  Tb - Base temperature: 5.2 
#>  To - Optimal temperature: 18.6 
#>  Tc - Ceiling temperature: 54.3 
#>  theta50 1 - Suboptimal thermal time (median): 33.52 
#>  Sigma of the suboptimal thermal time: 111.42 
#>  R2 of the suboptimal model: 0.39 
#>  theta50 2 - Supraoptimal thermal time (median): 449.8 
#>  Sigma of the supraoptimal thermal time: 460.87 
#>  R2 of the supraoptimal model: 0.23 
#> 

The result is a huidobro object, which can be plotted using plot():

plot(h)

The physiotime() function

As we said at the beginning, physiotime() is the basic function of seedr, and the only one you need to know for common usage of the package. physiotime() allows to centralize the operations described in the previous sections. It creates your physiodata object, divides it in the groups you indicate, and fits the model named in its method argument:

m <- physiotime(d = centaury, # the name of your data object
                t = "times", # the column with the scoring times
                x = "temperature", # the column with the experimental treatment,
                g = "germinated", # the column with the number of germinated seeds,
                pg = "germinable", # the column with the total number of viable seeds 
                groups = c("species", "population"), # the columns with grouping variables 
                method = "huidobro") # the model to fit, in this case a thermal time model
m
#> A list of physiological time germination models 
#>  calculated for the following 2 groups: 
#>  
#> species: Centaurium somedanum population: La Malva 
#> Garcia-Huidobro's thermal time model 
#>  Suboptimal temperature levels in experiment: 8.7 10.1 11.5 12.9 14.3 15.8 17.2 18.6 
#>  Supraoptimal temperature levels in experiment: 20 25 
#>  Tb - Base temperature: 5.2 
#>  To - Optimal temperature: 18.6 
#>  Tc - Ceiling temperature: 54.3 
#>  theta50 1 - Suboptimal thermal time (median): 33.52 
#>  Sigma of the suboptimal thermal time: 111.42 
#>  R2 of the suboptimal model: 0.39 
#>  theta50 2 - Supraoptimal thermal time (median): 449.8 
#>  Sigma of the supraoptimal thermal time: 460.87 
#>  R2 of the supraoptimal model: 0.23 
#>  
#> species: Centaurium somedanum population: Malconcecho 
#> Garcia-Huidobro's thermal time model 
#>  Suboptimal temperature levels in experiment: 7.3 8.7 10.1 11.5 12.9 14.3 15.8 17.2 
#>  Supraoptimal temperature levels in experiment: 18.6 20 25 
#>  Tb - Base temperature: 4 
#>  To - Optimal temperature: 17.5 
#>  Tc - Ceiling temperature: 32.4 
#>  theta50 1 - Suboptimal thermal time (median): 1.23 
#>  Sigma of the suboptimal thermal time: 155.22 
#>  R2 of the suboptimal model: 0.44 
#>  theta50 2 - Supraoptimal thermal time (median): 165.63 
#>  Sigma of the supraoptimal thermal time: 93.27 
#>  R2 of the supraoptimal model: 0.73 
#> 

physiotime() can accept further arguments to tune the fitting of the thermal time model, but generally you don’t need to touch these. You can read more about these extra arguments by running ?physiotime in R.

The result of physiotime() is a physiotime object. Using summary() on this object, you obtain a table of the fitted model parameters, which is easy to export for further uses:

summary(m)
#>                 species  population   method n.treatments.sub n.treatments.sup
#> 1: Centaurium somedanum    La Malva huidobro                8                2
#> 2: Centaurium somedanum Malconcecho huidobro                8                3
#>          Tb       To       Tc theta50.sub sigma.sub    R2.sub theta50.sup
#> 1: 5.179545 18.60158 54.26304   33.515734  111.4167 0.3923164    449.7998
#> 2: 3.998373 17.49968 32.44445    1.227839  155.2173 0.4381888    165.6308
#>    sigma.sup    R2.sup
#> 1: 460.86702 0.2342135
#> 2:  93.27052 0.7294838

Likewise, plot() will plot the fitted models:

plot(m)

Getting help

If you encounter a clear bug, please file an issue with a minimal reproducible example on GitHub.

References