Nitrogen- and phosphorus-starved Triticum aestivum show distinct belowground microbiome profiles

Many plants have natural partnerships with microbes that can boost their nitrogen (N) and/or phosphorus (P) acquisition. To assess whether wheat may have undiscovered associations of these types, we tested if N/P-starved Triticum aestivum show microbiome profiles that are simultaneously different from those of N/P-amended plants and those of their own bulk soils. The bacterial and fungal communities of root, rhizosphere, and bulk soil samples from the Historical Dryland Plots (Lethbridge, Canada), which hold T. aestivum that is grown both under N/P fertilization and in conditions of extreme N/P-starvation, were taxonomically described and compared (bacterial 16S rRNA genes and fungal Internal Transcribed Spacers—ITS). As the list may include novel N- and/or P-providing wheat partners, we then identified all the operational taxonomic units (OTUs) that were proportionally enriched in one or more of the nutrient starvation- and plant-specific communities. These analyses revealed: a) distinct N-starvation root and rhizosphere bacterial communities that were proportionally enriched, among others, in OTUs belonging to families Enterobacteriaceae, Chitinophagaceae, Comamonadaceae, Caulobacteraceae, Cytophagaceae, Streptomycetaceae, b) distinct N-starvation root fungal communities that were proportionally enriched in OTUs belonging to taxa Lulworthia, Sordariomycetes, Apodus, Conocybe, Ascomycota, Crocicreas, c) a distinct P-starvation rhizosphere bacterial community that was proportionally enriched in an OTU belonging to genus Agrobacterium, and d) a distinct P-starvation root fungal community that was proportionally enriched in OTUs belonging to genera Parastagonospora and Phaeosphaeriopsis. Our study might have exposed wheat-microbe connections that can form the basis of novel complementary yield-boosting tools.


Introduction
The spread of chemical fertilization practices was one of the main features of the Green Revolution. These yield-boosting methods are so efficient that they quickly became worldwide a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 pillars of intensive agriculture. The environmental costs associated with the currently generalized use of inorganic N and P fertilizers are, however, mounting [1][2][3][4]. Although it remains manageable now, this imbalance could eventually compromise the sustainability of our farming system [5][6]. It therefore constitutes a strong incentive for the development of complementary yield-boosting tools.
Farmers already have access to many such tools, including crop cultivars that possess improved nutrient use efficiencies, alternative agronomic practices, and microbial inoculants [7][8][9][10]. However, although appreciable, the reduction of global environmental impacts that these tools can collectively provide appears limited [9,11]. More will likely be needed to address the future challenges of sustainable food production. As many of the currently deployed complementary yield-boosting tools rest on natural N-or P-sharing plant-microbe partnerships, discovering new partnerships of these types could help us develop the necessary tools.
Targeting wheat may be a good choice for research endeavors that have such objectives. Several studies suggest that the plant can partner with N-and P-providing microorganisms [12][13][14][15][16][17][18][19][20], but its natural microbial associations are not well delineated. Since wheat farming currently consumes approximately 20% of the worldwide production of inorganic N and P fertilizers [21], this choice also opens up the possibility of producing highly valuable complementary yield-boosting tools. We thus focused the work presented herein on T. aestivum, the most widely grown wheat species.
Correlations between microbiome composition and soil N/P content have been reported for many plants [22][23][24][25][26][27][28][29][30][31][32]. This phenomenon is largely attributed to concomitant variations of soil microbiome composition, which presumably prompt different unsolicited plant colonization processes [33,34]. But the dynamics of N-and P-sharing plant-microbe partnerships are undoubtedly also at play. The N-related microbiome variations seen in legumes are, for example, largely tied to the plants' modulation of their rhizobia recruitment efforts [35][36][37]. The Prelated microbiome variations seen in many mycorrhizal embroyphytes are similarly tied to the recruitment of arbuscular mycorrhizae [38][39][40]. The prevalence of N/P-sharing plantmicrobe partnerships, and their contribution to N/P-related variations of plant microbiome composition, seem to extend much beyond the above examples [41][42][43][44][45]. To assess whether wheat may be among the plants that have undiscovered associations of these types, we thus first aimed to test if its root and rhizosphere microbiomes vary with soil N and P content.
The microbiomes of plant roots and rhizospheres are generally distinguishable from those of the adjacent bulk soils [46]. This phenomenon is again largely attributed to unsolicited microbial colonization. But plant recruitment of N-and P-providing microorganisms is also undoubtedly at play under conditions of low soil N and P. Legumes' activation of flavonoid release and nodule organogenesis mechanisms, for example, promotes such contrasts [36]. Similar patterns are also attributed to other microbial attraction and accommodation mechanisms [42,44]. We thus, as part of the aforementioned assessment, also aimed to test if wheat shapes its root and rhizosphere microbiomes under conditions of low N and P availability.
The Historical Dryland Plots, a unique long-term experiment that examines wheat production with and without the use of N and P fertilizers in Lethbridge, Canada [47], constituted the ideal setting to run the selected tests. Indeed, its Rotation A section (wheat without fallow since 1911) contains a plot that annually receives N and P fertilization (45 kg N ha -1 ammonium nitrate and 20 kg P ha -1 monocalcium phosphate, i.e. N 45 P 20 ) and plots that have not been amended with one or both for over 4 decades (i.e. N 45 P 0 , N 0 P 20 , N 0 P 0 ). It thus offers an access to T. aestivum that are both grown under sufficient N/P supplies and in conditions of extreme N/P-starvation on the same land. We therefore profiled the bacterial and fungal communities that were associated with T. aestivum roots, rhizospheres, and bulk soil collected from each Rotation A plot, and tested: 1) if the nutrient-starved T. aestivum showed root and rhizosphere microbial communities that were significantly different from those of the fertilized plants (hypothesis 1 , Fig 1), and 2) if starvation-specific T. aestivum communities were also significantly different from those of the adjacent bulk soils (hypothesis 2, Fig 1). As they may represent soil microbes that are specifically attracted towards the plant when it experiences Nor P-starvation, a behavior that is seen in known N-and P-providing plant partners [46, 48,49], we then identified all the OTUs that were proportionally enriched in one or more of the nutrient starvation-and plant-specific communities.

Sampling and DNA extraction
DNA samples analyzed in this study were provided by B. Helgason (Agriculture and Agri-Food Canada and University of Saskatchewan) and S. Hemmingsen (National Research Council Canada). Plant and soil samples (S1 Table) were collected from the Rotation A plot of the Rotation ABC study at the AAFC Lethbridge Research Center in 2013 and 2014, when the plants were in early vegetative growth and anthesis, following the procedure described by Siciliano and Germida [50]. DNA was extracted using DNeasy PowerSoil and PowerPlant Pro DNA Isolation Kits (Qiagen, Hilden, Germany) from bulk soil (250 mg per sample), rhizosphere soil (250 mg per sample), and washed root samples (50 mg per sample) collected by Helgason (Hemmingsen personal communication). Rotation A plots N 45 P 0 , N 0 P 20 , and N 45 P 20 were supplemented once a year with 45 kg N ha-1 ammonium nitrate and/or 20 kg P ha-1 monocalcium phosphate [47]. AC Lillian was planted during the period under study.

PCR amplification/sequencing of microbial taxonomic markers
Genes encoding for a portion of the V4 region of bacterial 16S rRNA were amplified from DNA extracts by PCR, with the following reagents and conditions: 1X HotStarTaq Plus Master Mix (Qiagen, Hilden, Germany), 0.4 mg/mL bovine serum albumin (Roche Diagnostics, Basel, Switzerland), and 0.6 uM each of forward (F520 5'-TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG AGC AGC CGC GGT AAT-3') and reverse (R799 5'-GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA GCA GGG TAT CTA ATC CTG TT-3', based on Huws et al. [51] MiSeq primers (Integrated DNA Technologies, Coralville, USA) per reaction; 5 min at 95˚C, followed by 25 cycles of 30 sec at 95˚C, 30 sec annealing at 45˚C, 45 sec at 72˚C, and a final 10 min at 72˚C. Fungal internal transcribed spacers were also PCR amplified using the same reagents and conditions as bacterial 16S rRNA genes, except for the forward (ITS_F 5'-TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG CTT GGT CAT TTA GAG GAA GTA A -3') and reverse (ITS_R 5'-GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA GCT GCG TTC TTC ATC GAT -3', based on Martin and Rygiewicz [52]) primers.
PCR amplicons were purified with Beckman Coulter's Agencourt AMPure XP system (Beckman Coulter Inc, Brea, USA) according to the manufacturer's instructions. They were then indexed by PCR with the following reagents and conditions: 5 uL of purified amplicon solution, 1X KAPA HiFi HotSart ReadyMix (Kapa Biosystems, Wilmington, USA), and 2.5 uL each of Nextera XT DNA Library Prep (Illumina, San Diego, USA) index 1 and index 2 primer; 3 min at 95˚C, followed by 8 cycles of 30 sec at 95˚C, 30 sec annealing at 55˚C, 30 sec at 72˚C, and a final 5 min at 72˚C. Indexed PCR amplicons were again purified using Beckman Coulter's Agencourt AMPure XP system, quantified with a Quant-iT PicoGreen dsDNA Assay Kit (ThermoFisher Scientific, Waltham, USA), pooled, size selected with a SPRISelect DNA Size Selection system (Beckman Coulter Inc, Brea, USA), and sequenced on a MiSeq DNA Sequencer (Illumina, San Diego, USA). Sequence data was deposited into the NCBI's sequence read archive under bioproject PRJNA343655 and biosamples SAMN05791727 to SAMN05791930.

Sequence quality control, collation of sequences into OTUs, normalization of OTU abundances, and assignment of taxonomic qualifiers
Sets of 16S rRNA gene and ITS sequence reads were all analyzed with a custom amplicon analysis pipeline [53]. Briefly: common sequence contaminants were removed from raw reads using a kmer matching tool (DUK v1.051, http://duk.sourceforge.net/), filtered reads were assembled using FLASH v1.2.11 [54], assembled amplicons were trimmed with custom Perl scripts to remove remaining primer sequences and filtered for quality (sequencing with >1 N, an average quality score lower than 33, or more than 5 nucleotides having a quality score lower than 10 were rejected).
OTU generation was conducted using a three-step clustering pipeline. Quality controlled sequences were dereplicated at 100% identity using a custom Perl script, denoised at 99% identity using USEARCH v.6.0.203 [55]. Clusters of less than three reads were discarded and remaining clusters were scanned for chimeras using UCHIME, first in de novo mode, then in reference mode, using the Broad Institute's 16S rRNA Gold reference database. Remaining clusters were clustered at 97% identity (USEARCH).
16S rRNA gene OTUs were assigned taxonomic qualifiers using the RDP classifier (v2.5) with a modified Greengenes training set built from a concatenation of the Greengenes database (version 13_8 maintained by Second Genome), Silva eukaryotes 18S r118, and a selection of chloroplast and mitochondrial rRNA sequences. For ITS OTUs, this task was performed using the ITS Unite database. Hierarchical tree files were generated with custom Perl scripts and used to generate training sets using the RDP classifier training set generator's functionality [56]. With taxonomic qualifiers in hand, OTU abundances were normalized with edgeR v3.10.2 [57].

Estimated microbial community coverages, comparisons of community structures, identification of OTU enrichments
The proportion of each microbial community that was captured through the sequencing efforts described above was estimated using R package Entropart (Chao estimator) [58].
We then sought to determine whether significant differences of mean community dissimilarity were found among the bacterial root data of the four Rotation A plots (N 0 P 0 , N 0 P 20 , N 45 P 0 , N 45 P 20 ). This was accomplished by running a permanova analysis on a matrix of community dissimilarity (Bray-Curtis index) created with the associated sample profiles (9999 permutations, R package Vegan version 2.4-1) [59]. This procedure was subsequently repeated with bacterial rhizosphere, bacterial bulk soil, fungal root, fungal rhizosphere, and fungal bulk soil data. Datasets that showed significant differences (p < 0.05) were further explored with post hoc permanova analyses that compared plot-specific data subsets against each other (e.g. bacterial root N 45 P 20 vs. N 0 P 20 for N starvation and N 45 P 20 vs. N 45 P 0 for P starvation, hypothesis 1- Fig 1, analyses parameters as above). The p values of these individual pairwise comparisons were adjusted with the Bonferroni correction [60] to account for multiple testing. The tests that revealed significant differences (p < 0.05 after adjustment) were followed by tests for multivariate homogeneity of group dispersion conducted with Vegan function betadisper and R stats package function tukeyHSD (default parameters, p values also adjusted with Bonferroni correction). Non-metric multidimensional scaling (NMDS) ordinations were performed with Vegan's metaMDS function to illustrate community differences (trymax = 100, all other parameters default).
We then sought to determine whether the identified starvation-specific communities were also significantly different from their bulk soil counterparts (e.g. bacterial root N 0 P 0 vs. bacterial bulk soil N 0 P 0 , hypothesis 2- Fig 1). This was also accomplished by running permanova analyses on matrices of community dissimilarity (parameters as above with Bonferroni corrections). The tests that revealed significant differences (p < 0.05 after adjustment) were again followed by tests for multivariate homogeneity of group dispersion (parameters as above), and community differences were again illustrated with NMDS ordinations. These analyses identified 6 N/P-starvation and T. aestivum-specific microbial communities (i.e. that were simultaneously different from the equivalent communities in fertilized T. aestivum and from the communities of their adjacent bulk soil). A conservative two-step process was used to identify the OTUs that were enriched in each one. All OTUs that were detected in the relevant data subsets (e.g. bacterial root N 0 P 0 , bacterial root N 45 P 0 , bacterial bulk soil N 0 P 0 ,) were individually submitted to a two-tailed Kruskal-Wallis Rank Sum Test using function kruskal.test from R stats package (default parameters). Those that had a significantly different average relative abundance among the three data subsets were then assessed for statistical enrichments in the starvation-and T. aestivum-specific data subset (e.g. bacterial root N 0 P 0 ). This was accomplished by fitting a linear model to each OTU's overall abundance in the three considered data subsets (using R stats package lm function), and repeatedly testing it for abundance differences in pairs of compartment-and plot-specific data subsets (e.g. bacterial root N 0 P 0 vs bacterial root N 45 P 0 , bacterial root N 0 P 0 vs bacterial bulk soil N 0 P 0 ) using R package Multcomp's general linear hypotheses testing function (glht) with the treatment difference set to "Tukey" [61]. To illustrate the results of these analyses, we presented the average relative abundances of each enriched OTU next to a dendrogram built with the average Bray-Curtis dissimilarities of the associated compartment-and plot-specific data subsets. Heatmap were drawn using R package pheatmap [62], average community similarities were calculated using Vegan's meandist function, and dendrograms were drawn with function of hclust of R's stats package.

Estimated microbial community coverages, comparisons of community structures, identification of OTU enrichments
The sequencing datasets associated with the 172 samples contained, on average, 50052 bacterial 16S rRNA gene sequences or 77115 fungal ITS sequences following quality control. Evaluations of community coverage calculated using a Chao estimator suggest that this sequencing effort allowed the identification of, on average, 98.1% of the bacterial OTUs and 99.4% of the fungal OTUs that were present in each sample.
Significant differences of mean community dissimilarity were found among the bacterial root data of the four Rotation A plots (p = 0.0001, R 2 = 0.13), and the difference found in a pair of N-starved/N-fertilized data subsets contributed to this signal (N 45 P 0 /N 0 P 0 , p = 0.0004, R 2 = 0.13, Fig 2A). An additional comparison revealed that the N 0 P 0 root data subset was also significantly different from that of N 0 P 0 's bulk soil (p = 0.0002, R 2 = 0.51, Fig 3E, Table 1). The latter subsets did, however, show significantly different group dispersions (p < 0.05).
The bacterial rhizosphere communities of the four Rotation A plots showed significant differences of mean dissimilarity (p = 0.0001, R 2 = 0. 27). The differences found in two pairs of N-starved/N-fertilized data subsets (N 45 P 20 /N 0 P 20 : p = 0.0004 and R 2 = 0.18, N 45 P 0 /N 0 P 0 : p = 0.0004 and R 2 = 0.27, Fig 2C) and one pair of P-starved/P-fertilized data subsets (N 45 P 20 / N 45 P 0 : p = 0.0176 and R 2 = 0.12; Fig 2C) contributed to this signal. Only two of the latter three nutrient-starved subsets were also significantly different from their bulk soil counterparts (N 0 P 20 : p = 0.2649 and R 2 = 0.07, N 45 P 0 : p = 0.0003 and R 2 = 0.15, N 0 P 0 : p = 0.0003 and R 2 = 0.17; Fig 3A, 3C and 3E; Table 1). None of the above comparisons were performed between data subsets that showed significant differences of group dispersion.
The fungal rhizosphere communities of the four Rotation A plots showed significant differences of mean dissimilarity (p = 0.0001, R 2 = 0. 21). The difference found in two pairs of Nstarved/N-fertilized data subsets (N 45 P 20 /N 0 P 20 : p = 0.0012 and R 2 = 0.18, N 45 P 0 /N 0 P 0 : p = 0.0008 and R 2 = 0.17, Fig 2D) contributed to this signal. However, neither of the nutrientstarved subsets were also significantly different from their bulk soil counterparts (N 0 P 20 : p = 0.2322 and R 2 = 0.06, N 0 P 0 : p = 0.1728 and R 2 = 0.06; Fig 3B and 3F; Table 1). None of the above comparisons were performed between data subsets that showed significant differences of group dispersion.
Significant differences of mean dissimilarity were, finally, also found among the bacterial and fungal bulk soil data of the four Rotation A plots (respectively p = 0.0001, R 2 = 0.21 and p = 0.0001, R 2 = 0.16). The differences found in pairs of N-starved/N-fertilized data subsets contributed to this signal in both cases (bacteria N 45 P 20 /N 0 P 20 : p = 0.0016 and R 2 = 0.13, bacteria N 45 P 0 /N 0 P 0 : p = 0.0004 and R 2 = 0.28, Fig 2E; fungi N 45 P 20 /N 0 P 20 : p = 0.0004 and R 2 = 0.12, fungi N 45 P 0 /N 0 P 0 : p = 0.0004 and R 2 = 0.12, Fig 2F). None of the above comparisons were performed between data subsets that showed significant differences of group dispersion.
Our sequencing effort globally uncovered 8692 bacterial OTUs and 3642 fungal OTUs. The 100 that had the highest global average relative abundances in each set are respectively showed in Fig 4 and Fig 5. The 9 bacterial orders that were the most represented among these OTUs were Saprospirales, Actinomycetales, Burkholderiales, Rhizobiales, Rubrobacterales, Xanthomonadales, Sphingomonadales, Cytophagales, and Flavobacterales. The 9 fungal orders that were the most represented among these OTUs were Pleosporales, Sordariales, Hypocreales, Agaricales, Mortierellales, Helotiales, Xylariales, Chaetothyriales, and Eurotiales.
Twenty-one of the 8692 identified bacterial OTUs had higher proportional representations in one of the 3 N/P-starvation and T. aestivum-specific bacterial communities. In 17 cases, the OTUs had higher average relative abundances in the roots of the N 0 P 0 wheat than in the roots of the N 45 P 0 wheat or the N 0 P 0 bulk soil (Fig 6A). These OTUs belonged to genera Lentzea (order Pseudonocardiales), Methylibium (order Burkholderiales), Cellulosimicrobium (order Micrococcales), Chitinophaga (order Chitinophagales), Nocardioides (order Propionibacteriales), Polaromonas (order Burkholderiales), families Actinosynnemataceae (order Pseudonocardiales), Caulobacteraceae (order Caulobacterales), Cytophagaceae (order Cytophagales), Streptomycetaceae (order Streptomycetales), Chitinophagaceae (order Chitinophagales), and order Actinomycetales. In 3 cases, the OTUs had higher average relative abundances in the rhizosphere of the N 0 P 0 wheat than in the rhizosphere of the N 45 P 0 wheat or the N 0 P 0 bulk soil (Fig 6B). These OTUs belonged to genera Erwinia (order Enterobacterales), Variovorax (order Burkholderiales), and family Chitinophagaceae (order Chitinophagales). In one case, the OTU had higher average relative abundances in the rhizosphere of the N 45 P 0 wheat than in the rhizosphere of the N 45 P 20 wheat or the N 45 P 0 bulk soil (Fig 6C). This OTU belonged to genus Agrobacterium (order Rhizobiales).
Eight of the 3642 identified fungal OTUs had higher proportional representations in one or two of the 3 N/P-starvation and T. aestivum-specific fungal communities. In 2 cases, the OTUs had higher average relative abundances in the roots of the N 0 P 20 wheat than in the roots of the N 45 P 20 wheat or the N 0 P 20 bulk soil (Fig 7A). These OTUs belonged to genera Lulworthia (order Lulworthiales) and class Sordariomycetes. In 5 cases, the OTUs had higher average relative abundances in the roots of the N 0 P 0 wheat than in the roots of the N 45 P 0 wheat or the N 0 P 0 bulk soil (Fig 7B). These OTUs belonged to genera Apodus (order Sordariales), Conocybe (order Agaricales), Crocicreas (order Helotiales), and phylum Ascomycota. In 2 cases, the OTUs had higher average relative abundances in the roots of the N 45 P 0 wheat than in the roots of the N 45 P 20 wheat or the N 45 P 0 bulk soil (Fig 7C). These OTUs belonged to genera Parastagonospora (order Pleosporales) and Phaeosphaeriopsis (order Pleosporales).

Discussion
Our analyses identified four microbial communities that were specifically found in the roots or rhizospheres of N-starved T. aestivum. These communities were, of course, most probably shaped largely by factors that lie outside of plant-microbe partnerships. Their distinctiveness from the corresponding communities of N-fertilized plots was for exmple likely linked to the differences that existed between N 45 and N 0 bulk soil communities. Soil microbiome disparities are notoriously seen in agricultural fields that receive different N fertilization [33, 34, 63,64]. They generally spread to plant microbiomes, and this seems to occur mostly through opportunistic microbial colonization [63,64]. The four communities' uniqueness could, however, also be partly tied to the activation of N-sharing plant-microbe partnerships in the N 0 plots. Plant-associated communities that are molded by such interactions generally do, indeed, show N-dependent variations. Ikeda et al. [26,27] reported on this phenomenon in soybeans.
Others have similarly reported on correlations between the plant colonization efficiency of known N-providing plant partners and N soil content [22,24]. This idea is also consistent with  aestivum's N-sharing partnerships could be analogous. The microorganisms that were potentially involved in the N 0 plots were, however, not the most common N-providing plant partners. Indeed, none of the OTUs that were enriched in the four N-starvation and plant-specific communities belong to orders Rhizobiales, Rhodospirillales, Cyanobacteriales, or Glomerales. Several of the enriched OTUs do, however, belong to bacterial taxa that contain known diazotrophic plant partners (i.e. genus Erwinia, family Enterobacteriaceae, order Burkholderiales, class Actinobacteria) [42,65]. Several others belong to taxa that were enriched in the microbiomes of N-deficient plants in other experiments (i.e. order Burkholderiales, families Cytophagaceae, Chitinophagaceae, Caulobacteraceae, and Comamonadaceae) [28,30]. These OTUs may represent soil microbes that are specifically attracted towards T. aestivum when it experiences N-starvation, a behavior that is seen in known N-providing plant partners [46,49].
The above observations do not constitute evidence for the existence of undiscovered Nsharing wheat-microbe partnerships, but they are coherent with it. Other such observations include links between the plant's cultured root microbiome and soil N content [15], links between the plant's metabolomic profile and soil N content [16,20], the fact that it can shape its root microbiome through jasmonic acid production [17], the fact that its rhizosphere contains putative N-fixing microorganisms [18], and the fact that T. aestivum can directly receive N from microbes in experimental conditions [19]. Interestingly, several taxa that were represented on the list of OTUs enriched in the N-starvation and plant-specific communities contain other types of known plant growth promoting microbes. This is the case of bacterial genera Methylibium, Polaromonas, Variovorax, Erwinia, bacterial family Chitinophagaceae, which all contain known ACC deaminase producers [67][68][69][70], and bacterial genera Cellulosimicrobium, bacterial order Actinomycetales, which contain known antibiotic producers [71][72][73]. N starvation may thus also trigger non-N-sharing T. aestivum-microbe partnerships.
Our analyses also identified two microbial communities that were specifically found in the roots or rhizospheres of P-starved T. aestivum. These communities were also most probably shaped largely by factors that lie outside of plant-microbe partnerships. Their associated bulk soil communities were not, however, different from those associated with their P 20 counterparts. In light of research that demonstrated correlations between plant microbiome composition and soil P content [32], correlations between the plant colonization efficiency of known P-providing plant partners and P soil content [23,25,29], and the root microbiome defining effects of microbial partner accommodation mechanisms in A. thaliana [44], these results suggest that the two communities may have been partly shaped by P-sharing plant-microbe partnerships. Wheat is a mycorrhizal plant [12][13][14], so the identification of a P-starvation-and T. aestivum-specific fungal root community was expected. The OTUs that were enriched in that community do not, however, belong to taxa known to contain fungal wheat partners (e.g. genera Glomus, Sclerocystis, Acaulospora, Scutellospora) or any other P-providing fungal plant partner (i.e. Glomeromycetes, Agaricomycetes) [41]. They rather belong to two taxa that contain notorious plant pathogens: genera Phaeosphaeriopsis and Parastagonospora [74,75]. This observation is consistent with a P starvation-specific depression of T. aestivum immunity, a phenomenon that could be similar to that described in A. thaliana by Castrillo et al. [44]. The OTU that was enriched in the P-starvation-and T. aestivum-specific bacterial rhizosphere community does, however, belong to a genus known to contain phosphate-solubilizing plant growth promoters [76,77].

Conclusions
Six N/P starvation-and plant-specific microbial communities were identified. The OTUs that were enriched in these communities may represent microbes that are specifically attracted towards T. aestivum when it experiences N-or P-starvation, a behavior that is seen in known plant partners [46, 48,49]. Many do, indeed, belong to taxa containing relevant N/P-providing plant partners. These results are consistent with the existence of undiscovered N/P-sharing wheat-microbe associations. Additional research will be needed to validate this interpretation. But the work presented here provides a way forward. The identification of potential T. aestivum partners gives us a target list for subsequent relationship-assessing studies, which is in line with modern microbiome research efforts that promote the identification of potentially beneficial microbes and their use in experimental system manipulations [78][79][80]. Wheat farming currently consumes approximately 20% of the worldwide production of inorganic N and P fertilizers, the latter experiments could thus pave the way for the development of valuable complementary yield-boosting tools.
Supporting information S1