The effect of wheat genotype on the microbiome is more evident in roots and varies through time

Crop breeding has traditionally ignored the plant-associated microbial communities. Considering the interactions between plant genotype and associated microbiota is of value since different genotypes of the same crop often harbor distinct microbial communities which can influence the plant phenotype. However, recent studies have reported contrasting results, which led us to hypothesize that the effect of genotype is constrained by growth stages, sampling year and plant compartment. To test this hypothesis, we sampled bulk soil, rhizosphere soil and roots of 10 field-grown wheat genotypes, twice per year, for 4 years. DNA was extracted and regions of the bacterial 16 S rRNA and CPN60 genes and the fungal ITS region were amplified and sequenced. The effect of genotype was highly contingent on the time of sampling and on the plant compartment sampled. Only for a few sampling dates, were the microbial communities significantly different across genotypes. The effect of genotype was most often significant for root microbial communities. The three marker genes used provided a highly coherent picture of the effect of genotype. Taken together, our results confirm that microbial communities in the plant environment strongly vary across compartments, growth stages, and years, and that this can mask the effect of genotype.


INTRODUCTION
Crop breeding has traditionally been carried out under high input conditions and at the same time, ignoring plant-associated microbial communities [1][2][3], leading to the view that modern genotypes have lost their ability to associate with nutritionbeneficial microbial communities when grown under low inputs [4][5][6][7][8][9]. Since many plant-associated microorganisms can positively affect the plant phenotype by increasing nutrition, deterring pathogens, promoting growth and reducing stress, their absence could significantly hamper efforts toward sustainable agriculture. However, plant-associated microorganisms are in turn influenced by seasonal variation in environmental conditions [10], plant developmental stage [11,12], plant compartment [11,13], and their interactions [14], which probably shape the strength of the association between microbial communities and particular plant genotypes.
Microbial communities were shown to vary with plant developmental stage but also with environmental conditions throughout the growing season. In the field, both sources of variation are confounded, making it difficult to tease apart the influence of development stages and plant and microbial responses to variation in environmental conditions. However, studies under controlled environmental conditions showed that plant developmental stage does influence the microbial community diversity [11,15,16] and activity [12], through changes in the composition of the root exudates [17,18]. Similarly, field studies have shown that soil microbial communities are influenced by time through changes in environmental conditions, such as precipitation patterns [19] nutrient availability [20], and temperature [21]. For instance, Wang et al. [19] showed that the effects on microbial communities of dry spells followed by rewetting, dwarfed the effects of reduced soil water content and of different wheat genotypes. Plants also respond to environmental conditions by modulating their rhizodeposition [22,23], providing an indirect way for seasonal variation in environmental conditions to affect microbial communities.
Plant compartment is often reported as the most dominant factor influencing the diversity of plant-associated microbial communities [11,24,25]. One of the most reported patterns is the difference between the rhizosphere and bulk soil microbial communities, dubbed the rhizosphere effect [6,[26][27][28]. Likewise, root, rhizosphere and aboveground plant microbial communities are distinct [11,29,30], these latter sharing similarities with the seed microbial communities [31]. These differences are due to a combined selective pressure of quantity and quality of nutrients [32], microbial capacity to invade plant tissues, and plant immune response to invaders [27,33,34]. For instance, assembly in the rhizosphere is linked to the presence of various chemicals exuded by the roots [18,35,36], to the capacity of the microbes to consume these exudates or react to them and to form biofilms [37,38]. Invasion of the plant tissue, such as the root endosphere, requires the microbes to evade plant defenses and to adapt to life inside tissue where nutrient sources are highly unbalanced.
A recent study from our group showed that spatial and temporal factors interact to modulate microbial communities [11]. This suggested that depending on the plant compartment, the effect of plant developmental stages on the microbial communities is not identical. It would therefore be expected that an effect such as plant genotype, that is generally reported to have a more subtle effect on microbial communities than environment or plant compartments [19,39], would be influenced and even constrained by these factors. One of our recent studies showed that the rhizosphere metagenomes of 10 different field-grown wheat genotypes were nearly identical [9]. However, previous work from our group did highlight significant differences between wheat genotypes in term of function, community structure and composition for wheat growing in pots in a growth chamber [39,40], in commercial fields [41], or in an experimental field [19]. We therefore hypothesized that the effect of wheat genotype might show some spatio-temporal variability, being only visible at certain growth stages, in certain plant compartments or under certain environmental conditions, which could explain the lack of significance observed in Quiza et al. [9]. Here, we expanded the analysis presented in Quiza et al. [9] by sampling the same field experiment, but twice per season, over four growing seasons, and by including root and bulk soil samples in addition to the rhizosphere soil. We used amplicon sequencing targeting the bacteria and archaea (16 S rRNA and CPN60 genes) and the fungi (ITS region) to characterize the soil microbial community. Our study sheds light on the interactive effects between plant genotype, sampling time, growth stage and compartment and suggests that the genotype effect is plant compartment and sampling time dependent.

MATERIALS AND METHODS Experimental design
A field experiment was conducted from 2013 to 2016 at the Nassar Crop Research Farm of the University of Saskatchewan, Saskatoon, Canada, that has been managed for more than 50 years to conduct experiments under low fertilization conditions. We selected 10 wheat genotypes that were introduced over the period 1845 through 2009. These included six Triticum aestivum or bread wheat genotypes of the Canada Western Red Spring , and CDC Verona (2008)) (https://grainscanada.gc.ca/ en/grain-quality/grain-grading/wheat-classes.html). The experiment was arranged in three randomized blocks, each consisting of ten 6.2 m 2 plots to which the genotypes were randomly assigned each year. Each plot contained eight rows spaced at intervals of 20 cm and were seeded at 320 seeds m -2 on May 25, 2013, June 2, 2014, May 20, 2015, and June 7, 2016. To minimize the effect of the seed source on plant performance, all genotypes were grown from seeds harvested from a common field under low fertilization. To minimize any erroneous effect on productivity measurements due to poor seedling establishment, 15 kg ha -1 of 11-55-0 (NPK) fertilizer was added at seeding. No other fertilized was added after that. Weather conditions (average daily temperature and total monthly precipitation) for the Saskatoon RCS station (World Meteorological Organization station ID: 71496, 52°10'25.000" N, 106°43'08.001" W) for the duration of the experiment (summer 2013-2016) are presented in Table 1. The 1981-2010 averages were collected from the nearby station at Saskatoon Diefenbaker International Airport (World Meteorological Organization station ID: 71866, 52°10'00.000" N, 106°43'00.000" W), as they were not available for the Saskatoon RCS station.

Sampling
The experiment was originally designed to compare the performance of wheat cultivar under low nitrogen. Sampling was therefore done at the stem elongation (June or July) and dough development (August or September) growth stages, which are crucial for wheat N nutrition. The plants were pooled before processing on a clean piece of bench cover. The shoots were cut with sterilized scissors 1 cm above the soil line. The roots were gently separated from the bulk soil by removing all the lumps and leaving the roots with very little soil adhering to them. The roots were cut off from the remaining stem and crown of the plant with sterilized scissors; they were immediately transferred, with the adhering soil, to a 500 ml Erlenmeyer flask containing 200 mL of sterile PBS buffer. All the roots from one plot were placed in one flask. The flask was placed on a shaker at 150 rpm, at 22°C for 25 min. The roots were removed from the PBS solution and rinsed with distilled water until completely clear, finishing the cleaning with a final rinse of sterile water. The PBS containing the remaining soil from the roots was then centrifuged at 2000 ×g for 5 minutes, after which the supernatant was discarded, and the rhizosphere soil collected. The remaining bulk soil was fragmented into pea-size pieces, homogenized and subsampled in 50 ml conical tubes. The roots, rhizosphere soil and bulk soil were immediately stored at -80°C until extraction.   [42] whereas for the fungal ITS1 region the primers ITS1F and 58A2R [43] were used. We also used the gene coding for the group I chaperonins (CPN60, also known as GroEL or hsp60) to have an independent confirmation of the 16 S rRNA gene results. The cpn60 UT was amplified using the type I chaperonin universal primer cocktail containing a 1:3 ratio of H279/H280:H1612/H1613 as previously described [44,45]. The three pools were then loaded on an Illumina MiSeq sequencer and sequenced in-house using a 600-cycles MiSeq Reagent Kit v3.

Bioinformatics and statistical analyses
We retrieved between on average 41,516, 42,757 and 25,134 quality filtered paired reads for the 16 S rRNA gene, the ITS region and the cpn60 gene, respectively (Supplementary Table S1). Sequencing data were analysed using Amplicon Tagger [46] as described previously [47]. All statistical analyses and figure generation were performed in R v.4.0.2. The Bray-Curtis dissimilarity between samples calculated from the ASV relative abundance was visualized by principal coordinate analysis ("cmdscale" function) based on Bray-Curtis dissimilarity matrices ("vegdist" function of the "vegan" package) [48]. The effects of genotype, year, growth stage and compartment on the community composition was tested by permutational multivariate analysis of variance (PERMANOVA) using the "adonis2" function of the "vegan" package. For all the significant (P < 0.05) and marginally significant (P < 0.10) year/compartments/growth stage combinations, pairwise comparisons of the genotypes was performed (function "pairwiseAdonis") [49]. ANOVA analyses were performed with the function "aov" followed by post hoc Tukey's honestly significant difference (HSD) tests ("agricolae" package) to detect differences in relative abundance of phyla and classes between genotypes through the years.

Microbial community structure
The PCoA ordination based on bacterial 16 S rRNA gene amplicon ASVs revealed the overriding effect of compartment on the bacterial community structure, with bulk and rhizosphere samples clustering away from roots samples (Fig. 1). Within each  (Fig. 1). The Permanova confirmed this visual interpretation, with Compartment having the strongest effect on the community structure (highest pseudo-F ratio) followed by Year and Growth stage, and finally by Genotype ( Table 2). All these main effect terms were having a significant effect on the community structure, together with many of the interaction terms, some of them including Genotype (Genotype:Year and Genotype:Compartment) ( Table 2). The ordinations based on the fungal ITS region amplicon ASVs showed a slightly different picture, with an overwhelming effect of sampling year (2013-2014 vs 2015-2016) and growth stage (Fig. 1). Within the clusters due to sampling year, a clear effect of compartment could be observed, with root samples separated from rhizosphere and bulk soil samples (Fig. 1). Here again, this visual interpretation was confirmed by the Permanova analysis, where Year had the strongest effect followed by Growth stage, Compartment and then Genotype ( Table 2). These main effect terms were all highly significant, except for Genotype that was marginally significant (Table 2). Some interaction terms were also significant, but none of them included Genotype.
The ordination based on the cpn60 gene amplicon ASVs showed a picture in between the ones observed for the 16 S rRNA gene and ITS region amplicons (Fig. 1). Compartment and Year clearly had a strong effect on the microbial community structure, but these two terms seem to interact, with a more distinct root community structure for the 2013 samples (Fig. 1). The 2015-2016 showed here again a tight clustering, but where also joined by the 2014 samples (Fig. 1). In Permanova analyses, Year had the strongest effect, followed by Compartment, Growth stage and Genotype, which were all highly significant (Table 2). Several interaction terms were also significant, notably Genotype:Year and Genotype:Growth stage ( Table 2).

Effect of genotype
Because of the significant Year x Compartment x Growth stage interaction for all three amplicon datasets groups and to explore more deeply our initial hypothesis about genotypes, Permanova analyses for the effect of genotype were performed separately for each Year x Compartment x Growth stage combination (Table 3). Since this reduced the number of samples in each analysis and consequently, the statistical power of our approach, we are also reporting and discussing test results that were 0.05 < P < 0.10. For bacterial communities (16 S rRNA gene), the effect of genotype in the roots was significant at all growth stages for year 2013 and 2014, but only at dough development for 2015 and at stem elongation for 2016 ( Table 3). The only occurrence where genotype was significant in the rhizosphere for the 16 S rRNA gene dataset was at stem elongation in 2015 (Table 3). In contrast, the effect of genotype on fungal communities (ITS region) was only significant for the 2013 samples, in the roots for all growth stages and in the rhizosphere for the dough development growth stage ( Table 3). The patterns of significance for the cpn60 gene dataset were like the ones observed for the 16 S rRNA gene dataset, with significant effects of genotype on root communities in 2013 (both growth stages), 2014 (dough development), 2015 (stem elongation) and 2016 (dough development) (Table 3). Additionally, the genotype significantly affected the rhizosphere microbial communities at the stem elongation stage in 2016 (Table 3).
We further explored the effect of genotype on the phylum/class level community composition in the roots for the growth stages where genotype was significant in Permanova. For the 16 S rRNA gene, ANOVA revealed that, in the roots, Actinobacteria  (Table S2 and Fig. 2). In Tukey HSD post hoc tests, we observed that many of the significant differences observed were between durum and aestivum wheat genotypes (Table S2). For fungi, wheat genotypes significantly affected the relative abundance of Ascomycota and Basidiomycota in the roots at stem elongation in 2013 (Table S3 and Fig. 3), and the only significant difference in Tukey HSD post hoc tests was between the relative abundance of Ascomycota between Marquis and CDC Verona (Table S3). Finally, for the cpn60 gene   (Table S4 and Fig. 4). Tukey HSD post hoc tests showed that some of the significant differences were between the older genotypes (Pelissier, Red Fife and Marquis) and some other more recent genotypes (Table S4).

DISCUSSION
Our multi-year field study of the microbiome of different wheat genotypes revealed an overwhelming influence of sampling year and growth stage on the bacterial and fungal communities associated with wheat root and rhizosphere. In addition, the microbiome of the roots was generally well differentiated from the one inhabiting the rhizosphere and bulk soil. Within these strong effects, we could still detect a significant effect of wheat genotype, which, in agreement with our hypothesis, varied with year and growth stages, and was mostly significant for the root bacterial microbiome. These trends were coherent for the two different bacterial amplicons used (16 S rRNA and cpn60 genes). We had previously constructed metagenomic libraries from the rhizosphere samples collected from this experiment at the stem elongation stage in 2013. The lack of effect from the wheat genotype on the metagenome presented in Quiza et al. [9] is highly coherent with the results presented here, as no significant effect of genotype could be observed for the 16 S rRNA gene, the ITS region nor the cpn60 gene in the rhizosphere samples collected at the stem elongation stage in 2013. However, when expanding our sampling, we realized that, for genotype, significance was mostly seen for roots samples, and even then, it varied through the years and wheat growth stages. These results reconcile Quiza et al. [9] with the literature that reported significant effects of genotype for wheat microbial community structure, composition and functions [19,[39][40][41], by suggesting that the samples taken for the Quiza et al [9] study were not affected, but that at other times or in different compartments, this effect would have been apparent. Alternatively, it could very well be that, because of the high functional redundancy of microbial communities, the changes observed here at the taxonomic level have little effect on the functional makeup of the community and would have been undetectable using shotgun metagenomics.
Here, we report for the first time, in a multi-year field study, that the effect of genotype on the wheat microbiome in the field is Table 3. Permanova analyses for the effect of genotype for each year/compartment/growth stage combinations for 16 S rRNA gene, ITS region and cpn60 gene amplicon datasets generated from bulk soil, wheat rhizosphere soil and wheat roots collected from 2013 to 2016 at the stem elongation and dough development stages.  highly variable across compartments, growth stage, and sampling year, with more significant effects in the root compartment. This variability in the genotype-specific effects on root-associated microbial communities challenges the view that modern genotypes would have lost their ability to associate with nutritionbeneficial microbial communities when grown under low inputs [1,2,4,9]. Indeed, for root-and rhizosphere-associated bacteria and fungi of wheat grown under low nutrient conditions, most of the sample categories did not show an influence of genotype, suggesting that if there are any differences between modern and ancient genotypes, it is highly transient. Although changes in function could have potentially happened without concomitant changes in marker gene data, our data is highly coherent with the shotgun metagenomic study previously carried out on a subset of our samples [9]. Alternatively, even though we used two different species of wheat (durum and bread wheat), the relatively short genotype development gradient (~165 years) might have precluded the observation of more significant trends. However, previous studies that found a clear genotype effect were mostly done using controlled growth conditions, or on a single sampling date or single compartment, and it is thus difficult to conclude if our observations would also apply to other crops or along a longer genotype development gradient (e.g. by including wild parental plants).
In sharp contrast with many previous reports, where fungi are often reported to be more intimately linked to plants, and therefore more influenced by variations between genotypes [50] whereas bacteria are more often primarily influenced by soil properties [39,51,52], we found here that bacteria were more strongly influenced by wheat genotype than fungi. In fact, in general Permanova tests, there was no main or interactive effect of genotype on fungal communities, and in more in-depth ANOVA analyses, we only found an effect of genotype on the fungal communities of the root and rhizosphere of wheat in 2013. Another evidence of the less intimate linkage between fungi and the plant in our study is the fact that fungal communities were overwhelmingly influenced by sampling year, and very little by compartment, whereas bacteria were more strongly influenced by plant compartment than by sampling year. Microbe-microbe interactions could partly explain these patterns, with microbes strongly influenced by a factor limiting the response of other microbes. It was also shown that the influence of genotypes on fungal communities increases with stress [39,50], suggesting that the low nutrient conditions under which wheat was planted in the current field experiment were not stressful. However, Quiza et al. [9] reported a decrease in yields as compared to expected yields for the different genotypes under high input conditions, which would indicate some nutrient limitation or an acclimation period.
Under controlled experimental conditions, plant exudation patterns were shown to vary with development stages [17], which results in shifts in microbial communities and activities [11,12]. Plant development stage could also have played a role in the seasonal patterns observed here, although it is difficult to disentangle from the influence of fluctuating environmental conditions under field conditions. Indeed, the changes observed with wheat growth stage could be due to the direct and indirect   (Table 1). Conversely, the rest of the summer was wetter in 2015 and 2016 (179.5 mm and 152.9 mm of rain, respectively) as compared to 2013 and 2014 (65.9 mm and 73.7 mm of rain, respectively) and to the 30-year average (137 mm) ( Table 1). Precipitation is often cited as an important factor shaping soil microbial communities, in view of its influence on soil water content, gas diffusion and redox conditions and the soil processes that they influence. Recently, we showed that microbial communities, and most especially archaeal ammoniaoxidizers, shifted dramatically following a large drying-rewetting event, whereas little change was observed in response to small changes in average soil water content [19]. Soil water stress history is also a major influencing factor for wheat microbial communities [40].
In our study, the effect of genotype did not extend past the rhizosphere, as in no case was there a genotype effect recorded in bulk soil samples. Recent studies had highlighted that an effect of genotype could be seen in the bulk soil [19,39,41], which led to the suggestion that because of the effect of gaseous compounds, the rhizosphere could in practice extend past the few mm of soil surrounding the roots [53]. Alternatively, fungal hyphae extending from the rhizosphere could also have played a role in these previous studies. However, here, as expected, the effect of genotype was more often significant for the root communities, as 13 out of the 16 significant effects of genotype were observed for root communities. Previous studies showed that genotype effect of wheat increased in the following order rhizosphere<roots<leaves [39], which makes sense as the strength of the plant selective pressure is likely to increase inside plant tissues.
Overall, although we did find some significant influence of genotype on root and rhizosphere microbial communities, this effect varied with sampling years and plant growth stages and was mostly visible for bacteria in root samples. This suggests that the development of modern wheat cultivars did not result in stable shifts in the root-associated microbial communities, and that these changes are anyway dwarfed by the shifts caused by changing environmental conditions and plant development stages that occur throughout and across the years. This knowledge is highly relevant for microbiome engineering approaches [13,54], as it highlights the overwhelming strength of environmental selection vs. plant selection.

DATA AVAILABILITY
Raw sequencing data was deposited in the NCBI database under BioProject accession PRJNA947932 (https://www.ncbi.nlm.nih.gov/bioproject/947932). Fig. 4 Community composition of the most abundant phylum/class based on the cpn60 gene amplicons for the growth stage/ compartment combinations that showed a significant (P < 0.10) effect of genotypes in Table 3.