Genetic Control of Rhizosphere Microbiome of Cotton Plant Under Field Conditions

doi:10.21203/rs.3.rs-840382/v1

Background: Understanding the extent of heritability of plant-associated microbiome (phytobiome) is critical for exploitation of phytobiome in agriculture. Two crosses were made between pairs of cotton cultivars with differential resistance to the Verticillium wilt. F₂ plants were grown together with the four parents to study the heritability of cotton rhizosphere microbiome. Amplicon sequencing was used to profile bacterial and fungal communities.

Results: F₂ hybrid offspring of both crosses had higher alpha diversity indices as well as greater variability than the two parents. Two types of data were used generated to study the heritability of rhizosphere microbiome: principal components (PCs) and individual top microbial OTUs. For the L1×Z49 cross, the variance among the 100 progeny genotypes (V_T) was significantly greater than the random variability (V_E) for 12 and 34 out of top 100 fungal and bacterial PCs, respectively. For the Z2×J11 cross, the corresponding values were 10 and 20 PCs. Out of the top 100 OTUs, 29 fungal OTUs and 10 bacterial OTUs had V_T that was significantly greater than V_E for the L1×Z49 cross; the corresponding value for the Z2×J11 cross were 24 and one. The estimated heritability was mostly in the range of 40% to 60%.

Conclusions: These results suggested the existence of genetic control in terms of rhizosphere microbiome in cotton.

General Microbiology

Cotton

microbiome

Rhizosphere soil

heritability

Verticillium wilt

Plant associated microbiome (phytobiome) may benefit their host plants in several aspects [1], including protection against pathogen and pests [2–5], increased tolerance to drought and nutrient stress [6], changes in flowering time [7], and enhanced plant productivity [8]. Microbiome structure in the rhizosphere often differ across plant species [9], as well as among genotypes within a single species [5]. Recent studies have shown that rhizosphere microbiome are shaped, to a certain extent, by host genetics [10–12]. A study comparing the root microbiomes of cereal crops shows that there is a strong correlation between host genetic differences and microbiome composition, which indicates that a subset of phytobiome may be affected by host genotypes across a series of plant hosts [9]. Microbiome can function as a phenotypically plastic buffer between host genotype and environment, and interact with environment to shape host phenotypes [13]. Therefore, expression of almost any host phenotype depends to some extent on phytobiome composition [13]. Although there is consistent evidence of interactions between host genotype and phytobiome composition, identifying specific genetic factors driving/controlling phytobiome acquisition and assembly remains a challenge.

The main focus of plant breeding has been on the improvement of crop productivity, quality and pest/disease resistance; modern breeding has been shown to reduce genetic diversity in modern crops [14,15]. This may have unintended consequences on phytobiome associated with the reduced crop diversities. The importance of the rhizosphere microbiome in the plant ecosystem functioning has been widely recognized, but plants have been bred by altering their genomic information with little consideration of their interaction with surrounding organisms. Recently, there has been a paradigm shift in considering plants as a holobiont, an ecological and evolutionary unit containing both the host and its microbiome [16]. Thus, we may need to include the ability of recruiting and interacting with beneficial microbes through root exudates as a selection criterion [17,18]. Mendes et al. [19] have shown that resistance breeding in common bean has unintentionally co-selected for plant traits that strengthen the rhizosphere microbiome network structure and enrich for specific beneficial bacterial genera that express antifungal traits involved in plant protection against infections by soilborne pathogens. Identifying these plant traits and microbial taxonomies will help breeders to select for plant traits that enrich desired microbial groups [19]. However, to pursue this route, we also need to understand the heritability of microbiome.

Heritable components of microbiome variation may reflect plant-microbe interactions, which are formed through natural selection acting on plant traits underlying fitness. Recently, Wagner et al. [20]found that inbred lines and hybrids differ consistently in composition of bacterial and fungal rhizosphere communities of maize. Most studies on microbiome composition in complex environments have observed that the interaction between genotype and environment are at least as strong as the main genotypic effects [11,20–22]. Currently, most studies on host genetic control of microbiome have shown that there are differences in microbiome composition among varieties or accessions with no well-defined genetic relationship to each other, often chosen to represent a breadth of diversity within the host species [23]. Thus, there is a lack of purposely designed studies to investigate heritability of phytobiome through designed crosses between plant genotypes.

Cotton Verticillium wilt, caused by Verticillium dahliae, is one of the devastating plant diseases worldwide. Because of the inaccessibility of V. dahliae during infection, long-term survival of its microsclerotia in soil, its broad host range and withdrawal of broad-spectrum soil fumigants, it is difficult to control Verticillium wilt [24]. Recently, we demonstrated that genotypic response to V. dahliae in cotton is associated with many microbial groups in rhizosphere and endophytic microbiomes [5], suggesting the potential of improved phytobiome through breeding for wilt resistance. In the present study, we selected two pairs of cotton cultivars with differential resistance to V. dahliae as parents [L1 (susceptible) × Z49 (resistant) and Z2 (resistant) × J11 (susceptible)] to investigate to what extent is the rhizosphere microbiome is heritable and hence open to breeding. We used amplicon sequencing to profile bacterial and fungal communities in 100 F₂ offspring from each cross hybrids and the four parental genotypes, grown together in the field.

Plant materials

Two resistant and two susceptible cotton cultivars were chosen as parents for crossing: Gossypium hirsutum cultivars, Zhongzhimain2 (Z2, resistant against V. dahliae) crossed with Jimian11 (J11, susceptible), and LocalVIR875-1 (L1, susceptible) with Zhongmiansuo49 (Z49, resistant) (Fig. 1A). F₂ segregating populations were obtained from the two crosses and used to study heritability of rhizosphere microbiome.

Site description and field experiment design

In 30 April 2019, a field experiment was set up at the Institute of Cotton Research, Chinese Academy of Agricultural Sciences (Anyang, China) (Fig. 1C). The field, was used for evaluating cotton cultivar resistance against V. dahliae (disease nursery), were artificially inoculated with the pathogen 20 years ago. The soil at the experimental site (36˚03′44″N, 114˚28′52″E) is classified as cambisol type soil [25]. A completely randomized block design with four blocks was used. Twenty-five F₂ plants from each cross were randomly allocated to block. Within each block, there were six plots, each was 5 m long with two rows (0.8 m between two rows); neighbouring plots were separated by 0.8 m; each of the six plots was randomly assigned to one of the six groups of plants: four parents and two F₂ populations. There were 100 F₂ plants for each cross. On 22 August 2019, wilt severity of all individual plants was recorded on a scale of 0 to 4 as described previously [5]. For cv. Z2, there were 68, 19, 11, 2 and 0 plants with severity score of 0, 1, 2, 3 and 4, respectively; the corresponding values for cv. J11 were 8, 8, 22, 51 and 11 plants, for cv. L1 were 15, 13, 25, 39 and 8 plants, and for cv. Z49 were 52, 29, 15, 4 and 0 plants (Table S1). For Z2×J11 cross, there were 23, 27, 28, 17 and 5 plants with severity score of 0, 1, 2, 3 and 4, respectively; the corresponding values for L1×Z49 cross were 31, 34, 19, 14 and 2 plants (Fig. 1B; Table S1).

Rhizosphere sample collection

During late August (at the boll-forming stage) (Fig. 1C), approximately 16 weeks after sowing, rhizosphere soil samples were collected and stored as described previously [5]. For F₂, all plants were sampled; for each of the four parents, only four plants were sampled from each block, giving 16 plants per parental genotype. Each plant was carefully removed from the soil with a spade. Root system of the plant were first vigorously shaken to remove loosely adhering soil particles. Plant fine root were cut into pieces of approximately 2 cm length with a pair of sterile scissors. Rhizosphere soil samples were harvested in 500-ml screw-cap bottles with ca. 20 g roots. Each bottle was filled up to 300 ml with 1:50 TE buffer (1 M Tris, 500 mM EDTA, and 1.2% Triton diluted in sterile distilled water) and shaken at 270 rpm for 1 h. The root-washing suspension was filtered with sterile cheesecloth and centrifuged at 4,000 ×g for 20 min [5]. The supernatant was discarded by pipetting. This procedure was repeated three times before the pellets were re-suspended in the remaining solution, transferred to a 2-ml Eppendorf tube and centrifuged at 14,000 ×g for 20 min. The pellets were immediately frozen and stored at -80˚C before DNA extraction.

DNA extraction and sequencing

Soil pellets were re-suspended in 500 μl MoBio PowerSoil bead solution, and DNA was extracted from the resulting pellets (250 mg) using the MoBio PowerSoil DNA Isolation Kit (MoBio Laboratories, Carlsbad, CA, USA) according to the manufacturer’s protocol. The extracts were checked on a 1% agarose gel and DNA concentration was estimated by NanoDrop ND-2000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). DNA was stored at -80°C until further analysis.

For bacteria, the V3-V4 hypervariable region of the 16S rRNA gene was amplified in triplicates for each sample using the 341F/805R primers [26]. For fungi, primers ITS5/ITS2 [27] were used to amplify the ITS1 region in triplicates for each sample. PCR reactions and the extraction and purification of amplicons followed previously established methods [5]. Sequencing libraries were generated with the Ion Plus Fragment Library Kit 48 rxns (Thermo Scientific, USA) following the manufacturer's recommendations. The quality of each library was assessed on a Qubit 2.0 Fluorometer (Life Technologies, USA). Finally, total DNAs were submitted to Novogene Co. Ltd. (Beijing, China) for Illumina sequencing, the libraries were sequenced on an Ion S5™ XL platform (Thermo Fisher Scientific, Waltham, MA) to generate single-end reads. In total 528 libraries were sequenced: 264 rhizosphere samples (2 crosses × 100 F₂ plants + 4 parents × 16 replicates) each for 16S rRNA gene and ITS rRNA gene.

Sequencing processing and taxonomy assignment
The clear reads from the service provider were used to generate OTUs and assign taxonomy, following an established pipeline [5]. Briefly, high-quality sequences were obtained by quality control and filtering of sequence quality with very stringent criteria followed our previous publication [28] and was carried out separately for the two type of data sets (16S and ITS). High quality sequences were first dereplicated and unique sequences with only 1 read were discarded. Then all unique sequence reads were sorted by their respective frequencies and then grouped into operational taxonomic units (OTUs) based on 97% or greater identity. All OTU were processed using the UPARSE pipeline (Version 10.0) [29] unless specified otherwise. The clustering algorithm also removed chimeras. The SINTAX algorithm (https://www.drive5.com/usearch/manual/sintax_algo.html) then assigned each OTU representative sequence to taxonomic ranks by alignment with the gene sequences against Unite V7 fungal database [30] and RDP training set (v16) bacterial database [31]. Then an OTU table (a sample-by-observation contingency table) was generated by aligning all sequences filtered with far less stringent criteria with the OTU representative sequences as described by Deakin et al. [28].

Statistical data analysis

The median-of-ratios method implemented in DESeq2 [31,32] was used to normalise the OTU counts before any statistical analysis. All statistical analyses were carried out in R 4.0.3. There were six groups of samples, classified into two hierarchical levels: two crosses (Z2×J11 and L1×Z49), and three groups within each cross (two parents, each with 16 replicates, and the 100 F₂ progeny plants).

Alpha (a) diversity (Shannon and Simpson) indices were calculated using the R vegan 2.3-1 package [33]. The rank of a diversity indices were subjected to ANOVA to assess the differences between two crosses and between three groups within each cross via a permutation test for significance. Beta (b) diversity indices were not calculated since the main objective of this study was to compare the variance among the 100 progeny samples with the estimated residual variance for estimation of heritability.

To determine whether there is significant genetic variability among the 100 F₂ progeny from each cross, a F-test was carried out to compare two variance estimates: variance [V_T] among the 100 progeny samples, and the random variance [V_E] measuring random environmental errors experienced in the experiment. This random variance was estimated as the residual variance of the ANOVA of the four parents, each with 16 replicates. If the F test showed that V_T is significantly greater than V_E, then the genetic variance [V_G] component was estimated as (V_T – V_E). Consequently, broad sense heritability was estimated as V_G/V_T.

Two datasets were used for analysis of genetic components: principal component (PC) scores and the normalised data of those OTUs with highest counts. To generate PC scores, only those OTUs with highest counts accounting for 99.99% of the total normalised counts were retained for PC analysis (PCA). Before PCA, the normalised counts data were first logarithm transformed on the natural base and then standardised. The 100 OTUs with the highest counts data were selected for genetic component analysis. For each PC or OTU, V_E was estimated from the parent samples and V_T from the 100 progeny samples. Similarly, two single degree contrasts were used to test whether the two parents of each cross differ significantly. Within the analysis of each data type (PC or OTU), the Benjamini-Hochberg (BH) adjustment was used to correct for the false discovery rate associated with the multiple testing. Statistical significance was determined at the 5% level (BH adjusted). BLASTn searches against the GenBank non-redundant database was then used to further characterise those OTUs with significant genetic components.

Sequence quality and generation of OTUs

For L1×Z49 samples (including the two parents), the number of fungal raw reads ranged from 52,163 to 92,596 per sample, with an average of 80,795. Of these fungal raw reads, the number of good quality reads ranged from 49,956 to 89,182 with an average of 77,718. For Z2×J11 samples (including the two parents), the number of fungal raw reads ranged from 52,009 to 89,799 per sample, with an average of 78,051. Of these fungal raw reads, the number of good quality reads ranged from 49,603 to 86,987 with an average of 75,119. Average read length is 221 and 222 bp for the L1×Z49 and Z2×J11 crosses, respectively. There were 1,487 fungal OTUs. Number of sequences classified into fungal OTUs ranged from 31,853 to 77,683 per sample, with an average of 59,781. Sequencing depth is sufficient for all samples (Additional file 1: Fig. S1A). Although there were 1,487 fungal OTUs, most fungal reads came from fewer than 100 OTUs. The most common OTU accounted for 16.7% of the total number of sequences; the top four and 51 fungal OTUs accounted for more than 50% and 90% of the total number of sequences, respectively (Fig. 2A). Nearly 53% of the fungal sequences cannot be assigned to the phylum level at the 90% confidence level (Fig. 2B); most of the other sequences were assigned to Ascomycota (25.5%) or Basidiomycota (20.9%) (Fig. 2B).

For L1×Z49 samples (including parents), the number of bacterial raw reads per sample ranged from 49,813 to 108,445 with an average of 71,309; the number of good quality reads per sample ranged from 39,759 to 88,548 with an average of 60,298. For Z2×J11 samples (including parents), the number of bacterial raw reads ranged from 48,575 to 100,987 per sample (average – 76,404); the number of good quality reads per sample ranged from 38,795 to 87,633 (average – 64,492). The average read length for both crosses is 418 bp. There were 4,196 bacterial OTUs. Number of sequences classified into bacterial OTUs ranged from 27,189 to 77,162 per sample, with an average of 47,371. Sequencing depth is sufficient for all samples (Additional file 1: Fig. S1B). The most common OTU accounted for 23.5% of the total sequences, with the top five and 252 bacterial OTUs accounting for more than 50% and 90% of the total number of sequences, respectively (Fig. 2A). The dominant bacterial phylum is Proteobacteria, accounting for 64.4% of the sequence reads, followed by Actinobacteria (22.0%) and Firmicutes (4.9%) (Fig. 2B).

Alpha diversity of fungal and bacterial communities in rhizosphere

There were significant differences among the two parents and the progeny within crosses for both Shannon and Simpson indices of fungal communities (P < 0.001). However, most of the observed variability in both Shannon and Simpson indices was unexplained – 86.6% (Shannon index) and 89.9% (Simpson index). Overall, the progeny group for both crosses had higher alpha diversity indices of fungal communities as well as greater variability than the two parents (Fig. 3A). The same trend was also present for the observed number of fungal OTUs per sample (Fig. 3A).

As for fungi, there were significant differences among the two parents and the progeny within crosses for both Shannon and Simpson indices of bacterial communities (P < 0.001) although most of the observed variability was not accounted for – 89.8% (Shannon index) and 91.7% (Simpson index). Overall, the progeny group for both crosses had higher alpha diversity indices of bacterial communities than the parents (Fig. 3B). The same trend was also present for the observed number of bacterial OTUs per sample (Fig. 3B). Overall, samples from cultivar J11 had smaller indices well as lower variability than the other groups.

Genetic control of rhizosphere fungal and bacterial microbiome

Of 1,487 fungal OTUs, the top 1,188 OTUs (accounting for 99.99% of the total number of sequences) were kept and subjected to PCA. Although the first PC only accounted for 4.4% of the total variability, the percentage of variability explained by individual PCs initially declined steeply (Fig. 4A). For L1×Z49, the variance of fungal microbiome among the 100 progeny plants (V_T) was significantly greater than the random variability (V_E) for 12 PCs (Table 1). The two parents did not differ significantly for any of the 12 PCs. Fig. 4B gives an example plot for fungal PC1, showing that the progeny displayed a greater variability, particularly towards the high scores. The fungal genetic component as the proportion of the total variability ranged from 0.470 to 0.935. For Z2×J11, V_T was significantly greater than V_E for 10 fungal PCs (Table 1). Only for two of the 10 PCs did the two parents differ significantly; interestingly for these two PCs, the significance level of V_T > V_E was also the greatest among the 10 fungal PCs. As an example, the progeny not only had a higher average fungal PC3 score than the two parents but also displayed a greater variability (Fig. 4C). The estimated heritability ranged from 0.473 to 0.740. The two crosses shared common 6 fungal PCs for which V_T was significantly greater than V_E (Table 1).

Of 4,196 bacterial OTUs, the top 3,493 OTUs (accounting for 99.99% of the total number of sequences) were kept and subjected to PCA. Although the first PC only accounted for 6.6% of the total variability, the percentage of variability explained by individual PCs initially declined steeply (Fig. 4A), more than the rhizosphere fungal data. For L1×Z49, V_T was significantly greater than V_E for 34 bacterial PCs (Table 2) although the two parents did not differ significantly for any of the 34 bacterial PCs. Fig. 4D plots bacterial PC3 scores as an example, showing that the progeny displayed a greater variability, particularly towards the high scores, as well as a higher average. The broad sense heritability ranged from 0.408 to 0.933. For Z2×J11, V_T was significantly greater than V_E for 20 PCs (Table 2), but the two parents did not differ significantly for any of the 20 PCs. Fig. 4E uses bacterial PC5 as an example to illustrate that the progeny not only had a higher average score than the two parents but also displayed greater variability (Fig. 4E). The broad sense heritability ranged from 0.477 to 0.784. The two crosses shared common bacterial 16 PCs for which V_T was significantly greater than V_E (Table 2).

Genetic control of top 100 fungal OTUs and bacterial OTUs

For L1×Z49, V_T was significantly greater than V_E for 29 out of the top 100 fungal OTUs (Table 3), five of which cannot be assigned to the phylum level. The two parents did not differ significantly for any of these OTUs. The broad sense heritability ranged from 0.414 to 0.788. For Z2×J11, V_T was significantly greater than V_E for 24 out of the top 100 fungal OTUs (Table 3), three of which cannot be assigned to the phylum level. The two parents did not differ significantly for any of these OTUs. The broad sense heritability ranged from 0.424 to 0.841. For 14 fungal OTUs, V_T was significantly greater than V_E for both crosses, including Trichoderma brevicompactum and V. dahliae (Table 3).

For L1×Z49, V_T was significantly greater than V_E for 10 out of the top 100 bacterial OTUs (Table 4). broad sense heritability ranged from 0.475 to 0.728. For Z2×J11, V_T was significantly greater than V_E for only one OTU (Achromobacter mucicolens) with its broad sense heritability being 0.615.

The dominant bacterial phylum in the cotton rhizosphere is Proteobacteria, accounting for 64.4% of the sequence reads, consistent with previous findings in the cotton rhizosphere [4, 5]. Proteobacteria are generally adapted to the plant rhizosphere and across different plant species [11] and respond to labile carbon sources, and are usually considered to be r-selected or rapidly growing microbiota of weeds [34]. With regard to rhizosphere fungi, nearly 53% of the fungal sequences cannot be assigned to the phylum level at the 90% confidence level; most of the other sequences were assigned to Ascomycota (25.5%) or Basidiomycota (20.9%), similar to the strawberry rhizosphere [35].

For rhizosphere bacteria and fungi, we demonstrated that F₂ offspring progeny of both crosses had higher alpha diversity indices as well as greater variability than the two parents. A previous study reported that fungal communities in the rhizospheres of hybrid maize plants had higher alpha diversity than inbred lines but no consistent differences in alpha diversity of rhizosphere bacteria between inbred lines and hybrids [20]. In addition, we found significant genetic component for specific rhizosphere microbiome components of cotton, often with moderate to high broad sense heritability (> 0.4), but these components are only a minor proportion of the entire microbiome. This can be seen from the small number of PCs or OTUs with significant genetic components. The maize inbred lines differ significantly in their rhizosphere microbiome but the heritability level was low and the genetic relationship among the inbred lines was not correlated with the diversity characteristics of the rhizosphere microbiome [11]. The overall low heritability of phytobiome may be due to strong environmental effects [23]. Plants are subjected to a variety of biotic and abiotic stresses, which can induce changes in transcriptomics and metabonomics, leading to changes in root and leaf exudates, thereby influencing phytobiomes [36]. Greenhouse or controlled environment growth chambers may be able to regulate climate fluctuations and reduce soil heterogeneity than under field conditions, and thus could lead to greater heritability estimates.

In the present study, transgressive segregation is common for those PCs or OTUs with significant genetic component, indicating a polygenic nature of the genetic control. Peiffer et al. [11] speculated that maize microbial community is controlled by several major genes, rather than many minor genes in the whole genome, which is used to explain the lack of correlation of maize genetic relationships with microbial diversity. Studies with Arabidopsis mutants and transgenes have shown that leaf cuticle characteristics affect microbial colonization of microbes [37], and that physiological traits such as hormone signaling and defensive secondary chemistry affect rhizosphere microbiome [38]. Further research is needed to assess both plant traits and phytobiome characteristics to determine their relationships, which may shed lights on the nature of genetic control of phytobiome by plants as observed in the present study.

The present study showed that the overall plant genetic component was stronger for cotton rhizosphere fungi than for rhizosphere bacteria. For many microbial groups, albeit still a small proportion of the entire microbiome, the variability among F₂ offspring is greater than random variability. These microbes include V. dahliae, Rhizoctonia solani, Alternaria alternate, and Fusarium, which are fungal pathogens of cotton, but also include well-known biocontrol agents of soilborne diseases, such as Trichoderma brevicompactum [39], Streptomyces [40], and Rhizobiales [41]. Thus, plant genotypes may affect the relative abundance of pathogens and beneficial microbes, influencing plant susceptibility and disease development.

There is significant genetic component in a number of PCs, which are determined by many microbes. This is most likely due to two reasons. Firstly, plant traits that are expected to affect microbiome composition and activity, such as root exudates [42] and root architecture [43], are complex in nature and may be controlled by many genes. Specific root exudates would selectively attract those microbes that could directly or indirectly use the exudates metabolically [44]. Root exudates and other root deposits secreted by host plants influence rhizosphere composition [45]. Not all root exudates are directly involved in plant nutrition and growth. Some of them act as signalling molecules to mediate interactions in root communities [46]. Secondly, many microbes may share similar environmental requirement for their development.

Although much progress in manipulating crop microbiomes has been made recently, further research is still needed to implement holobiont-level breeding effectively. Present results suggested that exploiting rhizosphere microbiome for sustainable agriculture via breeding could be difficult. This can summarised in one phrase: an one-to-many relationship. One selection criterion may select for or against many microbes (e.g., those contributing significantly to specific microbiome features as captured by PCs). These affected microbes may include beneficial and pathogenic microbes, as demonstrated by the present study. Moreover, many microbes that may be co-selected but with unknown identities and functions. Therefore, in the short term, soil amendment with specific microbial consortia may be the way forward to improve soil conditions and support intensive agriculture. This is similar to augmented application of biocontrol microbes to manage pathogens or inundative release of predator to control pests. To improve the persistence of these introduced beneficial microbes, research is needed to understand their ecological requirement and how they interact with resident phytobiome.

We characterized the rhizosphere composition of cotton F₂ offspring plants of two crosses and demonstrated that specific rhizosphere components are likely to be genetically controlled by plants when rhizosphere microbiome is characterized as PCs or individual top microbial groups were considered. More encouragingly, there are common microbiome components (i.e., PCs) and specific microbial groups with significant genetic component between the two crosses.

PCA: principal component analysis; PCs: principal components; V_T: variance among the 100 F₂ progeny genotypes; V_E: random variability; OTUs: operational taxonomic units

Acknowledgements

We would like to especially thank Dr. Greg Deakin (NIAB East Malling Research, East Malling, Kent, UK) for providing an abundance of helpful advice on my data analysis and writing, and Mr. Dezheng Zhang and Mrs. Caixia Wu (Institute of Cotton Research of Chinese Academy of Agricultural Sciences, Anyang, China) for their assistance in field management and sampling.

Funding

This work was supported by National Natural Science Foundation of China (grant number 31901938).

Availability of data and materials

The raw sequencing data reported in this paper are publicly available in the NCBI Sequence Read Archive (SRA) under the Bioproject number PRJNA756865.

Authors’ contributions

FW, HZ, and XX planned and designed the research and experiments. ZF, LZ and HF managed the field trial stations. FW and YZ collected samples. FW, CY, ZF and HF conducted the laboratory analyses. FW, and XX analyszed the data. FW and XX wrote the manuscript. FW and HZ acquired the funds for the study. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Authors' information

¹ State Key Laboratory of Cotton Biology, Institute of Cotton Research of Chinese Academy of Agricultural Sciences, Anyang 455000, China; ² School of Agricultural Sciences, Zhengzhou University, Zhengzhou 450001; ³ NIAB East Malling Research, East Malling, West Malling, Kent ME19 6BJ, UK.

Goh CH, Veliz Vallejos DF, Nicotra AB, Mathesius U. The impact of beneficial plant-associated microbes on plant phenotypic plasticity. J Chem Ecol. 2013;39:826–39.
Chapelle E, Mendes R, Bakker PAHM, Raaijmakers JM. Fungal invasion of the rhizosphere microbiome. ISME J. 2016;10:265–8.
Mendes R, Garbeva P, Raaijmakers JM. The rhizosphere microbiome: significance of plant beneficial, plant pathogenic, and human pathogenic microorganisms. FEMS Microbiol Rev. 2013;37:634–63.
Wei F, Feng H, Zhang D, Feng Z, Zhao L, Zhang Y, et al. Composition of rhizosphere microbial communities associated with healthy and Verticillium wilt diseased cotton plants. Front Microbiol. 2021;12:618169. doi:10.3389/fmicb.2021.618169.
Wei F, Zhao L, Xu X, Feng H, Shi Y, Deakin G, et al. Cultivar-dependent variation of the cotton rhizosphere and endosphere microbiome under field conditions. Front Plant Sci. 2019;10:1659. doi:10.3389/fpls.2019.01659.
Coleman-Derr D, Tringe SG. Building the crops of tomorrow: Advantages of symbiont-based approaches to improving abiotic stress tolerance. Front Microbiol. 2014;5:1–6.
Lu T, Ke M, Lavoie M, Jin Y, Fan X, Zhang Z, et al. Rhizosphere microorganisms can influence the timing of plant flowering. Microbiome. 2018;6:231.
Berendsen RL, Pieterse CMJ, Bakker PAHM. The rhizosphere microbiome and plant health. Trends Plant Sci. 2012;17:478–86.
Naylor D, Degraaf S, Purdom E, Coleman-Derr D. Drought and host selection influence bacterial community dynamics in the grass root microbiome. ISME J. 2017;11:2691–704.
Schlaeppi K, Dombrowski N, Oter RG, van Themaat EVL, Schulze-Lefert P. Quantitative divergence of the bacterial root microbiota in Arabidopsis thaliana relatives. Proc Natl Acad Sci. 2014;111:585–92.
Peiffer JA, Spor A, Koren O, Jin Z, Tringe SG, Dangl JL, et al. Diversity and heritability of the maize rhizosphere microbiome under field conditions. Proc Natl Acad Sci U S A. 2013;110:6548–53.
Edwards J, Johnson C, Santos-Medellín C, Lurie E, Podishetty NK, Bhatnagar S, et al. Structure, variation, and assembly of the root-associated microbiomes of rice. Proc Natl Acad Sci U S A. 2015;112:E911–20.
Mueller UG, Sachs JL. Engineering microbiomes to improve plant and animal health. Trends Microbiol. 2015;23:606–17.
Pérez-Jaramillo JE, Mendes R, Raaijmakers JM. Impact of plant domestication on rhizosphere microbiome assembly and functions. Plant Mol Biol. 2016;90:635–44.
Busby PE, Soman C, Wagner MR, Friesen ML, Kremer J, Bennett A, et al. Research priorities for harnessing plant microbiomes in sustainable agriculture. PLoS Biol. 2017;15:1–14.
Wei Z, Jousset A. Plant Breeding Goes Microbial. Trends Plant Sci. 2017;22:555–8.
Kroll S, Agler MT, Kemen E. Genomic dissection of host–microbe and microbe–microbe interactions for advanced plant breeding. Curr Opin Plant Biol. 2017;36:71–8.
Wille L, Messmer MM, Studer B, Hohmann P. Insights to plant–microbe interactions provide opportunities to improve resistance breeding against root diseases in grain legumes. Plant Cell Environ. 2019;42:20–40.
Mendes LW, Mendes R, Raaijmakers JM, Tsai SM. Breeding for soil-borne pathogen resistance impacts active rhizosphere microbiome of common bean. ISME J. 2018;12:3038–42.
Wagner MR, Roberts JH, Balint-Kurti P, Holland JB. Heterosis of leaf and rhizosphere microbiomes in field-grown maize. New Phytol. 2020;228:1055–69.
Lundberg DS, Rio TG, Tringe SG, Dangl JL, Mitchell-olds T, Wagner MR. Host genotype and age shape the leaf and root microbiomes of a wild perennial plant. Nat Commun 2016;7. doi:10.1038/ncomms12151.
Walters WA, Jin Z, Youngblut N, Wallace JG, Sutter J, Zhang W, et al. Large-scale replicated field study of maize rhizosphere identifies heritable microbes. Proc Natl Acad Sci U S A. 2018;115:7368–73.
Clouse KM, Wagner MR. Plant genetics as a tool for manipulating crop microbiomes: opportunities and challenges. Front Bioeng Biotechnol. 2021. doi:10.3389/fbioe.2021.567548.
Klosterman SJ, Atallah ZK, Vallad GE, Subbarao KV. Diversity, pathogenicity, and management of Verticillium species. Annu Rev Phytopathol. 2009;47:39–62.
FAO. World reference base for soil resources. 84 World Soil Resources Reports. International Society of Soil Science. Rome: International Soil Reference and Information Centre; 1998.
Herlemann DPR, Labrenz M, Jürgens K, Bertilsson S, Waniek JJ, Andersson AF. Transitions in bacterial communities along the 2000 km salinity gradient of the Baltic Sea. ISME J. 2011;5:1571–9.
White TJ, Bruns T, Lee S, Taylor J. Am-plification and direct sequencing of fungal ribosomal RNA genes for phylogenetics. PCR Protoc. 1990;315–22.
Deakin G, Tilston EL, Bennett J, Passey T, Harrison N, Fernández-Fernández F, et al. Spatial structuring of soil microbial communities in commercial apple orchards. Appl Soil Ecol. 2018;130. doi:10.1016/j.apsoil.2018.05.015.
Edgar RC. UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat Methods. 2013;10:996–8.
Kõljalg U, Nilsson RH, Abarenkov K, Tedersoo L, Taylor AFS, Bahram M, et al. Towards a unified paradigm for sequence-based identification of fungi. Mol Ecol. 2013;22:5271–7.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550. doi:10.1186/s13059-014-0550–8.
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11:R106.
Dixon P. VEGAN, a package of R functions for community ecology. J Veg Sci. 2003;14:927–30.
Ganguly A, Ghosh A. Quantitative diophantine approximation on affine subspaces. Math Zeitschrift. 2019;292:923–35.
Wei F, Passey T, Xu X. Amplicon-based metabarcoding reveals temporal response of soil microbial community to fumigation-derived products. Appl soil Ecol. 2016;103:83–92.
Liu H, Brettell LE, Qiu Z, Singh BK. Microbiome-mediated stress resistance in plants. Trends Plant Sci. 2020;25:733–43.
Ritpitakphong U, Falquet L, Vimoltust A, Berger A, Métraux JP, L’Haridon F. The microbiome of the leaf surface of Arabidopsis protects against a fungal pathogen. New Phytol. 2016;210:1033–43.
Lebeis SL, Paredes SH, Lundberg DS, Breakfield N, Gehring J, McDonald M, et al. Salicylic acid modulates colonization of the root microbiome by specific bacterial taxa. Science. 2015;349:860–4.
Shentu X, Zhan X, Ma Z, Yu X, Zhang C. Antifungal activity of metabolites of the endophytic fungus Trichoderma brevicompactum from garlic. Brazilian J Microbiol. 2014;45:248–54.
Niu G, Chater KF, Tian Y, Zhang J, Tan H. Specialised metabolites regulating antibiotic biosynthesis in Streptomyces spp. FEMS Microbiol Rev. 2016;40:554–73.
Erlacher A, Cernava T, Cardinale M, Soh J, Sensen CW, Grube M, et al. Rhizobiales as functional and endosymbiontic members in the lichen symbiosis of Lobaria pulmonaria L. Front Microbiol. 2015. doi:10.3389/fmicb.2015.00053.
Zhalnina K, Louie KB, Hao Z, Mansoori N, Da Rocha UN, Shi S, et al. Dynamic root exudate chemistry and microbial substrate preferences drive patterns in rhizosphere microbial community assembly. Nat Microbiol. 2018;3:470–80.
Saleem M, Law AD, Sahib MR, Pervaiz ZH, Zhang Q. Impact of root system architecture on rhizosphere and root microbiome. Rhizosphere. 2018;6:47–51.
Huang X, Chaparro JM, Reardon KF, Zhang R, Shen Q, Vivanco JM. Rhizosphere interactions: root exudates, microbes, and microbial communities. Botany. 2014;275:267–75.
Moe LA. Amino acids in the rhizosphere: from plants to microbes. Am J Bot. 2013;100:1692–705.
Olanrewaju OS, Ayangbenro AS, Glick BR, Babalola OO. Plant health: feedback effect of root exudates-rhizobiome interactions. Appl Microbiol Biot. 2019;103:1155–66.

Table 1 Summary of analysis of the first 100 principal components (PCs) for fungal community. These 100 PCs of the top 3,493 fungal OTUs (accounting for 99.99% of the total number of sequences) for the two crosses, testing whether the variance among 100 progeny samples (V_T) is greater than the random variability (V_E). The probability value was adjusted for multiple tests (Benjamini-Hochberg). *: the two parents also differed significantly.

PC	% variance explained by PC	L1×Z49		Z2×J11
PC	% variance explained by PC	P value: V_T > V_E	(V_T - V_E)/V_E	P value: V_T > V_E	(V_T - V_E)/V_E
1	4.390	< 0.001	0.935
2	2.964			0.011	0.538
3	2.796	< 0.001	0.727	< 0.001	0.660
4	2.240	< 0.001	0.810	< 0.001	0.644
17	0.876	< 0.001	0.785
18	0.848	< 0.001	0.703	< 0.001*	0.686
19	0.834	0.005	0.554	0.041	0.473
22	0.751	< 0.001	0.657	< 0.001*	0.740
25	0.703	0.036	0.475
26	0.692			0.007	0.559
29	0.659	0.014	0.517
30	0.655	< 0.001	0.657	0.002	0.599
31	0.634	< 0.001	0.665
36	0.586			0.022	0.508
51	0.498	0.036	0.470
61	0.456			0.026	0.496

    </tr>
</tbody>

Table 2 Summary of analysis of the first 100 principal components (PCs) for bacterial community. These 100 PCs of the top 3,493 bacterial OTUs (accounting for 99.99% of the total number of sequences) for the two crosses, testing whether the variance among 100 progeny samples (V_T) is greater than the random variability (V_E). The probability value was adjusted for multiple tests (Benjamini-Hochberg).

PC	% variance explained by PC	L1×Z49		Z2×J11
PC	% variance explained by PC	P value: V_T > V_E	(V_T- V_E)/V_E	P value: V_T > V_E	(V_T- V_E)/V_E
2	3.332	0.006	0.510
3	2.068	0.001	0.586
4	1.915	0.005	0.520	0.015	0.490
5	1.549	0.002	0.561	< 0.001	0.664
8				0.012	0.505
9				0.019	0.477
10	0.933	< 0.001	0.774	0.001	0.616
11	0.872	< 0.001	0.684	0.006	0.536
12	0.813	< 0.001	0.773	0.001	0.600
13	0.773	< 0.001	0.933	0.018	0.480
14	0.741	< 0.001	0.779
15	0.710	< 0.001	0.852	0.001	0.614
16	0.704	< 0.001	0.833	0.001	0.611
17	0.696	< 0.001	0.810
18	0.669	< 0.001	0.627	< 0.001	0.679
19	0.657	0.017	0.465
20	0.648	0.044	0.408	< 0.001	0.784
21	0.634	0.042	0.412	< 0.001	0.647
22	0.626	< 0.001	0.824
23	0.616	< 0.001	0.832	< 0.001	0.703
24	0.610	< 0.001	0.831	0.013	0.500
25	0.607	< 0.001	0.691
26	0.600	< 0.001	0.644
28	0.568	< 0.001	0.741	< 0.001	0.716
29	0.563	< 0.001	0.624	< 0.001	0.684
30	0.560	< 0.001	0.694
31				0.011	0.513
32	0.551	0.002	0.561
33	0.549	0.004	0.528
34	0.542	< 0.001	0.621
35	0.537	0.032	0.428
37	0.523	0.002	0.545
38	0.517	0.002	0.546
41	0.507	0.001	0.574
42				0.013	0.498
43	0.497	< 0.001	0.607	0.004	0.558
48	0.476	0.029	0.436
96	0.349	0.023	0.448

Table 3 Summary of analysis of the top fungal 100 fungal OTUs. Summary of analysis of the top fungal 100 OTUs (with the highest counts), testing whether the variance among 100 progeny samples (V_T) is greater than the random variability (V_E). The probability value was adjusted for multiple tests (Benjamini-Hochberg).

OUT_ID	Taxonomy	L1×Z49		Z2×J11
OUT_ID	Taxonomy	P value: V_T > V_E	(V_T- V_E)/V_E	P value: V_T > V_E	(V_T- V_E)/V_E
OTU106	Fusarium (g)	0.022	0.477	0.022	0.472
OTU12	Verticillium dahliae	0.005	0.550	0.046	0.424
OTU127	Trichoderma (g)	0.015	0.502
OTU1439	Fusarium (g)	0.046	0.414	0.004	0.559
OTU1481	Cryptococcus (g)	0.035	0.430
OTU1595	Microascaceae (f)	< 0.001	0.698	< 0.001	0.809
OTU21	Pseudogymnoascus pannorum	0.022	0.480
OTU239	Fungi (k)	0.035	0.431
OTU24	Phaeosphaeria fuckelii	< 0.001	0.668
OTU27	Trichoderma brevicompactum	0.026	0.463	0.034	0.445
OTU304	Fungi (k)	0.001	0.621
OTU32	Dothideomycetes (c)	0.031	0.449	0.007	0.527
OTU430	Mortierellales (o)	< 0.001	0.717
OTU445	Fungi (k)	0.032	0.442	< 0.001	0.701
OTU452	Alternaria alternata	0.031	0.449	0.007	0.532
OTU51	Corynespora cassiicola	< 0.001	0.684	0.022	0.474
OTU53	Fungi(k)	0.022	0.475	0.002	0.611
OTU55	Fungi(k)	< 0.001	0.788	0.003	0.57
OTU6	Fusarium (g)	0.035	0.432
OTU60	Mortierella alpina	0.031	0.452
OTU65	Neonectria (g)	0.001	0.607	0.005	0.546
OTU727	Cladosporium cladosporioides	0.010	0.522
OTU75	Mortierella alpina	0.020	0.486
OTU8	Fusarium solani	0.003	0.568
OTU81	Truncatella angustata	0.034	0.438
OTU82	Gymnoascus reesii	0.031	0.446
OTU83	Rhizoctonia solani	0.011	0.516	0.025	0.464
OTU88	Acremonium (g)	0.004	0.562	< 0.001	0.841
OTU91	Ascomycota (p)	0.002	0.593
OTU119	Hypocreales (o)			0.043	0.430
OTU1201	Ascomycota (p)			0.027	0.457
OTU1519	Metarhizium anisopliae			0.003	0.583
OTU22	Hannaella luteola			0.010	0.512
OTU28	Cercospora coniogrammes			0.003	0.570
OTU34	Dothideomycetes (c)			0.003	0.570
OTU37	Hannaella sinensis			0.009	0.516
OTU587	Chaetomium (g)			0.003	0.577
OTU61	Acremonium acutatum			0.002	0.604
OTU9	Ascomycota (p)			0.016	0.490

Table 4 Summary of analysis of the top 100 bacterial OTUs. Summary of analysis of the top bacterial 100 OTUs (with the highest counts), testing whether the variance among 100 progeny samples (V_T) is greater than the random variability (V_E). The probability value was adjusted for multiple tests (Benjamini-Hochberg).

OUT_ID	Taxonomy	L1×Z49		Z2×J11
OUT_ID	Taxonomy	P value: V_T > V_E	(V_T - V_E)/V_E	P value: V_T > V_E	(V_T - V_E)/V_E
OTU1396	Rhizobiales (o)	0.000	0.646
OTU187	Povalibacter uvarum	0.007	0.554
OTU24	Luteimonas cucumeris	0.031	0.489
OTU2479	Streptomyces (g)	0.000	0.728
OTU26	Ilumatobacter (g)	0.000	0.671
OTU51	Myxococcales (o)	0.017	0.518
OTU54	Steroidobacter denitrificans	0.000	0.698
OTU59	Actinomycetales (o)	0.004	0.576
OTU88	Catellatospora methionotrophica	0.000	0.653
OTU887	Bacillus niacini	0.039	0.475
OTU57	Achromobacter mucicolens			0.005	0.615

FigureS1.tif
Fig. S1 Rarefaction curves for the full study. Rarefication curves showing the accumulated sequence reads against the number of fungal (A) and bacterial (B) OTUs, which indicates that sequencing depth is sufficient for all samples.
TableS1.docx
Table S1 Wilt scores of individual plants in field.

Genetic Control of Rhizosphere Microbiome of Cotton Plant Under Field Conditions

Status:

Journal Publication

Version 1

Abstract

Figures

Background

Materials And Methods

Results

Discussion

Conclusions

Abbreviations

Declarations

References

Tables

Supplementary Files

Status:

Journal Publication

Version 1