In order to study the competition and cooperation mechanism of two species in the same environment, we, respectively, collected 100 strains of E. coli and S. aureus and randomly paired these strains to form a total of 100 independent interspecific combinations and numbered them. The source information of 200 strains can be found in Supplementary Table S1 and the relevant sequencing information has been uploaded to the National Genomics Data Center (NGDC) database. Two hundred original strains were cultured on the tryptose soya agar (TSA) solid medium (OXOID, Basingstoke, England) for activation treatment. A total of 100 conical flasks (50 mL) were, respectively, added to 25 mL of the same two-times dilution of brain-heart infusion (BHI) medium (OXOID, Basingstoke, England). We picked out the single colonies of two bacteria on the TSA culture plate and then inoculated them to 100 BHI culture media for co-culture according to the numbered pairing principle. The absorbance of the bacterial solution was measured with an enzyme-labeled instrument and converted to the bacterial solution concentration, and the bacterial solution was diluted to 5 × 10³ cells/mL, and the flask was cultured continuously for 36 h in a constant temperature shaker at 30°C and 130 r/min. This co-culture experiment was replicated at least three times. The method of randomized block experiment was adopted to eliminate the random error which may be produced by the environment. Faced with the pressure of limited resources, two strains from different species in a culture flask may choose to compete or cooperate.

According to the law of bacterial growth, the strain will show a trend logarithmic growth before the stable period, 1 mL of samples were taken every 0.5 h, 2 h, 4 h, and 6 h at 0–2 h, 2–12 h, 12–24 h, and 24–36 h, respectively, and then bacterial genomic DNA was extracted by using the TIANamp Bacterial Genome Extraction Kit. The specific primers were designed according to the uidA (1694260–1,696,071) gene encoding β-D-glucosidase in E. coli and the nuc (1397756–1,398,289) gene encoding thermostable nuclease in S. aureus. The specific fragments obtained by primer amplification were cloned into the pMD18-T vector for preparation of standards, and then these standards were further diluted by 10 times gradient to obtain standard curves. Through the Mx3005P real-time fluorescence quantitative system (Stratagene, La Jolla, United States), two specific primers were used to perform qPCR amplification on E. coli and S. aureus to obtain CT values. These values were then substituted into the standard curve to calculate the bacterial abundance of each bacteria in each flask at each time point, which was repeated three times to take the average value for genetic mapping. The total volume of the qPCR reaction system was 25 μL, including 2 × SuperReal PreMix Plus (SYBR Green I) (TIANGEN, Beijing, China), 300 nM forward primer, and 300 nM reverse primer. The thermal cycle reaction conditions were as follows: initial denaturation at 95°C for 10 min, followed by 40 cycles of 30 s at 95°C, 1 min at 55°C, and 1 min at 72°C. Fluorescent signals were detected and collected during annealing and extension.

Using E. coli str. K-12 substr. MG1655 and S. aureus subsp. aureus NCTC 8325 as reference strains, the whole genome sequencing of 100 strain combinations was performed on the Illumina HiSeq 2,500 (Novogene, Beijing, China) platform to obtain SNPs (single-nucleotide polymorphisms) genotype data of two species at the genome-wide level. The obtained effective sequencing data was performed for the sequence alignment with the reference genome using BWA (Li et al., 2009) software (version 0.7.17) and for quality control to ensure the SNPs data with high-quality scores (Q value ≥20) and enough supporting bases (≥ 4, with variation) using SAMtools (Li et al., 2009; Danecek et al., 2021) software (version 1.17). Ultimately, a total of 745,528,965 SNP combinations were obtained from the E. coli and S. aureus genomes, which should be sufficient enough in density to identify genomic regions in the QTL mapping.

To obtain high-quality SNP combinations and improve systems mapping efficiency, this study adopted a two-step approach to filter out low-quality SNP combinations: (1) Before combining the SNPs of two bacteria, the single marker analysis method was used to control the number of SNPs, and T-test was used to filter out SNPs below the threshold; (2) After pairing the screened SNPs, SNP combinations were further filtered by comparing the proportions of the four genotypes in each combination. SNP combinations with genotype proportions below 10% were filtered out. Based on the above method, 272,873 high-quality SNP combinations were ultimately obtained, greatly improving the efficiency of parameter estimation for systems mapping.

Population genetic structure analysis can provide information on the origin and composition of individuals and is a worthwhile tool for analyzing genetic relationships. Based on the posterior variational Bayes framework, fastStructure software (Raj et al., 2014) is able to calculate the estimated values faster than the classical population structure analysis software STRUCTURE, which is suitable for inferring the population structure from large SNP genotype data. In this study, we conducted population structure analysis on SNPs data of two species using fastStructure software (version 1.0). Based on the default convergence criteria and prior probability, K values ranging from 2 to 20 were tested with 10 replicates per K. Using the chooseK.py function in fastStructure, a reasonable range of K values was determined. Meanwhile, the prcomp function in the R software was performed for principal component analysis (PCA) to further investigate the structural components of the two bacteria and their relationship. The emma.kinship function in the R package emma was used to calculate the kinship matrix of two bacterial SNPs used in subsequent analysis (Kang et al., 2008). According to the analysis of kinship and population structure, we used the lmer function in the R package lme4 (Bates et al., 2015) to adjust the phenotype of two bacteria based on the linear mixed effect model, to avoid detecting the phenotypic data of the falsely associated population structure.

Assume that there are two QTL that affects the growth trait of two bacteria in co-culture, including two genotypes A and a by species E and two genotypes B and b by species S, respectively. These two QTLs formed four interspecific genotype combinations, expressed as A / B , A / b , a / B , and a / b . By correlating SNPs genotype data and dynamic phenotypic data of n interspecific strains combinations, we used maximum likelihood estimation (MLE) to construct a dynamic model for parameter estimation of specific genotypes. A mixture-based likelihood model has been widely used to map QTLs for complex traits. The mixture-based likelihood of time-dependent abundance data for n interspecific strains combinations is formulated as

where N i E = ( N i E ( 1 ) , … , N i E ( T ) ) and N i S = ( N i S ( 1 ) , … , N i S ( T ) ) are the vectors of abundance trajectories at T times for species E and S, respectively, and f · ( N i E ,N i S ; Θ , Ψ ) is a multivariate normal distribution with expected mean vector (6 parameters) and covariance matrix Σ (5 parameters) for pair i that belongs to a particular genotype–genotype combination. The covariance matrix can be simulated using the second-order structured antedependence (SAD(2)) model (Zimmerman et al., 2001; Zhao et al., 2005).

Based on the dynamic likelihood model (1), it becomes a crucial issue how to map the significant QTLs related to the interaction between two species. We can test whether there are significant interspecific QTLs involved in interspecific interactions, which can be done by formulating the following two assumptions:

where H 0 is the null hypothesis that each component uses the same ODE parameters in different genotypes. H 1 is an alternative hypothesis, which means that there is at least one component in the system that has different ODE parameters under different genotypes. Using the values of H 0 and H 1 to calculate the likelihood values L 0 and L 1 , respectively, and further calculate their log-likelihood ratio (LR), which can be expressed as

These LR values are compared to genome-wide critical thresholds determined by 1,000 permutation tests. If H 0 is rejected, it means that significant interspecific QTLs from both species have been mapped.

After mapping all the significant QTLs among species, we further implemented the corresponding functional annotations, so as to obtain the key genes that play a critical role in the interaction process of different species in the system, which can provide a theoretical basis for the study of competition and cooperation mechanism among species. Meanwhile, we proofread the results of GO annotation via MGI¹ and AureoWiki² databases.

We used a set of generalized Lotka-Volterra (LV) ordinary differential equations to describe the interaction between E. coli and S. aureus in co-culture (Fujikawa et al., 2014; Jiang et al., 2018a,b). The LV equations were integrated into the systems mapping model, which can screen out significant QTLs for the interaction between the two species in co-culture. We combined this new approach with community ecology so that patterns of interactions between two species can be quantified and explained. The LV equations could be divided into two different parts to describe microbial abundance differently, expressed as

where r E and r S are the Malthusian growth rates; K E and K S are the environmental capacities of the two different species; and α E|S and α S|E are dimensionless parameters used to simulate how one species affects the other through the interaction in co-culture; M ˙ E and M ˙ S represent the independent growth of each species, which is determined by its inherent attributes; N ˙ E|S and N ˙ S|E represent the interactive growth of each species, which is determined by how the species interacts with the other in co-culture. If the interactive growth of a species is positive or negative, it indicates that the species is beneficially or detrimentally affected by the other. If the interactive growth is zero, it means that the two species will not affect each other. Therefore, by estimating ODEs parameters Θ = ( r E , K E , α E|S ; r S , K S , α S|E ) , the LV equations can not only specify the dynamic pattern of each species abundance, but also describe two interaction patterns of two species in co-culture.

According to quantitative genetic theory (Li and Wu, 2009), four interspecific genotypic values can be partitioned into different components and then converted into various genetic effect components, expressed as

Where μ A B l ( t ) , μ A b l ( t ) , μ a B l ( t ) and μ a b l ( t ) are the genotypic values for the abundance of each species l ( l = E o r S ) at time t ; μ l ( t ) is the time-dependent population means of abundance of each species l ( l = E o r S ) ; if l = E , a A l ( t ) is the time-dependent direct genetic effect of species E’s QTL on its own abundance, while if l = S , a A l ( t ) is the time-dependent indirect genetic effect of species E’s QTL on the abundance of its coexisting species S; if l = S , a B l ( t ) is the time-dependent direct genetic effect of species S’s QTL on its own abundance, while if l = E , a B l ( t ) is the time-dependent indirect genetic effect of species S’s QTL on the abundance of its coexisting species E; and I A × B l ( t ) is the time-dependent genome–genome epistatic effect between the QTLs of two species on the abundance of each species l ( l = E o r S ) .

After mapping significant interspecies QTLs, the next step is to test whether these QTLs exert significant direct genetic effects, indirect genetic effects, and interspecies genome-genome epistasis effects. The null hypotheses for these tests are expressed as

For a given QTL, the systems mapping models can estimate its time-dependent direct genetic effects on independent growth, time-dependent indirect genetic effects on interactive growth, and time-dependent genome-genome epistasis effects on interactive growth. The procedure of constructing the corresponding genetic network of all QTLs based on their genetic effects is as follows: Let q m ( t ) denotes the genetic effects of QTL m ( m = 1 , … , M ) at time t . All these M QTLs interact with each other in a network to jointly affect growth. Through a set of ordinary differential equations, q m ( t ) is used to denote the genetic effects of other QTLs, expressed as

where w m ( q m ( t ) ) is the parametric or non-parametric function of q m ( t ) , and u m ′ ( q m ′ ( t ) ) is the parametric or nonparametric function of the genetic effects of QTL m ′ ( m ′ ≠ m , m ′ = 1 , … , M ) . The latter presents the interactive effects of other QTLs on the focal QTL. Non-parametric methods (Wu et al., 2014) are used to estimate the parameters of these functions, so that the corresponding genetic effect network can be constructed among these QTLs, and how a QTL interacts with others through activation or inhibition, thereby affecting the expression of the latter.

According to the marginal possibility, we identified eight and nine subpopulations in E. coli and S. aureus strains, respectively (Supplementary Figure S1A). In terms of PCA analysis, the first two PCs of E. coli accounted for 50.4 and 11.5% of the variation in the genotypic data respectively, and PC1 and PC3 of S. aureus accounted for 47.5 and 10.0%, respectively (Supplementary Figure S1B). The results of fastStructure and PCA indicated significant population structure differences among individuals of the natural population of these two bacteria. The confounding effects of population structure and kinship were well removed to a certain extent (Supplementary Figures S1C, S2, S3). The adjusted phenotypic data can be used for systems mapping of the two species in co-culture.

Taking the number of microbes (the number or weight of the living bacteria) as the vertical coordinate, and the cultivation time as the horizontal coordinate, we can draw a microbial growth curve. A typical microbial growth process includes three successive growth phases: lag, logarithmic, and stationary phases. At the lag phase, the population adapts to a new environment, in which the specific growth rate gradually increases from zero and then accelerates the maximum value at the logarithmic phase. At the stationary phase, however, the rate decreases to zero and the population growth terminates, thereby reaching the culture-carrying capacity. A number of growth equations have been derived to capture these phases (Zwietering et al., 1990; Tonner et al., 2020). In this study, a set of generalized LV ordinary differential equations was used to fit the growth curve of E. coli and S. aureus in co-culture during the first 36 h after culture (Figure 1A). For the growth curve of each strain, the actual individual growth curve and the fitted growth curve of 100 pairs of interspecific strains of E. coli and S. aureus were shown in Supplementary Figures S4, S5. In the average growth curve, we can observe that the fitting of the two bacteria met the analysis requirement, with a coefficient of variation R² ≥ 0.94, which indicated that the growth curve fitted by LV can well represent the dynamic phenotype data of the interaction between the two bacteria. Overall, the microbial abundance of E. coli in the unit time is generally higher than S. aureus. For example, at 0.5 h, the abundance value of E. coli is over e 16 , while that of S. aureus is less than e 14 . Selecting 0.5 h, 6 h, 16 h, and 36 h to draw the phenotype distribution of the two bacteria (Figure 1B), we could see that the vast majority of the E. coli strains grow faster than S. aureus strains, but there were also some outliers indicating the existence of some opposite situations. We selected four kinds of representative growth curves (Figure 1C) and found that, in addition to the fact that E. coli has consistently been leading the growth of S. aureus or vice versa, there were also two bacteria with similar initial growth conditions, but during the logarithmic period, the growth of E. coli exceeded that of S. aureus and even S. aureus was severely inhibited. This indicates that there may be some genes that affect the actual growth of the two bacteria when they interact, and the degree of influence varies greatly between different bacterial pairs.

The correlation analysis between the phenotypic data of two bacteria in co-culture and SNP genotype data by using the systems mapping theory can excavate the QTL closely related to the interspecific interaction. The LR threshold at the genome-wide level was determined by 1,000 permutation tests, and we could draw the Manhattan diagram (Figure 2) and search out significant QTL combinations affecting the microbial abundance through the interspecific interaction, whose LR values are greater than the LR threshold. The LR threshold of this study was 450.0346, and 244 vital QTL combinations were obtained through comparison, including 69 SNPs from E. coli and 59 SNPs from S. aureus (Supplementary Tables S2, S3). These SNPs were further mapped to genes, and 57 and 43 genes were obtained in E. coli and S. aureus, respectively (Supplementary Tables S4, S5).

Based on the results of 3.1, we successfully annotated 57 genes in E. coli with 44 effective GO annotation results and 43 genes in S. aureus with 38 effective GO annotation results. Supplementary Tables S4, S5, respectively, listed the detailed functional annotation information of the significant interspecific interaction QTLs of the two bacteria.

GO annotations can be divided into three categories: Biology Process (BP), Cellular Component (CC), and molecular function (MF). Each large category can be subdivided into subcategories, which can help us understand the biological significance behind each gene (Figures 3, 4).

In the gene annotation of E. coli (Figure 3), the BP category involved 37 genes, accounting for the highest proportion (65%), followed by the CC category (31 genes, 54%). The BP category mainly included DNA-templated transcription (5 genes), anaerobic respiration (3 genes), phosphorus metabolic process (3 genes), and sulfur metabolic process (3 genes); the CC category mainly included peptidoglycan-based cell wall (13 genes), cytoplasm (12 genes), and organelle inner membrane (8 genes); the MF category includes transcription activator activity (3 genes) and transcription repressor activity (2 genes).

In the gene annotation of S. aureus (Figure 4), the CC category involved 29 genes, accounting for the highest proportion (67%), followed by the MF category (23 genes, 53%). The CC category was mainly related to membranes, such as membrane (14 genes) and plasma membrane (12 genes); the MF category included oxidoreductase activity (4 genes), ATP binding (3 genes), and catalytic activity (3 genes); the BP category mainly included pathogenesis (2 genes) and transmembrane transport (2 genes).

Any species in co-culture may choose to cooperate or compete with other species from the same system, depending on the level at which common resources can be shared for their respective growth. According to the community ecology theory, the interaction types of two different species can be formulated by a strategy matrix (Supplementary Figure S6), including mutualism, competition, neutralism, commensalism, predation/parasitism, and amensalism.

Different strategies used between two species in one system resulted in six distinct types of interactions: (++) → Mutualism, which means that two species benefit from each other; (00) → Neutralism, which means that any one species does not depend on or affect the other; (−-) → Competition, which means that two species fight each other; (+0) → Commensalism, which means that one species is beneficial to the other, while the latter does not affect the former; (+−) → Predation/Parasitism, which means that one species helps the other, whereas the latter is harmful to the former; (−0) → Amensalism, which means that one species harms the other, while the latter does not affect the former.

Based on the LV equations, the overall growth of each species in co-culture can be decomposed into independent growth and interactive growth. By estimating the ODEs parameters Θ = ( r e , K e , α e|s ; r s , K s , α s|e ) , we can characterize the interaction mode of two species in co-culture by comparing the interactive growth between species. We plotted the overall growth curves (solid lines), independent growth (long dashed lines), and interactive growth curves (short dashed lines) of E. coli and S. aureus over time in co-culture, according to ODE parameters (Figure 5A). On the whole, the interactive growth was always negative because two bacteria in the same environment compete with each other for limited resources, and the independent growth is much greater than the actual overall growth observed, indicating that these two bacteria have a competitive relationship in co-culture.

Each interspecies QTL combination inherited two bacterial alleles 0 and 1, respectively, forming four interspecies genotype–genotype combinations, denoted as 0/0, 0/1, 1/0, and 1/1. Taking the selected QTL combination E3414424 × S383182 (LR value = 505.2964) as an example to analyze the interaction between two bacteria in co-culture, we drew the microbial growth curves of four genotype combinations of two bacteria in co-culture (Figure 5B). By combining the relevant parameters of ODEs, the overall growth curve (solid line) of microbial abundance for each genotype combination can be decomposed into independent growth (long dashed line) and interactive growth curve (short dashed line). For each genotype combination, it can be observed that the competition between E. coli and S. aureus was intense, and the overall growth was lower than the independent growth, but the overall growth, independent growth, and interactive growth were significantly different in scale and pattern among all four combinations. These results can be summarized as follows: (1) The two bacteria in the same system were in an antagonistic relationship, but the degree of antagonism of S. aureus to E. coli was greater than that of E. coli to S. aureus; (2) When the same genotype was paired, the antagonistic was expressed more strongly, for example, 0/0 and 1/1 exhibited larger negative interactive growth than 0/1 and 1/0; (3) Combination-dependent differences in independent growth and interactive growth was more pronounced than in overall growth, suggesting that bacteria exhibited more variation in their internal machinery than what can be phenotypically observed. These results suggested that E3414424 and S383182 were antagonistic QTLs that participated in determining and shaping the antagonistic relationship between E. coli and S. aureus when they were in the same system.

Interspecies QTL combinations can affect the independent growth of two bacteria as a growth potential and also affect their interactive growth determined by interacting with other species in the ecosystem. To some extent, such QTLs play a radical role in regulating microbial communities’ dynamic structure and behavior. The systems mapping model can also dissect the genetic structure of community traits into its direct, indirect, and genome-genome epistatic genetic effects. The direct genetic effect of an interspecific QTL combination describes how it affects the independent growth of one species (such as E. coli) where this QTL is mapped. The indirect genetic effect specifies the genetic effect of an interspecific QTL combination from one species (such as S. aureus) on the interactive growth of the other (such as E. coli) competing with this species in the same system. And the genome-genome epistatic effect depicts the genetic effect of the interaction between the alleles of two species on the interactive growth of one species.

After estimating and mapping these three genetic effects of the interspecies QTL combination (E3414424 × S383182) on the growth trajectory of two bacteria, it was found that the indirect effects and genome-genome epistasis effects, which were ignored in previous studies, were significant, and even their influence range was greater than the direct effect (Figure 6). E3414424 from E. coli appeared to be a more “aggressive” QTL, because its indirect effect on the abundance dynamics of S. aureus was greater than that of S383182 from S. aureus on the abundance dynamics of E. coli. Meanwhile, its direct effect on the abundance of its home species at the logarithmic phase was greater than the indirect effect of S383182 on its abundance. From the analysis of the genetic effect curve, we can see how an SNP affects the growth of two microbes over time. For this particular QTL combination, the direct, indirect, and genome-genome epistatic effects on E. coli reached their maximum at 10–14 h, and these effects on S. aureus showed a similar pattern.

The bacterial abundance of E. coli in co-culture was jointly determined by the direct effects of its own 60 QTLs, the indirect effects of 49 QTLs from S. aureus, and the interspecific epistatic effects between 230 QTL combinations. The bacterial abundance of S. aureus in the same system was determined by the direct and indirect effects triggered by its own 49 QTLs and 60 E. coli QTLs, respectively, and the interspecific epistatic effects among 230 QTL combinations. We employed an ODE-based approach to map genetic networks to characterize how these significant QTLs interacted in the network to affect the microbial abundance through three different types of effects (Figure 7). By comparing the three types of QTL networks, we found that the structure and organization of these networks differed greatly between the two bacteria. The core feature of the QTL network is its ability to identify the hub QTLs, which play a vital role in the genetic architecture of microbial growth in the same system.

For E. coli, there were a major hub QTL SNP 50 (E4323969) and two minor hub QTLs SNP 28 (E2357245) and 54 (E4544464) in the direct effect network that directly affected its own growth (Figure 7A), while S. aureus’s direct effect network had SNP 21 (S301259), 38 (S590022), 11 (S134253), and 27 (S321552) as the hub QTLs directly affecting its own growth. In the indirect effect network, there were 3 hub QTLs from S. aureus which exerted an indirect effect on the growth performance of E. coli, including SNP 47 (S1731046), 22 (S301667), and 24 (S306067). However, there were 5 hub QTLs from E. coli in the network that indirectly affected the growth performance of S. aureus, including SNP 16 (E2015702), 41 (E4125943), 8 (E555930), 30 (E2683642), and 28 (E2357245; Figure 7B). Interestingly, whether directly affecting its own growth or being indirectly affected by competing strains, the number of hub QTLs in S. aureus was greater than that in E. coli. SNP 28 (E2357245) from E. coli acted as a hub QTL in both direct-to-self and indirect-to-S. aureus effect networks, although the network structures of these two effects were very different. In terms of the genome-genome epistasis effect network, E. coli had 10 QTL combinations acting as hubs, and S. aureus had 8. Surprisingly, the 34th QTL combination (E2015702 × S1385806) acted as a major hub in E. coli, whereas it was a minor hub for S. aureus (Figure 7C).