GWAS, MWAS and mGWAS provide insights into precision agriculture based on genotype-dependent microbial effects in foxtail millet – Nature.com

GWAS identifies genetic variations associated with agronomic traits in foxtail millet

A total of 827 foxtail millet cultivars collected from China were sequenced and genotyped using common single-nucleotide polymorphisms (SNPs) based on a ~423Mb Setaria italica cv. Zhanggu reference genome (v.2.3)27. In total, 161,562 SNPs were detected after stringent steps of quality control, including population stratification and pedigree filtering, individual- and site-level call-rate filtering, and minor allele frequency (MAF) filtering. The SNPs were evenly distributed along chromosomes and the genetic distance for linkage disequilibrium (LD) decay to its half maximum was 9kb (Supplementary Fig.1A, B). Phylogenetic analysis based on the genetic SNPs revealed three main groups in the tested foxtail millet cultivars (Supplementary Fig.1C).

In addition, we planted these 827 foxtail millet cultivars for a field trial in Yangling, China, and measured their agronomic traits (Supplementary Data1). Twelve agronomic traits were used for further analysis, including six growth traits and six yield traits. The growth traits were mainly composed of top second leaf length (TSLL), top second leaf width (TSLW), main stem height (MSH), main stem width (MSW), panicle diameter of the main stem (MSPD) and fringe neck length (FNL) while the yield traits were represented by panicle length of the main stem (MSPL), per plant grain weight (PGW), main stem panicle weight (MSPW), hundred kernel weight (HKW), spikelet number of the main stem (MSSN) and grain number per spike (SGN). Genotypephenotype analysis showed that all 11 traits were significantly heritable except the trait HKW (H2=0.006, P=0.15). Growth traits exhibited higher heritability than yield traits, for example, MSPD showed the highest heritability (H2=0.46, the broad sense heritability) while PGW showed the lowest heritability (H2=0.16) (Supplementary Fig.2). GWAS on phenotypes was performed to identify the SNPs associated with the growth and yield traits. In total, 86 significant SNP loci and 91 associations for 10 traits (except MSPW and TSLW) were identified under suggestive P-value thresholds (P<2.01e5), some of which were for multiple traits (Fig.1, Supplementary Data2). Among these, 15, 16, 11, 10 and 16 significant SNPs co-located on chromosomes 2, 4, 6,7 and 9, respectively. The candidate genes located around the significant signal were analyzed for known molecular functions (Supplementary Data3). Firstly, several candidate genes responsible for growth and development regulation were observed such as ATG8C, ERF1B, PRR37 and Cyclin-like F-box. For example, the peak SNP signal si7:30050703 of MSW, located within the genic region of a homolog of ATG8C (autophagy-related protein 8C, Fig.1B), which functions in the early development of xylem and phloem tissues28. Additionally, SNP si2:6562955 was associated with MSW and near ERF1 (ethylene-responsive transcription factor) (Fig.1B and Supplementary Data3). ERF1 is implicated in cambium proliferation29, which might influence main stem width. Interestingly, the candidate gene PRR37 near the peak SNP si2:49328133 of the MSPL (Fig.1D), suppressed heading and showed shorter panicle length than its mutant in rice30, which might directly regulate the panicle length of foxtail millet. Besides, the peak SNP si2:6646016 of PGW was located within the genic region of Cyclin-like F-box (Fig.1C), which controls many crucial processes such as embryogenesis, hormonal responses, seedling development, floral organogenesis, senescence, and pathogen resistance31,32.

Manhattan plots showing the genome-wide associations between host genetic SNPs and A panicle diameter of the main stem (MSPD), B main stem width (MSW), C per plant grain weight (PGW) and D panicle length of the main stem (MSPL). The dotted line corresponds to a significance threshold of 2.01e5. Genes with significant SNPs are marked in red; genes near the significant SNPs are marked in green. NADKs: NAD+ kinase; PP2C:Phosphoinositide phospholipase C 2; WAT1:WAT1-related protein; ERF1: ethylene-responsive transcription factor 1; ATG8C:autophagy-related protein 8C; PRR37: two-component response regulator-like PRR37.

Secondly, numerous drought stress-responsive (PP2C, ARR12, NPF1.2, NPF4.6, WDR26, Plastocyanin-like protein, CPK2a, PIP5K1) and tolerant genes (APX, DTX12, bHLH3, Thioredoxin fold domain containing protein, SAPK9, Ca2+-transporting ATPase, InsP3, E3 ubiquitin-protein ligase, MIEL1) whose expression are frequently upregulated and contribute to drought resistance in drought-stressed seedlings, were found to be associated with the growth and yield traits (Supplementary Data3). For example, the SNP si2:49320133 that was associated with MSPD was located within the genetic region of PP2C (phosphatidylinositol-specific phospholipase C) (Fig.1A, Supplementary Data3), a stress and ABA-responsive gene that is involved in many physiological responses, including salt, drought and osmotic stress, carbon fixation in C4 plants, and inducible plant responses to pathogen33,34. In addition, NPF1.2 (protein NRT1/ PTR FAMILY 1.2) near the peak SNP (si4:27590764) of MSPD, which functioned as an ABA importer, is important for the regulation of stomatal aperture in inflorescence stems of Arabidopsis35. Another candidate gene SAPK9 (serine/threonine-protein kinase SAPK9) near the significant SNP (si3:44029863) of PGW improves drought tolerance and grain yield in rice by modulating cellular osmotic potential, stomatal closure and stress-responsive gene expression36 (Supplementary Data3).

Thirdly, a great number of plant immune responsive genes and pathogen defense genes were also found to be associated with traits, mainly including RPP13, RGA2, RPS2, LRR-RLKs, EF-Tu SYP22, NOG1, BBE, NB-ARC, and WAK2. Finally, several candidate genes responsible for nutrient uptakes such as iron transporter (IRT1, IRT2) and phosphate transporter (PT) were also observed (Supplementary Data3). Most of the candidate genes related to the significant SNPs were mainly involved in abiotic and biotic stress responses, implying that the host genotype and environment interaction might co-contribute to plant adaption and modulate the traits of foxtail millet.

To explore the contributions of genetic variations to plant performance, linear regression models were used to calculate the role of host genotypes on key growth (TSLW, MSPD, MSW)- and yield (MSPW, PGW, MSPL)-related traits of the 827 different foxtail millet cultivars. Considering no SNPs associations with phenotypes TSLW and MSPW under suggestive thresholds, we extended the candidate SNPs (adjusted P<1.0e4) as inputs of the linear regression models34,35 (Supplementary Data4). After performing thirty rounds of five-fold cross-validation, the genetic SNP markers in predicting model could explain an average of 32.82%, 28.55%, 47.27%, 15.02%, 38.89% and 64.60% of the variances in TSLW, MSPD, MSW, MSPW, PGW and MSPL in the testing data, respectively (Supplementary Fig.3).

Root associated microbiota are thought to promote resistance to pathogens and tolerance to specific environmental constraints, and also contribute to plant performance37,38. Firstly, linear regression models were performed to calculate the effect of the rhizoplane microbiota on growth- and yield-related traits of the 827 different foxtail millet cultivars. The 1004 rhizoplane operational taxonomic units (OTUs) with a 70% occurrence in all samples (here defined as common OTUs), covering an average of 61.30% of total abundances were used as the input data as these OTUs commonly exist in the root zone of foxtail millet cultivars. The common sub-community (1004 common OTUs) showed higher evenness and correlations with the growth traits than the whole microbial community (Kruskal-Wallis test (one-way analysis), Pevenness=2.58e30, Supplementary Fig.4AC). The average variation degree (AVD) index from the common sub-community, 0.5 and 0.3 sub-community (OTUs with 50 and 30% occurrence), were calculated to assess the microbiota stability. The common sub-community had a lower AVD value than the other two sub-communities, indicating that it has a more stable microbiota (Supplementary Fig.4D). In the common sub-communities, the moderate OTUs (covered 83.67% of OTU numbers) were abundant, followed by abundant OTUs (12.85%) and rare OTUs (3.48%) (Supplementary Table1). The network analysis was used to disentangle the ecological role and co-occurrence patterns of 1004 OTUs in the common sub-community. Abundant OTUs (ATs) had significantly higher values of the degree, closeness, betweenness centrality and hub scores than both rare (RTs) and moderate OTUs (MTs) in the network (Kruskal-Wallis test (one-way) with P<0.001, Supplementary Fig.4E), indicating their important roles in sustaining the stability of the microbial community. Thus, the candidate OTUs that were significantly correlated with the traits (adjusted P<0.05) were selected from the common sub-community and used as the input of the predicting models (Supplementary Data5). A five-fold cross-validation approach was repeated thirty times for each trait to reduce the noise in the estimated model performance. The candidate OTU markers in the OTU-predicting models explained an average of 32.47%, 17.43%, 56.06%, 30.36%, 35.17% and 12.61% of the variances in TSLW, MSPD, MSW, MSPW, PGW and MSPL in the testing data, respectively (Supplementary Fig.3).

To explore the contributions of genetic variations and environmental microbiota to plant performance, we used linear mixed models to calculate the role of host genotypes and rhizoplane microbiota on the aforementioned growth and yield traits. We used the candidate SNPs (adjusted P<1.0e4) as inputs of the linear regression models39,40 to predict phenotypic variations (Supplementary Data4). Then the candidate OTUs markers from the above models were also added to the linear regression model (Supplementary Data5). After performing thirty rounds of five-fold cross-validation, the genetic SNP and OTU markers in predicting model could explain an average of 46.50%, 59.08%, 65.69%, 38.45%, 43.04% and 44.31% of the variances in TSLW, MSPD, MSW, MSPW, PGW and MSPL in the testing data, respectively (Supplementary Fig.3). The correlation coefficients only using genotype as variables were obviously higher than that only using root microbiota as variables in several agronomic traits such as MSW and MSPL.However, in the trait MSPD and MSPW, the contribution of root microbiota to phenotypic plasticity was higher when root microbiota variables were used instead of genotype variables alone, indicating a different contribution of host genotype and root microbiota to phenotypic plasticity. The combination of host genotype and root microbiota significantly promoted the explanation of variations in all six traits than genotype and root microbiota alone (Wilcox rank test, P<0.001) except for the trait MSPL (Supplementary Fig.3). Consistently, the panicle length has been proven to be directly impacted by gene PRR3725, similar to our observed data. The predictive models with the best prediction accuracy for the phenotypes using the SNP and OTU variables were retained, which explained 53.42%, 63.73%, 70.54%, 50.16%, 55.88%, and 54.82% variations for TSLW, MSPD, MSW, MSPW, PGW and MSPL trait, respectively, resulting in a final set of 257 marker OTUs (Fig.2AF, Supplementary Data5). Network analysis of 257 marker OTUs showed that the abundant marker OTUs (AMTs) had a significantly higher value of the degree, closeness and betweenness centrality than both rare (RMTs) and moderate marker OTUs (MMTs), indicating the abundant marker OTUs have more important roles in community structure (Kruskal-Wallis test (one-way) with P<0.05, Supplementary Fig.5 and Supplementary Table1).

AF The variation of growth (TSLW, MSW, MSPD) and yield (MSPW, PGW, MSPL) traits explained by the genetic SNPs and microbial OTUs combined. Each panel shows observed values on the x-axis and model-predicted values on the y axis, with a fitted linear regression. Specifically, the predicted value of TSLW, MSW, MSPD, MSPW, PGW, and MSPL is calculated based on 136, 100, 117, 126, 110 and 106 samples in the testing dataset, respectively. The dark trend line illustrates the predicted effect in the linear model (LM). The gray shading around the line represents a confidence interval of 0.95. TSLW, top second leaf width; MSW, main stem width; MSPD, panicle diameter of the main stem; MSPW, main stem panicle weight; PGW, per plant grain weight; MSPL, panicle length of the main stem.

Among the 257 marker OTUs identified by MWAS, 145 and 128 marker OTUs were significantly correlated with growth and yield traits, respectively (Supplementary Data5). Taxonomic profiling of these marker OTUs revealed 86 genera distributed across 15 phyla. The top five abundant phyla were Proteobacteria (with 68 OTUs), Actinobacteria (54 OTUs), Bacteroidetes (36 OTUs), Acidobacteria (35 OTUs), and Firmicutes (33 OTUs) (Fig.3A). In particular,17 marker OTUs were shared by growth and yield traits. Unexpectedly, no marker OTU or genus was shared by all six traits (Supplementary Fig.6), suggesting that the microbial markers may function in different development stages or different processes of foxtail millet.

A Phylogenetic tree of the 257 microbial markers associated with agronomic traits of foxtail millet. The outer circle represents the phylum level. The beta estimates of the microbial OTUs to growth and yield traits are plotted in the inner circles, respectively. The arrows indicate the strains tested in planta (B, C), including strains responded to six positive marker OTUs: Acid550 to Acidovorax OTU_46, Baci299 to Bacillaceae OTU_22228, Kita594 to Kitasatospora OTU_8, Baci154 to Bacillus OTU_19414, Baci312 to Bacillus OTU_25704, Baci429 to Bacillales OTU_381, and strains responded to four negative marker OTUs: Shin228 to Shinella OTU_37, Baci81 to Bacillus OTU_54, Baci173 to Bacillaceae OTU_19835 and Baci554 to Bacillaceae OTU_28133. The strains predicted to affect growth traits are validated by plate (B) and sterilized soil (C). Significance is determined within each pair of treatment and control via one-tailed t-test and the P-values are adjusted by Benjamini-Hochberg (BH) method. n=41, 40, 38, 34, 43, 43, 44, 32, 21, 34 and 27 (from left to right) biological replicates in plate experiment. From Kita594 to Baci554, adjusted P(plant height)=8.68e07, 1.02e09, 0.07, 0.001, 0.49, 0.29, 0.03, 0.40, 0.001, 5.39e06, adjusted P(root length)=0.012, 0.10, 0.40, 0.18, 0.01, 0.11, 1.14e06, 0.02, 1.20e09, 5.68e12. n=20, 32, 31, 40, 23, 20 and 23 (from left to right) biological replicates in sterilized soil experiment. From Kita594 to Baci554, adjusted P(plant height)=1.6e07, 0.003, 0.018, 0.011, 0.98, 0.06, adjusted P(root length)=0.046, 4.4e04, 8.43e05, 0.016, 0.058, 0.15. *, ** and *** represented the adjusted P<0.05, 0.01 and 0.001, respectively. The box depicts the interquartile range (IQR) between the 25th and 75th percentiles, respectively and the line within the box represents the median. The whiskers extend 1.5 times the IQR from the top and bottom of the box, respectively.

To validate the predicted effects of these microbial markers on foxtail millet growth, we isolated a range of taxonomically different bacterial strains from root microbiota of the foxtail millet varieties grown in the field. A total of 644 bacterial strains were collected, and 257 bacterial isolates with complete 16S rRNA gene sequences were retained, representing four bacterial phyla and 25 genera (Supplementary Data6).

A cultured strain was considered a representative OTU if its 16S rRNA gene had 97% similarity with the rhizoplane microbiota OTU (Supplementary Data6). Representative cultivated strains of six positive marker OTUs (Acid550 to Acidovorax OTU_46, Baci299 to Bacillaceae OTU_22228, Kita594 to Kitasatospora OTU_8, Baci154 to Bacillus OTU_19414, Baci312 to Bacillus OTU_25704 and Baci429 to Bacillales OTU_381) and four negative marker OTUs (Shin228 to Shinella OTU_37, Baci81 to Bacillus OTU_54, Baci173 to Bacillaceae OTU_19835 and Baci554 to Bacillaceae OTU_28133) with top beta estimation in the regression model were selected for the validation experiments (Fig.3A and Supplementary Fig.7). We co-cultivated these 10 biomarker strains with foxtail millet Huagu12 (a bred cultivar of foxtail millet (Setaria italica) at Shenzhen, China) for 7-days in sterilized plates, and observed altered root lengths and plant heights compared with the control (Fig.3B and Supplementary Fig.7A). The positive biomarker strains representing OTUs with top beta estimation showed significant growth-promoting abilities. Specifically, positive biomarker strain Kita594 (Kitasatospora OTU_8) promoted both root and stem growth, whereas Baci299 (Bacillus OTU_22228) and Acid550 (Acidovorax OTU_46) only promoted shoot growth compared to the control (one-tailed t-test with adjusted P<0.05, Fig.3B and Supplementary Fig.7A). The negative marker strain Baci173 (Bacillaceae OTU_19835) and Baci554 (Bacillaceae OTU_28133) suppressed the shoot and root growth of Huagu12 (one-tailed t-test with adjusted P<0.05, Fig.3B and Supplementary Fig.7A). While the negative marker strains Shin228 (Shinella OTU 37) and Baci81 (Bacillus OTU 54) exhibited growth-promoting effects, they may only function in special root microbial flora in collaboration with other strains or be mistakenly identified as representative strains due to high 16S rDNA sequence similarities with negative marker OTU 37 and 54.

Next, we validated the effects of four positive marker strains (Kita594, Baci299, Baci154 and Acid550) and two negative marker strains (Baci173 and Baci554) with good promoting or suppressing performances on plant growth in plate experiment by watering millet seedlings grown in sterilized soil with these bacterial suspensions separately. Consistently, the seedlings watered with suspensions of the promoting bacterial strains Kita594, Baci299, Baci154 and Acid550 showed significantly increased plant height and root length compared with the control, whereas the seedlings watered with suspensions of the suppressing bacteria Baci173 showed shorter roots (one-tailed t-test with adjusted P<0.05, Fig.3C and Supplementary Fig.7B). These results validated the plant growth promoting (PGP) traits of marker microbes in foxtail millet.

To shed light on how bacterium regulates the growth of foxtail millet, we analyzed the transcriptomes of seedlings colonized for 14 days with the growth-promoting strains Baci299, Acid550, Kita594 or with the growth-suppressing strain Baci173. The differentially expressed genes from biomarker strain-inoculated versus non-inoculated samples were enriched in different pathways (Fishers exact test, q<0.05, Fig.4A). For example, the differentially expressed genes caused by growth-promoting strains were mainly enriched in the pathways such as Phenylalanine, tyrosine and tryptophan biosynthesis (ko00400), Biosynthesis of amino acids (ko01230), Phenylalanine metabolism (ko00360), Carbon fixation in photosynthetic organisms (ko00710), Photosynthesis-antenna proteins (ko00196), Photosynthesis (ko00195), MAPK signaling pathway-plant (ko04016), Plant-pathogen interaction (ko04626), Diterpenoid biosynthesis (ko00904), Monoterpenoid biosynthesis (ko00902), alphaLinolenic acid metabolism (ko00592) and Selenocompound metabolism(ko00450), while the differentially expressed genes caused by suppressing strain were mainly involved in the pathways such as Arginine and proline metabolism (ko00330) and Valine, leucine and isoleucine degradation (ko00280) (Fig.4A).

A KEGG enrichment analysis of differentially expressed genes in Baci173-, Baci299-, Acid550-, Kita594- inoculated seedlings. The differentially expressed genes represent the genes that were significantly upregulated or downregulated in seedlings inoculated with marker strain compared with control. Circle size represents the number of genes within the pathway and color represents the significance of the pathway. B Venn diagram showing the overlap of the significantly upregulated genes under different inoculations. C Transcript abundance of genes that were induced only in Baci173-, Baci299-, Acid550- and Kita594-inoculated seedlings, respectively.

Interestingly, the growth-promoting strains displayed strain-specific induction of genes involved in nutrient transformation, pathogen defense, anti-abiotic stresses and growth-promoting processes (Fig.4B, C, Supplementary Data7). For instance, the ammonia producing gene (K01455_fomamidase) and terpenoids synthase (K15803, germacrene D) were highly induced by strain 299; ethylene synthase (K05933, aminocyclopropanecarboxylate oxidase) and plant immunity responsive genes (K18834, WRKY1; K20538, MPK8; K00430, peroxidase; K13422, MYC2; and K04079, HSP90A) were abundantly induced by strain 550, and photosynthesis-related genes (K02692, psaD, K01092, IMPA; and K08916, LHCB5), anti-oxidant gene (K00434, Lascorbate peroxidase) and pterostilbene biosynthesis gene (K16040, ROMT) were highly induced by strain 594 (Fig.4C). Intriguingly, the expansin gene that mediates cell wall loosening and increased root and shoot growth in rice41, was induced by all of the growth-promoting strains.

Similarly, 39 genes were significantly induced only by the growth-suppressing strain Baci173, including auxin synthetase (K01426, amidase; K00128, aldehyde dehydrogenases ALDH), auxin-responsive protein IAA (K14484, auxin-responsive protein), L-glutamine synthetase (K01915) and branched-chain amino acid synthetase (BCAT, K00826) (Fig.4B, C, Supplementary Data7), which all have well-documented roles in inhibiting root growth42,43,44. Thus, the plant growth mechanisms mediated by microorganisms were strain-dependent.

To explore the relationship between the host genotype and rhizoplane microbial composition, Mantels test was used to evaluate the correlation between host phylogenetic distances and rhizoplane microbiota distance, exhibiting a significant Mantels correlations (r=0.06, P=0.0003, 9999 permutations). Subsequently, to investigate host genotype-dependent variation in the foxtail millet rhizoplane microbiota, the heritable microbes were identified based on a common rhizoplane OTUs data set, which covered 17 phylum and 52 orders. Using an SNP-based approach, the heritability for individual OTU was calculated. 281 OTUs with H2 (the broad sense heritability) more than 0.15 were defined as highly heritable and the others as lowly heritable (Supplementary Data8). Bacillales and Gp4 orders enriched greater numbers of highly heritable OTUs when compared with the lowly heritable fraction (Fishers exact test, q<0.05, Supplementary Fig.8A), implying that these bacterial orders were more easily impacted by host genotypes of foxtail millet. To explore whether there are similarities in heritable microbes across Poaceae family, we compared the top 100 most heritable OTUs from foxtail millet, sorghum45 and maize datasets46,47. After removing the order with a total number of OTUs less than 4, seven bacterial orders such as Bacillales, Actinomycetales, Burkholderiales, Rhizobiales, Myxococcales, Sphingobacteriales and Xanthomonadales were identified, which shared and covered more than half of the most heritable OTUs from foxtail millet, sorghum, and maize datasets, respectively (Supplementary Fig.8B, C). These results hence indicated that the microorganisms in these bacterial orders were more sensitive to genetic variations across both sorghum, maize and foxtail millet.

To further assess the association of host genetic variations and root microbial abundance, we ran mGWAS on 1004 common rhizoplane OTUs of foxtail millet. We identified significant associations of 2108 SNP loci with 838 microbial OTUs (here called SNPs-associated OTUs) at the genome-wide suggestive significance threshold of P<2.01e5 (Supplementary Data9). To identify how the host genetic variations drove abundance variations of the specific microbial taxonomies, especially the bacterial orders that were more sensitive to genetic variations, the SNP-associated genes for each order were enriched into pathways (Supplementary Fig.9). However, only four bacterial orders associated genes were significantly enriched into different pathways. Taking Bacillales for example, the associated genes were mainly enriched in the monoterpenoid biosynthesis pathway (Fishers exact test, q=0.05, Supplementary Fig.9). The GP4 associated genes were significantly enriched in producing D-galacturonic acid (Fishers exact test, q=0.08, Supplementary Fig.9). GP4 from Acidobacteria phylum, which has been reported with the capability of utilizing galacturonic acid, a characteristic component of the cell wall in higher plant48, might be recruited to rhizoplane by plant-secreted galacturonic acid. The genes associated with plant pathogen-containing order Xanthomonadales were significantly enriched into the pathway such as peroxisome and MAPK signaling pathway (Fishers exact test, q=0.01 and 0.04, Supplementary Fig.9), which are involved in disease and abiotic resistance49,50. These results provide key insights into how the host genetic mechanism drive plant-associated microbiota.

In addition, significant SNP loci located in the generic region were also deeply analyzed (Fig.5, Supplementary Data9). For example, the peak SNP signal si7:13687399 located in the genic region of bHLH35 was associated with 39 common OTUs from different microbial taxonomies such as Acidobacteria (28), Proteobacteria (8) and Bacteroidetes (3). bHLH35 proteins are transcription factors induced by effector-triggered immunity (ETI), and also involved in tolerance to abiotic stresses51. The SNP si1:32157654 located in the generic region of WAK2 (wall-associated receptor kinase 2) was associated with 30 common OTUs, including Acidobacteria (21), Bacteroidetes (4), Proteobacteria (4) and Actinobacteria (1). The WAK2 protein bound to pectin, is required for cell expansion and is induced by a variety of environmental stimuli, including pathogens and wounding52. Similarly, the 50 common OTUs were found to be associated with FLS2 (si7:2994337, Supplementary Data9), a flagellin sensor that perceives conserved microbial-associated molecular patterns (MAMPs) in the extracellular environment53. Clostridia OTU_19207 and Nocardioides OTU_26357 associated si8:20598566 located within the gene of NPF1.2 (Fig.5 and Supplementary Data9), which is involved in ABA importing and nitrate utilization, regulates plant development and influences the root microbiota14,35,54. An NPF1.2 homologue in loci si1:20064466 was significantly associated with Bacillaceae OTU_28839. Collectively, host genes related to plant immunity (RPM1, RGA2, HSL1, CRKs, LRR-RLKs), metabolites (Flavonoids, Diterpenes, amyA, alpha-N-arabinofuranosidase, beta-glucuronosyltransferase), nutrient uptake (Acid phosphatase, Mg2+ transporter, H+-transporting ATPase), plant hormone signal transduction (BRI1, DELLA protein, EFR3, PI-PLC, SDR, ARR1) and others (E3 ubiquitin protein ligase) are perhaps common host genetic factors with function to modulate root microbial composition assembly (Fig.5 and Supplementary Data9).

Manhattan plots show the significant SNPs for microbial abundance. SNPs located in gene coding regions are labeled with numbers. Details of the associations between the host genes and microbial species are given in the table below. All of these associations of SNP loci and microbial OTUs were significantly lesser than 2.01e5.

Plants primarily influence their microbiomes through targeted interactions with key taxonomic groups or diffuse interactions with entire communities55. To further investigate the mode of host-microbe interactions, the hub microbial taxa and non-hub microbial taxa and their associated genes were identified. Firstly, we defined hub taxa as OTU with high values of degree (>400) and closeness centrality (>0.5) in the network as described in a previous study56, resulting in 102 hub OTUs. We identified that 90 hub OTUs and 748 non-hub OTUs had significant associations with the host genetic SNP loci (Supplementary Fig.10A, Supplementary Data9), indicating host plant might interact with these hub microbes and diffusely interact with these non-hub microbes. We aggregated these SNP-associated hub OTUs (90 hub OTUs) and non-hub OTUs (748 OTUs) into 12 and 36 microbial orders, respectively. Comparative analysis showed that one order GP7 was only composed of SNP-associated hub OTUs, and 25 orders such as Sphingobacteriales, Bacillales, Ohtaekwangia, Sphingomonadales and Acidimicrobiales were only composed of SNP-associated non-hub OTUs, and 11 orders were composed of both SNP-associated hub and non-hub OTUs (Supplementary Fig.10B). These data indicated that the foxtail millet employed two modes to structure the rhizoplane microbiota: targeted interaction with several hub microbes and diffused interaction with most of the microbes. To decipher the potential mechanism of the interaction between plant and microbe, the candidate host genes around the SNP loci associated with the hub and non-hub OTUs were extracted separately. The networks showed that the host immune genes FLS2 and transcription factor bHLH35 are widely associated with the hub and non-hub taxa (Supplementary Fig.11A, B). However, the host plant still employed different genes to interact with different taxa (Supplementary Fig.11C), suggesting a taxa-dependent regulation model.

To determine if the genotype-dependent rhizoplane microbiota influence agronomic traits in foxtail millet, we compared the 838 SNP-associated OTUs (mGWAS identified) with the 257 marker OTUs (MWAS identified). We discovered that 219 of the SNP-associated OTUs overlapped with the marker OTUs in our data sets, covering 85.2% of 257 marker OTUs. (Supplementary Fig.12A, 219 out of 257=85.2%). 682 SNP loci were significantly associated with 219 marker OTUs (here called marker OTU-associated SNPs). However, for the 682 marker OTU-associated SNPs, only 4 overlapped with the 45 non-redundant marker SNPs (GWAS identified) that were associated with the aforementioned agronomic traits of foxtail millet (Supplementary Fig.12B). Most of the genetic variations that were associated with marker OTUs were not directly associated with agronomic phenotypes. These genetic variations might affect agronomic phenotypes indirectly, only in the presence of environmental factors such as marker microbes. Moreover, the Mantel test also showed that SNP-associated marker OTUs had higher correlations with the growth trait (MSPD and MSW) than non SNP-associated marker OTUs, while having no difference in correlations with trait TSLW, MSPW, PGW and MSPL (Supplementary Table2). It means that the genotype-dependent marker OTUs might explain more variances in plant growth traits.

To decipher host plant genetic mechanisms for marker microbe selection, KEGG pathway enrichment analysis revealed that the genes within or nearby the significant SNP loci were enriched in pathways related to plant-pathogen interaction (ko04626), MAPK signaling (ko04016), Steroid biosynthesis (ko00100) and so on (Supplementary Data10). Specifically, the genes enriched in plant-pathogen interactions included microbial pattern-recognition receptors (PRRs), disease-resistant genes RPM1 and RPS2, an activator of pathogenesis-related genes PTI1 and PTI6, key regulators of plant immune responses CALM and transcription factor WRKY25. These results suggest that the plant defense genes may also underpin the microbial ecology in the root habitat in addition to protecting from pathogens.

Among the 219 SNP-associated marker OTUs, 77 were highly heritable (Fig.6A). The association between host genetic variation, the abundance of specific marker microbes and phenotypes, especially for 77 genomic heritable marker OTUs were closely examined (Fig.6A, Supplementary Data8). Remarkably, plant defense-related genes and transcription factors, such as the plant immune receptor FLS2 (si7:2994337), transcription factor bHLH35 (si7:13687399) and WAK2 (si:2:5642650) had a dominant impact on the marker OTUs from the phylum of Acidobacteria (Fig.6A). In contrast, genes involved in nutrient uptake, metabolites and abiotic stress response, such as magnesium transporter (si7:19232862), triterpene synthase (si:7:11346839, achilleol B synthase), BGLU12 (si:3:4780643, Beta-glucosidase 12) and RXW8 (si:3:39749463, CSC1-like protein RXW8), mainly associated with marker OTUs from Actinobacteria, Bacteroidetes and Proteobacteria, which mostly contribute positively to the growth and yield traits of foxtail millet (Fig.6A). Other genes involved in plant growth and development processes, such as SUZ12 (Polycomb protein SUZ12) and WAT1 (WAT1-related protein), impacted the marker OTUs from Firmicutes (Bacillaceae OTU_19835) and Proteobacteria (Xanthomonadaceae OTU_10146) respectively, but these marker OTUs have opposite effects on the growth of foxtail millet (Fig.6A). Additionally, we observed strong associations between the positive marker Acidovorax OTU_46 and EREBP-like factor (si7:27291504, dehydration-responsive element-binding protein 1B-like), and between positive marker Kitasatospora OTU_8 and FaQR (si2:36177507, 2-methylene-furan-3-one reductase) (Fig.6A). To explore the host genetic mechanisms that might drive the associations of the plant host gene and rhizoplane microbiota, we examined the specific expression pattern of candidate genes from the RNA-seq datasets obtained from the sterilized soil experiments. Obviously, the genes FaQR, vWA (von Willebrand factor, type A), SUZ12 and EREBP-like factor (ethylene response element binding protein) exhibited significant variation after being inoculated with strain Kita594 (Kitasatospora OTU_8), Baci299 (Bacillaceae OTU_22228), Baci173 (Bacillaceae OTU_19835) and Acid550 (Acidovorax OTU_46) compare to control, respectively (Supplementary Fig.13), implying that the candidate host genes likely interacted with specific bacterial strains.

A Venn diagram displaying the overlaps among 838 SNP associated OTUs, 257 marker OTUs and 281 highly heritable OTUs. B A network of associations between the candidate genes and marker microbial OTUs. Edges between the marker OTUs and host genes were colored according to the correlation coefficients. The pink color represents the positive correlations while the green color represents the negative correlations. The circle represents the OTUs colored according to the phylum taxonomy information, the triangle represents the genes colored according to the function module information, the square represents the growth traits colored in green and the yield traits are shown in yellow. Plant height (C) and root length (D) of seedlings of FaQR reference cultivars (C494 and C1631) and allele cultivars (C571 and C1119) grown axenically (no bacteria, control) or with growth-promoting Kita594. n=22, 26, 36, 37, 50, 46, 45 and 41(from left to right) biological replicates. Plant height (F) and root length (G) of seedlings of SUZ12 reference cultivar (C946 and C306) and allele cultivar (C1296 and C1021) grown axenically (no bacteria, control) or with growth-suppressing Baci173. n=39, 32, 46, 40, 50, 24, 45 and 31 (from left to right) biological replicates. The deviation of promoting and suppressing effect of marker strain Kita594 (E) and Bci173 (H) were calculated separately. n=26, 37, 46, 41, 26, 37, 46 and 41 (from left to right) biological replicates for the treatment with marker strain Kita594 (E). n=32, 40, 24, 31, 32, 40, 24 and 31 (from left to right) biological replicates for the treatment with marker strain Bci173 (H). Different letters in C to H indicate statistical significance (adjusted P<0.05) among the treatments according to one-way ANOVA and LSD test at the 5% level. In C, df=7, F=19.73, adjusted P<2.0e16; in D, df=7, F=13.76, adjusted P=2.09e-15; in E, df=3, F=2.10, adjusted P=0.102; df=3, F=18.29, adjusted P=4.04e10; In F, df=7, F=18.11, adjusted P<2.0e16; in G df=7, F=14.83, adjusted P<2.0e16; in H, df=3, F=19.57, adjusted P=2.0e10; df=3, F=6.751, adjusted P=2.90e4. The box edges depict the 75th and 25th percentiles, respectively and the line within the box represents the median. The whiskers extend 1.5 times the IQR from the top and bottom of the box, respectively.

Finally, based on cultivars with different genotypes, the influence of functional SNPs on marker OTU abundance was thoroughly examined. The abundance of marker OTUs shifted among the different genotypes at the most strongly associated SNPs (Supplementary Fig.14). We hypothesize that host gene-regulated promotion/suppression microbes could establish genotype-dependent microbe-mediated growth phenotypes. To test this hypothesis, we germinated the FaQR and SUZ12 reference and allele foxtail millet cultivars on sterile plates inoculated with a growth-promoting or suppressing strain that corresponds to each cultivar: the growth-promoting strain Kita594 to FaQR reference (C494 and C1631) and allele (C1119 and C571) genotype cultivars, the growth-suppressing strain Baci173 to SUZ12 reference (C946 and C306) and allele (C1021 and C1296) genotype cultivars. Intriguingly, we found that strain Kita594 had a statistically significantly shoot-promoting effect only on the allele cultivars, but not on reference cultivars (Fig.6CE, adjusted P<0.05 by ANOVA-LSD), supporting that plant-growth promoting rhizobacteria support genotype-dependent cooperation with the plant. We observed strong root growth inhibition in seedlings inoculated with the growth-suppressing strain Baci173 (Fig.6FH, adjusted P<0.05 by ANOVA-LSD), and a more significant suppressing effect on root length was observed in the allele cultivars (C1296 and C1021) compared to the reference cultivars (C946 and C306). Significant effects of the interaction between the genotype and strain Kita594 and strain Baci173 on the shoot and root length were also detected by PERMANOVA, respectively (genotypes*Kita594: R2=13.048, P<0.001; genotypes*Baci173: R2=0.07, P<0.001, Supplementary Table3). Together, these results suggest that host genetic variation might impact the interactions between marker strains and host plants, finally affecting the plant phenotypes.

The rest is here:
GWAS, MWAS and mGWAS provide insights into precision agriculture based on genotype-dependent microbial effects in foxtail millet - Nature.com