A Holling Functional Response Model for Mapping QTLs Governing Interspecific Interactions

Genes play an important role in community ecology and evolution, but how to identify the genes that affect community dynamics at the whole genome level is very challenging. Here, we develop a Holling type II functional response model for mapping quantitative trait loci (QTLs) that govern interspecific interactions. The model, integrated with generalized Lotka-Volterra differential dynamic equations, shows a better capacity to reveal the dynamic complexity of inter-species interactions than classic competition models. By applying the new model to a published mapping data from a competition experiment of two microbial species, we identify a set of previously uncharacterized QTLs that are specifically responsible for microbial cooperation and competition. The model can not only characterize how these QTLs affect microbial interactions, but also address how change in ecological interactions activates the genetic effects of the QTLs. This model provides a quantitative means of predicting the genetic architecture that shapes the dynamic behavior of ecological communities.

Introduction

Understanding the internal workings of ecological communities is of fundamental importance to predict community dynamics and improve ecosystem services (Whitham et al., 2006; Vellend, 2020). Mounting evidence shows that genes play a pivotal role in shaping community structure, organization, and function (Whitham et al., 2006; Bailey et al., 2009; Hersch-Green et al., 2011; Miner et al., 2012; Crutsinger, 2016; Salazar et al., 2019; Wimp et al., 2019). For example, in response to predator-borne kairomones, some genes in Daphnia are activated, showing a higher-level expression than usually (Schwarzenberger et al., 2009; Miyakawa et al., 2010). To systematically characterize specific community genes, their number, chromosomal locations, and effect size, a powerful mapping approach has emerged through designing a community-ecological experiment using a mapping population (Jiang et al., 2018; Jiang et al., 2020).

Genetic mapping is a statistical approach widely used to map and identify genes, known as quantitative trait loci (QTLs), which control complex traits (Flint et al., 2005; Atwell et al., 2010; Consortium et al., 2015; Zeng et al., 2018). By integrating the mathematical aspect of trait formation, functional mapping has been developed to reveal the spatiotemporal pattern of the genetic architecture underlying phenotypic variation and evolution (Ma et al., 2002; Wu et al., 2005; Wu and Lin, 2006; Zhao et al., 2012; Li and Sillanpaa, 2015). The interpretable advantage of functional mapping has been leveraged to capture the biological rule governing how the components constituting a complex trait are interconnected, interdepended and interacted to mediate trait variation. A so-called system mapping approach has been assembled to map the genetic machineries underlying such component-component interconnections (Gai et al., 2011; Bo et al., 2014; Sun and Wu, 2015).

More recently, Jiang et al. (2018) integrated systems mapping and evolutionary game theory into a unified model to study the genetic control of species-species interactions (including cooperation and competition) in ecological communities. This cooperation-competition mapping (CoCoM) model differs from conventional systems mapping, in that the former needs to model genetic effects from two or multiple interacting species, whereas the latter only needs to consider one single genome of the species studied. By introducing the notion of evolutionarily stable strategy proposed by Smith and Price, CoCoM incorporates a system of nonlinear Lotka-Volterra (nLV) predator-prey equations (May, 1975) to partition the phenotypic value of each species in communities into its independent and dependent components. The independent component is one that occurs when this species is assumed to be in isolation, and the dependent component reflects the effect due to the influences of other species on this species. Thus, CoCoM is equipped with a capacity to reveal how much these two components contribute to phenotypic variation and how their relative contributions are controlled reciprocally by the genes of a pair of species. However, as the first model of its kind, CoCoM uses a simple from of nLV equations that may not adequately capture the chaotic complexity of ecological communities (May, 1975; Cushing, 1980; Arneodo et al., 1982; Panetta, 1996; Huisman and Weissing, 1999; D’Onofrio, 2002; Zeng et al., 2005).

The motivation of this study is to expand the utility of CoCo to a broader context that is allowed to be periodically oscillated. To model the perturbation of the community, we introduce a Holling-type functional response model to characterize impulsive perturbations of the nLV system. The Holling-type model has been widely used to study community behavior and dynamics at different scales (Liu and Chen, 2003; Zhang et al., 2005; Leeuwen et al., 2007). By analyzing a published dataset from an ecological experiment including monocultures and co-cultures of two bacterial species, the new model identifies previously unidentified genetic loci for microbial cooperation and competition. Comparing the difference of genetic architecture in socially isolated monocultures and socialized co-cultures gains new insight into the genetic mechanisms underlying microbial interactions. Computer simulation has been performed to validate the statistical properties of the new model.

Materials and Methods

Mapping Materials

We used a published experimental data (Jiang et al., 2018) to validate the utility of our model. The experiment, conducted with 45 strains from each of two bacterial species, Escherichia coli and Staphylococcus aureus, includes 45 strain-specific monocultures of each species and 45 interspecific co-cultures using independent strain pairs of two species. The abundance of each strain in monoculture and co-culture was measured repeatedly at 16 time points during growth process.

Holling Equations

Trait formation is a biological process that includes the increase of trait value with time, i.e., growth. Based on the biological rule governing growth, several growth equations, including logistic, Richards, and Gompertz, have been proposed to quantify the pattern of growth (Zwietering et al., 1990; West et al., 2001; Palacios et al., 2014). Functional mapping incorporates the growth equation to map QTLs for trait formation and development, showing increasing biological relevance and statistical power (Ma et al., 2002; Wu et al., 2005; Wu and Lin, 2006; Zhao et al., 2012; Li and Sillanpaa, 2015). However, these growth equations describe the growth of an organism without considering its biological surrounding where the other organisms would exert effects on its growth. By treating interactive organisms as a complex system, several nLV-based ordinary differential equations (ODEs) have been developed to model how different species interact with each other to determine the system (Fujikawa et al., 2014). nLV equations, widely applied to study prey-predator relationships in ecological communities (Hernández-Bermejo and Fairén, 1997; Terborgh and Estes, 2010), have been implemented to map cooperation and competition QTLs (Jiang et al., 2018; Jiang et al., 2020).

Beginning in the late 1950s, Holling conducted an experiment to investigate how a predator’s rate of prey capture is related to prey density (Holling, 1959a; Holling, 1959b). In the resulting series of seminal articles, Holling identified three general categories of functional response: Types I, II, and III (Holling, 1965; Morozov, 2010). Type I is the simplest: capture rate increases in direct proportion to prey density until it abruptly saturates. Type II is similar in that the rate of capture increases with increasing prey density, but in contrast to the linear increase of Type I, Type II approaches saturation gradually. Type III is similar to Type II except at a low prey density, where the rate of prey capture accelerates. Many experiments have proved that Holling-type functional response plays a key role in understanding predator-prey relationships (Ji et al., 2009; Tewa et al., 2013) and can fill the gap of nLV equations. Most studies view type II functional response as the basis for choosing an optimal foraging model because it is more likely to occur in communities where population dynamics is driven by predation (Ji et al., 2009).

Here, we integrate Holling type II functional response and nLV-based ODEs (HollinLV) to model the interaction and coordination mechanisms of two populations in communities. Assume that A and B are two species used to conduct an ecological experiment in which each species is grown in socially isolated monocultures and, also, both are grown in a socialized co-culture. Let E and S be the abundance of species A and species B, respectively, in a community. The HollinLV model is expressed as follows:

$\{\begin{array}{c}\frac{dE}{dt}=r_{e}E\left(1-\frac{E}{K_e}\right)+r_{e}E\left(\frac{{\alpha }_{E\leftarrow S}}{1+E}\right)S\\ \frac{dS}{dt}=r_{s}S\left(1-\frac{S}{K_s}\right)+r_{s}S\left(\frac{{\alpha }_{S\leftarrow E}}{1+S}\right)E\end{array}(1)$

where r_e and r_s represent the Malthusian growth rates of species A and species B, respectively; K_e and K_s are an intrinsic-carrying capacity of two different species; ${\alpha }_{E\leftarrow S}$, $and {\alpha }_{S\leftarrow E}$ represent the scalar parameter values that describe how one species affects the other through competition or cooperation in the co-culture. According to Eq. 1, the growth trajectories of each species is the summation of two parts:

$\{\begin{array}{c}\frac{dE}{dt}=r_{e}E\left(1-\frac{E}{K_e}\right)\\ \frac{dS}{dt}=r_{s}S\left(1-\frac{S}{K_s}\right)\end{array}(\mathrm{2A})$

and

$\{\begin{array}{c}\frac{dE}{dt}=r_{e}E\left(\frac{{\alpha }_{E\leftarrow S}}{1+E}\right)S\\ \frac{dS}{dt}=r_{s}S\left(\frac{{\alpha }_{S\leftarrow E}}{1+S}\right)E\end{array}(\mathrm{2B})$

Following Jiang et al. (2018), we define the first part as the independent growth of a species that occurs when this specie is assumed to grow in isolation, which corresponds to the growth of this species in monoculture, and the second part as the dependent growth of a species, determined by the interaction of its co-existing species through a certain mechanism. The degree of this dependence is described by interaction scalar parameters ${\alpha }_{E\leftarrow S}\mathrm{and}{\alpha }_{S\leftarrow E}$.

According to population ecology theory, there are mainly three types of interaction among species (Svenning et al., 2014): 1) neutral interaction, i.e., there is no interaction between species and they are independent of each other; 2) positive interaction, including commensalism and mutualism. If only one side is favorable and there is no influence on the other side, it is called partial benefit symbiosis; if both sides are favorable, it is called mutualism; and 3) negative interaction, including antagonism, predation and amensalism. Antagonism is opposite to mutualism; predator preys on other species, which is the traditional relationship between preying and being preyed on; amensalism means that the existence of one species has an inhibitory effect on the other, but the second no effect on the first.

By estimating the interaction and scale parameters, we can quantify the degree of each possible interaction. While the interaction parameters determine whether and in which direction the dependent growth occurs, scale parameters, ${\alpha }_{E\leftarrow S}\mathrm{and}{\alpha }_{S\leftarrow E}$, determine how the interaction strengthens or weakens with growth development. A positive or negative, ${\alpha }_{E\leftarrow S}\mathrm{or}{\alpha }_{S\leftarrow E}$ value indicates that this species is benefitted or harmed by another species. If the value of ${\alpha }_{E\leftarrow S}\mathrm{or}{\alpha }_{S\leftarrow E}$ is zero, this means that the interaction does not depend on another species. If they are also zero, which means that the two species are not affected by one another.

A Mapping Framework

In the co-culture environment, both HollinLV and CoCoM models consider the interactions between species. Here, we incorporate HollinLV into systems mapping, an approach of mapping complex traits by treating complex traits as a dynamic system (Gai et al., 2011; Bo et al., 2015; Sun and Wu, 2015). Consider two mapping populations each from a different species, A or B. Each population has n genome-wide genotyped members, grown individually in monoculture and in co-culture in interspecific pairs. Let $E_i=\left(E_i\left(1\right),\cdots ,E\left(T\right)\right)$ and $S_i=\left(S_i\left(1\right),\cdots ,S_i\left(T\right)\right)\left(i=1,\cdots ,n\right)$ denote the abundance of n pairs from species A and B measured at time T in co-culture. Consider a genetic locus with two genotypes A and a for species A and two genotypes B and b for species B. Pairing species A and B generates four interspecific genotype combinations AB, Ab, aB, and ab, with observations denoted as n_AB, n_Ab, n_aB, and n_ab, respectively. We formulate the likelihood function for the abundance data of interspecific pairs at the locus as

$L={\displaystyle \prod }_{i=1}^{n_{AB}}{f}_{AB }\left(E_i,S_i;{\theta }_{AB}\right)\times {\displaystyle \prod }_{i=1}^{n_{Ab}}{f}_{Ab }\left(E_i,S_i;{\theta }_{Ab}\right)\times {\displaystyle \prod }_{i=1}^{n_{aB}}{f}_{aB }\left(E_i,S_i;{\theta }_{aB}\right)\times {\displaystyle \prod }_{i=1}^{n_{ab}}{f}_{ab }\left(E_i,S_i;{\theta }_{ab}\right)(3)$

where θ is the HollinLV parameter ($r_e,K_e, {\alpha }_{E\leftarrow S}, r_s, K_s, {\alpha }_{S\leftarrow E}$) describing abundance, $f_j\left(E_i,S_i;{\theta }_j\right),\mathrm{j}\in \left(AB,Ab,ab,aB\right)$ is a bivariate longitudinal normal distribution specified by genotype-specific mean vectors,

${\mu }_j=\left({\mu }_{j1};{\mu }_{j2}\right)=\left({\mu }_{j1}\left(1\right),\cdots ,{\mu }_{j1}\left(T\right);{\mu }_{j2}\left(1\right),\cdots ,{\mu }_{j2}\left(T\right)\right)(4)$

and covariance matrix,

$\text{Σ}=\left(\begin{array}{cc}{\text{Σ}}_{11}& {\text{Σ}}_{12}\\ {\text{Σ}}_{21}& {\text{Σ}}_{22}\end{array}\right),(5)$

with variance matrices of abundance over time on the middle diagonal and covariance matrices of abundance between two species on the non-diagonal.

We model the vector structure of Eq. 4 by HollinLV Eq. 1 and the structure of covariance matrix Eq. 5 by a first-order structured antedependence [SAD(1)] statistical model (Zimmerman and Nunez-Anton, 1997; Zhao et al., 2005a; Zhao et al., 2005b). We implement a hybrid of the Nelder-Mead simplex algorithm and the fourth-order Runge-Kutta algorithm to solve the likelihood Eq. 3. The maximum likelihood estimates (MLEs) of ODE parameters and SAD(1) parameters can be obtained.

Hypothesis Testing

The existence of significant QTLs associated with interspecific interactions can be tested by a log-likelihood ratio (LR) approach. To do so, we formulate two hypotheses expressed as

$\begin{array}{l}{H}_0: {\Theta }_j\equiv \Theta ,\\ H_1: {\Theta }_j\ne \Theta ,\end{array}(6)$

$\mathrm{for}\mathrm{j}\in \left(AB,Ab,ab,aB\right)$, under which the likelihood values L₀ and L₁ are calculated, respectively. The test statistic LR is calculated as

$\mathrm{LR}=2\left(\log L_1-\log L_0\right)(7)$

which is compared with the critical threshold. Note that the null hypothesis states that these ODE parameters are genotype invariant. The rejection domain of the likelihood ratio statistics is, $W=\left\{\mathrm{LR}\ge c\right\}$, where c satisfies

$P_{\theta }\left(\mathrm{LR}\ge c\right)\le \alpha (8)$

When the value of the test statistic LR is in the rejection domain, the null hypothesis is rejected, which suggests that there are differences in the abundance from different interspecific combination genotypes. We can also test whether a significant QTL affects the independent growth by formulating

$\begin{array}{l}{H}_0:\left(r_{ej},k_{ej},r_{sj},k_{sj}\right)\equiv \left(r_e,k_e,r_s,k_s\right)\\ H_1:\left(r_{ej},k_{ej},r_{sj},k_{sj}\right)\ne \left(r_e,k_e,r_s,k_s\right)\end{array}(9)$

or the interaction growth by formulating

$\begin{array}{l}{H}_0:\left(r_{ej},{\alpha }_{E\leftarrow Sj};r_s,{\alpha }_{S\leftarrow Ej}\right) \equiv \left(r_e,{\alpha }_{E\leftarrow S},r_s,{\alpha }_{S\leftarrow E}\right),\\ H_1: \left(r_{ej},{\alpha }_{E\leftarrow Sj};r_s,{\alpha }_{S\leftarrow Ej}\right)\ne \left(r_e,{\alpha }_{E\leftarrow S},r_s,{\alpha }_{S\leftarrow E}\right),\end{array}(10)$

The LRs for the above two pairs of hypotheses are calculated and compared with the critical thresholds, respectively.

According to Jiang et al. (2018), we partition the overall genotype values of significant SNP pairs into direct effects, indirect effects, and genome-genome epistatic effects. The direct effect describe how a QTL directly affects the abundance dynamics of each species from its own genome. The indirect effect describes how a QTL indirectly affects the abundance dynamics of each species from the genome of its interactive genome. The epistatic effect characterizes how the interaction of alleles from different species affects the abundance dynamics of each species. Traditional genetic mapping can only estimate and test the indirect genetic effects, whereas our model can estimate and test all these three effects, which thus can gain new insight into the genetic architecture of interspecific interactions in ecological communities.

Results

Fitting Growth Curves

We used Gompertz, Logistic and Richards growth equations, Jiang et al. (2018) CoCoM, and our HollinLV to fit the mean growth curves of microbial abundance for 45 strains from each species in co-culture. We compared the performance of these models by calculating AIC (Akaike information criterion), BIC (Schwarz criterion), and HQ (Hannan Quinn criterion) (Table 1). As shown, CoCoM and HollinLV perform better than three classical growth equations, suggesting that it is crucial to consider the impact of interactive species on the growth of a given species in a community. Our HollinLV shows a better performance than CoCoM. In order to show the fitting effects specifically, the actual values and fitting curves of 45 pairs of co-culture samples are plotted, and the goodness of fit ($R^2$) are calculated in Supplementary Figures S1, S2, indicating the samples achieve good fitting effects.

TABLE 1

TABLE 1. The estimated parameters of various growth equations, the goodness of fit $R^2$, adjusted $R^2$, and their evaluation information based on akaike information criterion (AIC), bayesian information criterion (BIC), and hannan-quinn information criterion (HQIC). Gompertz, Logistic, Richards, DDHR, and CoCoM models were used to fit the average abundance of E. coli and S. aureus under co-culture.

Figure 1 illustrates the abundance trajectories of the same species in co-culture and monoculture. We find that the growth of two species in co-culture is affected by each other to a certain extent. The most significant change of growth trajectories displays in the timing of maximum relative growth rate; microbes from both species tend to enter fast-growing stages in co-culture than in monoculture. However, the degree to which S. aureus growth is affected by E. coli is greater than that to which E. coli growth is affected by S. aureus. This suggests that, as compared to E. coli, the growth of S. aureus is more sensitive to the biotic environment. All these biotically induced changes imply the existence of specific genes that are activated by species coexistence.

FIGURE 1

FIGURE 1. Growth curve of microbial abundance under monoculture and co-culture. Microbial abundance of individual strains from each species was observed in monoculture and co-culture during the first 36 h after culture. The upper: the abundance curve of E. coli under monoculture and co-culture, respectively. The thick blue line is the average growth curve of E. coli. The lower: the abundance curve of S. aureus under monoculture and co-culture, respectively. The thick red line is the average growth curve of S. aureus.

How QTLs Are Activated by Species Coexistence

We apply HollinLV to map QTLs for the microbial abundance of each species in co-culture, in a comparison with those detected in monoculture. To map growth trajectories in co-culture, we implement two ODEs of Eq. 1, each for a different species, into a mapping framework, whereas the corresponding mapping of growth in monoculture is based on the implementation of two ODEs in Eq. 2A.

Figure 2A shows the genomic distribution of QTLs detected for each species in different cultures. We identify more growth QTLs for S. aureus than E. coli in both monoculture and co-culture, but for both species, a number of new QTLs are activated from monoculture to co-culture (Figure 2B). In monoculture, there are a total of 81 QTLs that influence the growth trajectory of E. coli, but this number increases to 132 in co-culture, of which 32 are the common QTLs for both types of cultures. In other words, 49 QTLs are specifically expressed in monoculture, whereas 100 only function in co-culture. In S. aureus, 123 QTLs are monoculture-specific, 363 are co-culture-specific, and 23 are common to both cultures.

FIGURE 2

FIGURE 2. Identification of significant QTLs and comparative analysis for E. coli and S. aureus under different culture. (A) Manhattan plots in co-culture (the upper part of A) and monoculture (the lower part of A) for E. coli and S. aureus, respectively. From Manhattan plots A, significant QTLs are detected through significance tests. The horizontal lines are the genome-wide critical threshold at 0.01significance level. (B) Comparison of significant QTLs for microbial growth under different culture conditions. The blue color represents co-culture condition and the red color represents monoculture condition.

We estimate the genetic effect curves of each QTL in monoculture and co-culture for both species. Figure 3 depicts how three representative common QTLs from each species change their genetic effects during growth. Overall, the temporal effects of QTLs cyclically change with time; i.e., a QTL increases its effect after culture and then monotonically decrease its effect at a certain time point. For the culture-common QTLs, the pattern of time-varying genetic effects depends on culture type. It is interesting to note that the same QTL displays a peak of its effect earlier in co-culture than in monoculture. For example, QTL E4304609 reaches a maximum effect when E. coli was cultured for 5 h in co-culture, but the timing of its maximum effect occurs at time 10 h after culture in monoculture. For QTL S2217304 from S. aureus, its maximum effect occurs at 3 and 8 h in co-culture and monoculture, respectively. Furthermore, this QTL changes the direction of its genetic effect from monoculture to co-culture.

FIGURE 3

FIGURE 3. The genetic effect curves regulated by common QTLs which are significant under both single culture and co-culture conditions. The left: the effect curves of E. coli of E4489505, E4304609, and E509544, respectively. The right: the effect curves of S. aureus of S2217304, S293354, and S1308, respectively. The red lines represent the genetic effect under co-culture, and the blue lines represent the genetic effect under monoculture.

Mapping the Two-Dimensional Genetic Architecture of Microbial Growth in Co-culture

Unlike a case of monoculture in which the phenotype of a species is only controlled by its own genes, the growth of one species may be controlled epistatically by its own genome and the genome of its co-existing species in co-culture. We find that in co-culture, E. coli and S. aureus, tend to be antagonistic; i.e., one species grows at a cost of the other species, with the strength of such antagonism being quantified by the HollinLV model Eq. 1 (Figure 4A). Our mapping model can characterize significant interspecifically epistatic SNP pairs that affect a species’ growth trajectory in co-culture. By pairwise scanning loci, each from a different species, throughout the E. coli and S. aureus genomes, we plot a two-dimensional Manhattan plot that demonstrates 2,245 significant interspecific SNP pairs comprising of 249 SNPs from E. coli and 182 SNPs from S. aureus (Figure 4B). We determined the critical threshold at the significance level of 10⁻⁶ after Bonferroni correction. The functional annotation of genes from NCBI’s database shows that almost all QTLs detected residue within candidate genes of known biological functions. For example, for E. coli, E19056 resides in gene nhaR activating the distal promoter, osmCp1, of expression of osmC. It can stimulate osmCp1 in response to an osmotic signal (Sturny et al., 2003; Katarzyna et al., 2010). E49080 resides in gene kefC, which is a potassium transport system regulated by glutathione metabolites (Miller et al., 2000). E62059 is relevant to the rapA controlling an additional mechanism, which involved in wild-type biofilm resistance (Lynch et al., 2007). For S. aureus, S188004 exerts pronounced genetic interactions with many SNPs distributed over the E. coli genome, it resides in gene ggt associating with the gammaglutamyltranspeptidase facilitating glutathione utilization (Hussain and Roslinah, 2008).

FIGURE 4

FIGURE 4. The two-dimensional genetic analysis results of interspecies interactions between E. coli and S. aureus in community growth. (A) The overall microbial abundance (blue line) of E. coli and S. aureus in co-culture, and the underlying independent (red line), and interactive (green line) growth components. (B) Two-dimensional Manhattan plots of significant QTL pairs on two microbial genomes. The genome-wide significance threshold is determined by ${10}^{-6}$ significance level after Bonferroni correction. (C) Microbial growth curves of four genotype combinations T/T, T/C, C/T, and C/C formed by QTL E4320079 from E. coli and S188004 from S. aureus.

At each SNP pair, there are four possible genotype combinations between two species. To show how a SNP from different genomes interacts with each other to affect microbial abundance, we decompose the net growth curves of each genotype combination into independent curves and dependent curves for each species. Figure 4C illustrates such an example from SNPs E4320079 with two genotypes T and C residing within the region of phnI gene and S188004 with two genotypes T and C residing within the region of ggt gene, producing combinations TT, TC, CT, and CC. We find that net, independent, and dependent growth curves for each species vary among four genotype combinations; for example, CT displays less independent growth in E. coli than the other combinations, whereas there is less independent growth in S. aureus for TT than the other combinations. For all combinations, E. coli and S. aureus are antagonistic, but with the strength of antagonism depending on genotype combination, suggesting that SNP pair E4320079 and S188004 plays a role in mediating the microbial growth of two coexisting species.

We implemented Jiang et al. (2018) model to partition the genotypic values of each significant interspecific SNP pair into direct genetic effects, indirect genetic effects, and genome-genome epistatic effects on the growth trajectory of each species. Based on the estimates of these effect values, we calculate the genetic variance curves due to each effect. It is interesting to see that the growth of each species in co-culture is not only affected directly by its own genes, but also, to a similar extent, indirectly by genes from its co-existing conspecific and epistatically by cross-genome interactions of the gene between two species (Figure 5).

FIGURE 5

FIGURE 5. The genotypic value of each genotype combination is partitioned into its direct (red), indirect (blue), and genome–genome (G–G) epistatic effects (green) on the growth, with the time-varying genetic variance explained by each of these effects of E. coli (A) and S. aureus (B) in co-culture. The thick lines on the variance curves represent the average variance of direct, indirect, and epistatic effects, respectively.

Computer Simulation

To investigate the utility and the statistical properties of HollinLV, we performed simulation studies by mimicking the real example as described above. The data were simulated by assuming that two species were reared in monoculture and co-culture. The phenotype was determined by a set of QTLs among 1,000 simulated markers, plus a residual error following a bivariate normal distribution. We simulated n = 45, 100, and 200 interspecific pairs in co-culture under heritabilities of H²fn2 = 0, 0.05, and 0.10. The size of heritability was used to adjust the magnitude of innovative variance.

Four genotypic curves at a joint locus of two genomes can be reasonably well estimated even using a small sample size under a modest heritability level (H² = 0.05) (Supplementary Figure S3A; Supplementary Table S1). As expected, the accuracy and precision of curve estimation can be increased by increasing sample size and heritability. In general, the level at which estimation precision improves from H² = 0.05 to 0.10 under n = 0.45 is similar to the level of improvement from n = 45 to 100 under H2 = 0.05 (Supplementary Figure S3B; Supplementary Table S1). Because Jang et al. (2018) experiment was well controlled in a uniform environment, producing minimum noises, we can expect that abundance data contains so adequately large a heritability that small sample size n = 45 can provide reasonably precise estimation of genotypic curves by our HollinLV model.

We performed computer studies to investigate the empirical power of HollinLV for QTL detection under different samples and heritabilities (Table 2). When both heritability and sample size are small (H² = 0.05 and n = 45), the model’s power is quite low (0.37), but it increases dramatically when H² increases to 0.10 (0.76) or when n increases to 100 (0.85). As discussed above, Jiang et al. (2018) experiment might have increasing heritability because of well-controlled environment. However, the results from such a small sample size should be interpreted with caution. When sample size increases to 200, the power of QTL detection can increase to >0.96 even with a small heritability. Therefore, in practice, a sample size of n = 200 is recommended for our HollinLV model to detect interaction QTLs. In general, the model has a low false positive rate, especially when sample size is 100 or higher (Table 2).

TABLE 2

TABLE 2. The positioning accuracy and false positive probability of significant QTLs identified by DDHR. The experiments simulated the co-culturing 45 pairs (the same as the real example), 100, pairs and 200 pairs under heritability levels of 0, 0.05 and 0.10. The positioning accuracy and false positive probability were evaluated by computer simulation.

Discussion

Community genetics has emerged as a subdiscipline of genetics that combines community ecology and ecological genetics to gain insight into the genetic mechanisms underlying phenotypic diversity and evolution within and between species (Whitham et al., 2006; Hersch-Green et al., 2011). In this article, we develop and implement a computational model to address three fundamental questions in community genetics: 1) how a given species genetically adapts to its coexisting conspecifics, 2) which genetic machineries mediate a species’ phenotypes expressed in ecological communities, and 3) what is the genetic architecture underlying interspecific interactions. This model is the integration of systems mapping and community ecology through a system of HollinLV equations. Systems mapping views and maps a complex trait as a system whose interconnected components act together to mediate the trait, following a similar principle of community ecology, i.e., different species interact with each other to shape community dynamics.

There has been a rich body of literature on studying the genetic control of the phenotypic response of an organism to abiotic environment (Wang et al., 2013; Jiang et al., 2018; Diouf et al., 2020), but knowledge about how genes act differently in response to biotic environment is very limited. It is impossible to answer this question without the integration of quantitative genetic models and ecological data. Our HollinLV model has made this integration possible. An ecological experiment generates data from socially isolated environments (monoculture) and socialized environment (co-culture). HollinLV can test the change of genetic control in terms of the number and action of gens involved from monoculture to co-culture, providing insight into how a given species genetically responds to the co-existence of other species. By analyzing a real dataset, HollinLV identifies a set of specific QTLs that are activated by the existence of other species. In the co-culture of two bacterial species E. coli and S. aureus, a high percentage of QTLs are detected to be different from those expressed in monoculture. Such a percentage is highly species dependent; e.g., 75.76% for E. coli and 91.30% for S. aureus. Many more new QTLs activated for S. aureus in response to E. coli than for E. coli in response to S. aureus may be due to a higher degree of inhibition for S. aureus by E. coli than for E. coli by S. aureus. These discoveries could have immediate implications; for example, if altering the expression of genes of S. aureus that are related to with interactions with E. coli, the abundance of S. aureus can be controlled. A similar genetic manipulation can be made to control the abundance of E. coli.

Classic quantitative genetic theory suggests that the phenotype of an organism is controlled by its genes and surrounding environment (Lynch and Walsh, 1997). However, this theory does not interpret how a phenotypic trait is genetically controlled in ecological communities in which intra- and interspecific interactions are also a driver of phenotypic variation (Whitham et al., 2006; Bailey et al., 2009; Hersch-Green et al., 2011; Miner et al., 2012; Crutsinger, 2016; Wimp et al., 2019). By integrating community ecology theory, our HollinLV model equips quantitative genetic theory with new power to more comprehensively decipher the genetic landscape of biological diversity and community behavior. HollinLV models the phenotypic variation of a species through a joint genetic action of different genomes each from a co-existing species. This modeling framework allows us to dissect genetic control into direct genetic effects (due to genes from a given species’ own genome), indirect genetic effects resulting from genes from the genomes of other species coexisting with this species, and genome-genome epistatic effects (arising from the interaction between genes from co-existing species). Classic quantitative genetic approaches can only estimate the direct genetic effects of a complex trait, whereas our HollinLV model leverages these approaches to characterize the impacts of indirect and trans-genome interaction effects. By analyzing Jiang et al. (2018) real data, we find that the phenotype of a bacterial species is determined not only by its own genes, but also by the genes of its coexisting conspecific and, remarkably, by the interaction of genes from two species. In general, E. coli’s genes are found to exert more impactful indirect effects on the abundance of S. aureus than those of S. aureus’ genes on E. coli abundance (Figure 5), showing that E. coli plays a dominant role in shaping the coexistence of two species. Furthermore, trans-genome interactions are an important component of the genetic architecture underlying each species’ abundance. The precise dissection of overall genetic control over community phenotypes by HollinLV can greatly enhance our understanding of community behavior, dynamics and evolution.

HollinLV divides the net (observed) phenotype of a species into two components, the independent component arising from this species’ intrinsic capacity and the dependent components due to the influence of other species on it, and quantifies each of these two components. The sign and size of the dependent components for a pair of species characterize all possible patterns of ecological interactions from mutualism to commensalism to amensalism to predation to antagonism (Diouf et al., 2020). HollinLV allows us to characterize the genetic control mechanisms of each type of interaction and compare how the same genes differently govern these interactions. In Jiang et al. (2018) data, two bacterial species were found to be antagonistic to each other, although the strength of antagonism is larger for E. coli towards S. aureus than for S. aureus towards E. coli. Specific QTLs from each species have been identified to mediate their antagonism in co-culture. With these antagonism-related QTLs and the corresponding genes, microbial geneticists can destroy this antagonistic relationship or change it to other types of interactions with aid of gene editing techniques.

HollinLV is an ecology-oriented mapping model that can potentially gain new insight into community genetics. Its application to a microbial cultural experiment shows some promise to better understand the genetic control of interspecific interactions between two bacterial species. However, results from our analysis should be interpreted with caution because the data contains a small sample size. The future study should increase sample size to 200 interspecific pairs by which parameter estimation and detection power can reach a satisfactory level. Also, the model derivation is based on a pair of species in co-culture. This obviously is not adequate to be applied to practical ecological communities where many species are coexisting and interact with each other to form a complex ecological network. Holling-typed response models only characterize an aspect of community systems, and more comprehensive models are required to better capture multifaced features of communities. Regardless, our HollinLV provides a start point for integrating quantitative genetic theory and community ecology to disentangle the complexity of community genetics.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: The raw sequence data were deposited in the NCBI short reads archive under accession number SRP074089 and SRP074912.

Author Contributions

X-YZ contributed to methodology and supervision, and wrote the article, HG performed analysis and visualization, QF contributed to analysis and project administration, XZ supported the scientific research funding for the publication of this study, LJ made data curation, investigation, and methodology, RW conceived, and designed the analysis.

Funding

The work is supported by the Fundamental Research Funds for the Central Universities (No. BLX201912) (XZ) and the National Natural Science Foundation of China (Grant No. 11501032) (X-YZ). The work is partially supported by Japan Society for the Promotion of Science (Grant No. 19K03613) (QF).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.766372/full#supplementary-material

References

Arneodo, A., Coullet, P., Peyraud, J., and Tresser, C. (1982). Strange Attractors in Volterra Equations for Species in Competition. J. Math. Biol. 14 (2), 153–157. doi:10.1007/BF01832841

CrossRef Full Text | Google Scholar

Atwell, S., Huang, Y. S., Vilhjálmsson, B. J., Willems, G., Horton, M., Li, Y., et al. (2010). Genome-wide Association Study of 107 Phenotypes in Arabidopsis thaliana Inbred Lines. Nature 465 (7298), 627–631. doi:10.1038/nature08800

PubMed Abstract | CrossRef Full Text | Google Scholar

Bailey, J. K., Schweitzer, J. A., Úbeda, F., Koricheva, J., LeRoy, C. J., Madritch, M. D., et al. (2009). From Genes to Ecosystems: a Synthesis of the Effects of Plant Genetic Factors across Levels of Organization. Phil. Trans. R. Soc. B 364 (1523), 1607–1616. doi:10.1098/rstb.2008.0336

PubMed Abstract | CrossRef Full Text | Google Scholar

Bo, W., Fu, G., Wang, Z., Xu, F., Shen, Y., Xu, J., et al. (2014). Systems Mapping: How to Map Genes for Biomass Allocation toward an Ideotype. Brief. Bioinform. 15 (4), 660–669. doi:10.1093/bib/bbs089

PubMed Abstract | CrossRef Full Text | Google Scholar

Consortium, U., Walter, K., Min, J. L., Huang, J., Crooks, L., Memari, Y., et al. (2015). The UK10K Project Identifies Rare Variants in Health and Disease. Nature 526 (7571), 82–90. doi:10.1038/nature14962

PubMed Abstract | CrossRef Full Text | Google Scholar

Crutsinger, G. M. (2016). A Community Genetics Perspective: Opportunities for the Coming Decade. New Phytol. 210 (1), 65–70. doi:10.1111/nph.13537

PubMed Abstract | CrossRef Full Text | Google Scholar

Cushing, J. M. (1980). Two Species Competition in a Periodic Environment. J. Math. Biol. 10 (4), 385–400. doi:10.1007/BF00276097

CrossRef Full Text | Google Scholar

Diouf, I., Derivot, L., Koussevitzky, S., Carretero, Y., Bitton, F., Moreau, L., et al. (2020). Genetic Basis of Phenotypic Plasticity and Genotype × Environment Interactions in a Multi-Parental Tomato Population. J. Exp. Bot. 71, 5365–5376. doi:10.1093/jxb/eraa265

CrossRef Full Text | Google Scholar

D’Onofrio, A. (2002). Stability Properties of Pulse Vaccination Strategy in the SIR Epidemic Model. Bull. Math. Biol. 60, 1–26.

Google Scholar

Flint, J., Valdar, W., Shifman, S., and Mott, R. (2005). Strategies for Mapping and Cloning Quantitative Trait Genes in Rodents. Nat. Rev. Genet. 6 (4), 271–286. doi:10.1038/nrg1576

PubMed Abstract | CrossRef Full Text | Google Scholar

Fujikawa, H., Munakata, K., Sakha, M. Z., and Sakha, M. Z. (2014). Development of a Competition Model for Microbial Growth in Mixed Culture. Biocontrol Sci. 19 (2), 61–71. doi:10.4265/bio.19.61

PubMed Abstract | CrossRef Full Text | Google Scholar

Hernández-Bermejo, B., and Fairén, V. (1997). Lotka-Volterra Representation of General Nonlinear Systems. Math. Biosciences 140 (1), 1–32. doi:10.1016/S0025-5564(96)00131-9

CrossRef Full Text | Google Scholar

Hersch-Green, E. I., Turley, N. E., and Johnson, M. T. J. (2011). Community Genetics: what Have We Accomplished and where Should We Be Going? Phil. Trans. R. Soc. B 366 (1569), 1453–1460. doi:10.1098/rstb.2010.0331

PubMed Abstract | CrossRef Full Text | Google Scholar

Holling, C. S. (1959b). Some Characteristics of Simple Types of Predation and Parasitism. Can. Entomol. 91 (7), 385–398. doi:10.4039/Ent91385-7

CrossRef Full Text | Google Scholar

Holling, C. S. (1959a). The Components of Predation as Revealed by a Study of Small-Mammal Predation of the European Pine Sawfly. Can. Entomol. 91 (5), 293–320. doi:10.4039/Ent91293-5

CrossRef Full Text | Google Scholar

Holling, C. S. (1965). The Functional Response of Predators to Prey Density and its Role in Mimicry and Population Regulation. Mem. Entomol. Soc. Can. 97 (45), 5–60. doi:10.4039/entm9745fv

CrossRef Full Text | Google Scholar

Huisman, J., and Weissing, F. J. (1999). Biodiversity of Plankton by Species Oscillations and Chaos. Nature 402 (6760), 407–410. doi:10.1038/46540

CrossRef Full Text | Google Scholar

Hussain, R. M. (2008). Analysis of the Role of Glutathione and Stress Resistance in Staphylococcus aureus. Diss. [Sheffield]: University of Sheffield.

Google Scholar

Ji, C., Jiang, D., and Shi, N. (2009). Analysis of a Predator-Prey Model with Modified Leslie-Gower and Holling-type II Schemes with Stochastic Perturbation. J. Math. Anal. Appl. 359 (2), 482–498. doi:10.1016/j.jmaa.2009.05.039

CrossRef Full Text | Google Scholar

Jiang, L., He, X., Jin, Y., Ye, M., Sang, M., Chen, N., et al. (2018). A Mapping Framework of Competition-Cooperation QTLs that Drive Community Dynamics. Nat. Commun. 9 (1), 3010. doi:10.1038/s41467-018-05416-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, L., Liu, X., He, X., Jin, Y., Cao, Y., Zhan, X., et al. (2020). A Behavioral Model for Mapping the Genetic Architecture of Gut-Microbiota Networks. Gut Microbes 13, 1–15. doi:10.1080/19490976.2020.1820847

CrossRef Full Text | Google Scholar

Kennedy, M. S., and Estes, J. A. (2012). Trophic Cascades: Predators, Prey and the Changing Dynamics of Nature. Austral Ecol. 37 (2), e9–e10. doi:10.1111/j.1442-9993.2012.02361.x

CrossRef Full Text | Google Scholar

Leeuwen, E. V., Jansen, V. A. A., and Bright, P. W. (2007). How Population Dynamics Shape the Functional Response in a One-Predator-Two-Prey System. Ecology 88 (6), 1571–1581. doi:10.1890/06-1335

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., and Sillanpää, M. J. (2015). Dynamic Quantitative Trait Locus Analysis of Plant Phenomic Data. Trends Plant Sci. 20, 822–833. doi:10.1016/j.tplants.2015.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X., and Chen, L. (2003). Complex Dynamics of Holling Type II Lotka-Volterra Predator-Prey System with Impulsive Perturbations on the Predator☆. Chaos, Solitons & Fractals 16 (2), 311–320. doi:10.1016/S0960-0779(02)00408-3

CrossRef Full Text | Google Scholar

Lynch, M., and Walsh, B. (1997). Genetics and Analysis of Quantitative Traits. Sunderland, MA: Sinauer Associates.

Google Scholar

Lynch, S. V., Dixon, L., Benoit, M. R., Brodie, E. L., Keyhan, M., Hu, P., et al. (2007). Role of the Rapa Gene in Controlling Antibiotic Resistance of escherichia Coli Biofilms. Antimicrob. Agents Chemother. 51 (10), 3650–3658. doi:10.1128/aac.00601-07

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, C.-X., Casella, G., and Wu, R. (2002). Functional Mapping of Quantitative Trait Loci Underlying the Character Process: a Theoretical Framework. Genetics 161 (4), 1751–1762. doi:10.1093/genetics/161.4.1751

PubMed Abstract | CrossRef Full Text | Google Scholar

May, R. H. (1975). Stability and Complexity in Model Ecosystems. 2nd edition. Princeton: Princeton University Press.

Google Scholar

Miller, S., Ness, L. S., Wood, C. M., Fox, B. C., and Booth, I. R. (2000). Identification of an Ancillary Protein, Yabf, Required for Activity of the Kefc Glutathione-Gated Potassium Efflux System in escherichia Coli. J. Bacteriol. 182 (22), 6536–6540. doi:10.1128/JB.182.22.6536-6540.2000

CrossRef Full Text | Google Scholar

Miner, B. E., De Meester, L., Pfrender, M. E., Lampert, W., and Hairston, N. G. (2012). Linking Genes to Communities and Ecosystems: Daphnia as an Ecogenomic Model. Proc. R. Soc. B. 279 (1735), 1873–1882. doi:10.1098/rspb.2011.2404

PubMed Abstract | CrossRef Full Text | Google Scholar

Miyakawa, H., Imai, M., Sugimoto, N., Ishikawa, Y., Ishikawa, A., Ishigaki, H., et al. (2010). Gene Up-Regulation in Response to Predator Kairomones in the Water Flea, daphnia Pulex. BMC Dev. Biol. 10 (1), 45. doi:10.1186/1471-213X-10-45

PubMed Abstract | CrossRef Full Text | Google Scholar

Morozov, A. Y. (2010). Emergence of Holling Type III Zooplankton Functional Response: Bringing Together Field Evidence and Mathematical Modelling. J. Theor. Biol. 265 (1), 45–54. doi:10.1016/j.jtbi.2010.04.016

CrossRef Full Text | Google Scholar

Palacios, A. P., Marín, J. M., Quinto, E. J., and Wiper, M. P. (2014). Bayesian Modeling of Bacterial Growth for Multiple Populations. Ann. Appl. Stat. 8 (3), 1516–1537. doi:10.1214/14-AOAS720

CrossRef Full Text | Google Scholar

Panetta, J. C. (1996). A Mathematical Model of Periodically Pulsed Chemotherapy: Tumor Recurrence and Metastasis in a Competitive Environment. Bltn Mathcal Biol. 58 (3), 425–447. doi:10.1007/BF02460591

CrossRef Full Text | Google Scholar

Potrykus, K., Murphy, H., Chen, X., Epstein, J. A., and Cashel, M. (2010). Imprecise Transcription Termination within Escherichia coli greA Leader Gives Rise to an Array of Short Transcripts, GraL. Gral. Nucleic Acids Res. 38 (5), 1636–1651. doi:10.1093/nar/gkp1150

PubMed Abstract | CrossRef Full Text | Google Scholar

Salazar, G., Paoli, L., Alberti, A., Huerta-Cepas, J., Ruscheweyh, H.-J., Cuenca, M., et al. (2019). Gene Expression Changes and Community Turnover Differentially Shape the Global Ocean Metatranscriptome. Cell 179 (5), 1068–1083. doi:10.1016/j.cell.2019.10.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwarzenberger, A., Courts, C., and von Elert, E. (2009). Target Gene Approaches: Gene Expression in Daphnia magna Exposed to Predator-Borne Kairomones or to Microcystin-Producing and Microcystin-free Microcystis Aeruginosa. BMC Genomics 10 (1), 527. doi:10.1186/1471-2164-10-527

PubMed Abstract | CrossRef Full Text | Google Scholar

Sturny, R., Cam, K., Gutierrez, C., and Conter, A. (2003). NhaR and RcsB Independently Regulate the osmCp1 Promoter of Escherichia coli at Overlapping Regulatory Sites. J. Bacteriol. 185 (15), 4298–4304. doi:10.1128/JB.185.15.4298-4304.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, L., and Wu, R. (2015). Mapping Complex Traits as a Dynamic System. Phys. Life Rev. 13, 155–185. doi:10.1016/j.plrev.2015.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Svenning, J.-C., Gravel, D., Holt, R. D., Schurr, F. M., Thuiller, W., Münkemüller, T., et al. (2014). The Influence of Interspecific Interactions on Species Range Expansion Rates. Ecography 37 (12), 1198–1209. doi:10.1111/j.1600-0587.2013.00574.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tewa, J. J., Djeumen, V. Y., and Bowong, S. (2013). Predator-Prey Model with Holling Response Function of Type II and SIS Infectious Disease. Appl. Math. Model. 37 (7), 4825–4841. doi:10.1016/j.apm.2012.10.003

CrossRef Full Text | Google Scholar

Vellend, M. (2010). Conceptual Synthesis in Community Ecology. Q. Rev. Biol. 85 (2), 183–206. doi:10.1086/652373

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., Pang, X., Lv, Y., Xu, F., Zhou, T., Li, X., et al. (2013). A Dynamic Framework for Quantifying the Genetic Architecture of Phenotypic Plasticity. Brief. Bioinform. 14 (1), 82–95. doi:10.1093/bib/bbs009

PubMed Abstract | CrossRef Full Text | Google Scholar

West, G. B., Brown, J. H., and Enquist, B. J. (2001). A General Model for Ontogenetic Growth. Nature 413 (6856), 628–631. doi:10.1038/35098076

PubMed Abstract | CrossRef Full Text | Google Scholar

Whitham, T. G., Bailey, J. K., Schweitzer, J. A., Shuster, S. M., Bangert, R. K., Leroy, C. J., et al. (2006). A Framework for Community and Ecosystem Genetics: from Genes to Ecosystems. Nat. Rev. Genet. 7 (7), 510–523. doi:10.1038/nrg1877

PubMed Abstract | CrossRef Full Text | Google Scholar

Wimp, G. M., Tomasula, J., and Hamilton, M. B. (2019). Putting the Genes into Community Genetics. Mol. Ecol. 28 (19), 4351–4353. doi:10.1111/mec.15209

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, R., Cao, J., Huang, Z., Wang, Z., Gai, J., and Vallejos, E. (2011). Systems Mapping: How to Improve the Genetic Mapping of Complex Traits through Design Principles of Biological Systems. BMC Syst. Biol. 5 (1), 84–11. doi:10.1186/1752-0509-5-84

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, R., and Lin, M. (2006). Functional Mapping - How to Map and Study the Genetic Architecture of Dynamic Complex Traits. Nat. Rev. Genet. 7 (3), 229–237. doi:10.1038/nrg1804

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, R., Ma, C.-X., Hou, W., Corva, P., and Medrano, J. F. (2005). Functional Mapping of Quantitative Trait Loci that Interact with thehgMutation to Regulate Growth Trajectories in Mice. Genetics 171 (1), 239–249. doi:10.1534/genetics.104.040162

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, G. Z., Chen, L. S., and Sun, L. H. (2005). Complexity of an SIR Epidemic Dynamics Model with Impulsive Vaccination Control. Chaos, Solitons & Fractals 26 (2), 495–505. doi:10.1016/j.chaos.2005.01.021

CrossRef Full Text | Google Scholar

Zeng, J., de Vlaming, R., Wu, Y., Robinson, M. R., Lloyd-Jones, L. R., Yengo, L., et al. (2018). Signatures of Negative Selection in the Genetic Architecture of Human Complex Traits. Nat. Genet. 50, 746–753. doi:10.1038/s41588-018-0101-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, S., Dong, L., and Chen, L. (2005). The Study of Predator-Prey System with Defensive Ability of Prey and Impulsive Perturbations on the Predator. Chaos, Solitons & Fractals 23 (2), 631–643. doi:10.1016/j.chaos.2004.05.044

CrossRef Full Text | Google Scholar

Zhao, W., Chen, Y. Q., Casella, G., Cheverud, J. M., and Wu, R. (2005a). A Non-stationary Model for Functional Mapping of Complex Traits. Bioinformatics 21 (10), 2469–2477. doi:10.1093/bioinformatics/bti382

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, W., Hou, W., Littell, R. C., and Wu, R. (2005b). Structured Antedependence Models for Functional Mapping of Multiple Longitudinal Traits. Stat. Meth. Mol. Genet. Biol. 4 (1). doi:10.2202/1544-6115.1136

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, X., Tong, C., Pang, X., Wang, Z., Guo, Y., Du, F., et al. (2012). Functional Mapping of Ontogeny in Flowering Plants. Brief. Bioinform. 13 (3), 317–328. doi:10.1093/bib/bbr054

PubMed Abstract | CrossRef Full Text | Google Scholar

Zimmerman, D. L., and Núñez-Antón, V. (1997).Structured Antedependence Models for Longitudinal Data Modelling Longitudinal and Spatially Correlated Data. New York, NY: Springer, 63–76. doi:10.1007/978-1-4612-0699-6_6

CrossRef Full Text | Google Scholar

Zwietering, M. H., Jongenburger, I., Rombouts, F. M., and van 't Riet, K. (1990). Modeling of the Bacterial Growth Curve. Appl. Environ. Microbiol. 56, 1875–1881. doi:10.1128/aem.56.6.1875-1881.1990

PubMed Abstract | CrossRef Full Text | Google Scholar