Wheat dwarfing influences selection of the rhizosphere microbiome

The development of dwarf wheat cultivars combined with high levels of agrochemical inputs during the green revolution resulted in high yielding cropping systems. However, changes in wheat cultivars were made without considering impacts on plant and soil microbe interactions. We studied the effect of these changes on root traits and on the assembly of rhizosphere bacterial communities by comparing eight wheat cultivars ranging from tall to semi-dwarf plants grown under field conditions. Wheat breeding influenced root diameter and specific root length (SRL). Rhizosphere bacterial communities from tall cultivars were distinct from those associated with semi-dwarf cultivars, with higher differential abundance of Actinobacteria, Bacteroidetes and Proteobacteria in tall cultivars, compared with a higher differential abundance of Verrucomicrobia, Planctomycetes and Acidobacteria in semi-dwarf cultivars. Predicted microbial functions were also impacted and network analysis revealed a greater level of connectedness between microbial communities in the tall cultivars relative to semi-dwarf cultivars. Taken together, results suggest that the development of semi-dwarf plants might have affected the ability of plants to recruit and sustain a complex bacterial community network in the rhizosphere.

Wheat domestication originated in the Near East and the crop has undergone a series of crosses and modifications leading to the hexaploid bread species, Triticum aestivum L 1 . Wheat is a key crop for global food security, providing 20% of dietary requirements for calories and protein 2 . To achieve food demand of an increasing global population, projections forecast a need to increase wheat production by 11% by 2026 with just 1.8% increase in cultivation areas 3 . During the Green Revolution, Reduced height (Rht) dwarfing genes were introduced in modern wheat cultivars 4 , resulting in high-yielding wheat plants, which, when combined with agrochemical management and optimal conditions, increased yields, without productivity losses caused by lodging 5 . Plant height is a complex trait controlled by several Quantitative Trait Loci (QTLs) and more than 20 Rht genes have been described, however, most of them have low potential for successful breeding programs 6 . Most breeding programs take into consideration the improvement of above ground parts of plants, with little attention to below ground parts due to the inherent challenges of root analysis coupled with the influence of soil type on root traits 7,8 . As such, the effect of domestication and modern breeding on belowground traits in wheat largely remains unclear. However, Bai et al. 9 assessed a diverse set of wheat germplasm [199 double-haploid progeny derived from a cross between Avalon and Cadenza (Triticum aestivum L.)] and found that plant height was positively correlated to some evaluated root traits such as: total root length, seminal laterals length, seminal axes length, seminal laterals surface area, seminal laterals volume and root dry weight. In addition, Figueroa-Bustos et al. 10 reported a positive correlation between plant height and specific root length (SRL). The selection of desirable root traits through "root breeding", would enable plants to explore soil more efficiently, thus improving water and nutrient acquisition 11 . It has also become clear that future breeding programs should consider plant-microbe interactions 12 as plants are complex holobiont organisms influenced both positively and negatively by their microbial communities 13 Table 1. Cultivars chosen for the current study and some characteristics, such as year of release, pedigree and height. a Height was based on field data collected by Shewry et al. 29 and was categorised as described by Pask et al. 30 .
Metataxonomic sequence analysis pipeline. 16S rRNA gene sequences were analysed using the pipe-  22 . Differences in bacterial community structure were investigated by Permutational Analysis of Variance (PERMANOVA) 37 in Paleontological Statistics Software Package for Education and Data Analysis (PAST) 38 . PCoA plots were obtained using the same software. The number of observed OTUs and diversity based on Shannon index were calculated in QIIME. Statistical analyses of alpha diversity indexes were performed in R as described by Pérez-Jaramillo et al. 22 .
Briefly, normality and homogeneity of variances were checked using Shapiro-Wilk test and Levene's test, respectively, and One-way ANOVA and post-hoc test (Tukey HSD) were used to assess differences.

Results and Discussion
Effect of wheat breeding on root morphological traits. Breeding process appears to have altered some root traits, especially root diameter and specific root length (SRL). Tall plants (old cultivars), when compared to semi-dwarf plants (recent cultivars), showed statistically higher mean SRL (tall = 113.63 m.g −1 and semidwarf = 91.62 m.g −1 ) (One-Way ANOVA, p < 0.05). The opposite general trend is observed for root diameter (tall = 0.44 mm and semi-dwarf = 0.49 mm) (One-Way ANOVA, p < 0.05) ( Supplementary Fig. S2). Tall cultivars tend to have longer thinner roots, whereas shorter thicker roots are observed for semi-dwarf cultivars. Pérez-Jaramillo et al. 22 also observed higher SRL values for wild common bean accessions, which could be useful for enhanced water and nutrient uptake in stressed soils, than the lower SRL of short varieties, suggesting semi-dwarf cultivars could be more dependent on fertile soils. Indeed, modern Italian bread wheat cultivars increased their N demand over time when compared to old tall cultivars 46 . Plant height reduction was achieved during the 20th century wheat breeding 47 , though little attention was given to assessing the effects of these programmes on wheat root traits. In an attempt to study whether Rht genes control both shoot height and seedling root growth, Bai et al. 9 evaluated QTLs for plant height, root and seed traits from different wheat lines and identified some coincident QTLs for roots and height, concluding that the introduction of some known Rht dwarfing genes reduced both plant height and root proliferation. Breeding for higher yields in Australian wheat cultivars led to a reduced root length density (RLD) and total root length and increased efficiency of N uptake 48 . For other crops, such as barley, common bean, rice and maize, it has also been shown that both domestication and breeding caused changes in root architectural traits, leading to a differential spatial arrangement of roots 49 .
The rhizospheres of semi-dwarf cultivars are enriched in members belonging to the phyla Acidobacteria and Verrucomicrobia, which are typically more associated with bulk soil rather than the wheat rhizosphere 16 . A closer analysis at genus level, indicates that tall cultivars are enriched with OTUs assigned to genera normally associated with plant growth promotion (PGP), especially members of the phylum Proteobacteria (Table 2)  www.nature.com/scientificreports www.nature.com/scientificreports/ such as Brevundimonas (indole-3-acetic acid (IAA) production and P solubilisation) 50 ; Devosia (N 2 fixation) 51 ; Rhizomicrobium (reduction of nitrate to nitrite) 52 ; Sphingomonas (IAA and gibberellin production) 53 ; Sphingopyxis (IAA production and N 2 fixation) 54 ; Massilia (IAA, siderophore and proteolytic enzyme production) 55 ; Nitrosospira (ammonia-oxidizing bacterium (AOB)) 56 and Bradyrhizobium, Rhizobium, Methylobacterium, Variovorax, Klebsiella and Pseudomonas which typically possess genes contributing to plant-beneficial functions 57 . Although some differentially recruited OTUs within tall rhizospheres have been assigned to genera commonly related to plant growth promotion, these abilities are normally strain-specific, thus caution must be taken when analysing the results. Further work should evaluate these differences based on metagenomic datasets as well as functional characterisation of microbial isolates.
The major metabolic consequence of plant dwarfing is a reduction in response to the plant phytohormone gibberellin 58 and an increase of active endogenous gibberellin levels when compared to wild-type plants not carrying Rht dwarfing genes 59,60 . The correlation that the rhizosphere of short plants with decreased gibberellin sensitivity is associated with an increased colonisation of bacteria normally associated with bulk soil and a decrease of microbial phyla associated with plant growth promotion, suggests that gibberellin sensitivity influences the composition of the root microbiome. It is known that some microbes, such as PGP members of Rhizobiaceae family 61 as well as some Sphingomonas spp 53 . produce this hormone, and we find these bacteria to be differentially more abundant in tall cultivar rhizospheres. It could follow that plant gibberellin insensitivity leads to a reduction in plant-microbe communication in the root environment and could reduce selection of gibberellin-producing microbes in the root zone. Apart from controlling the growth and development of plants, gibberellins also play a role in plant signalling 60 , affecting defence response mechanisms dependent on jasmonic acid (JA) or salicylic acid (SA) 62 , rhizobial infection of legumes 63 mediated via DELLA proteins and can suppress arbuscular colonisation and development of mycorrhizal symbiosis 64 . Work should be conducted to determine if the root microbiome of these crop accessions diverges further when tall and shorts cultivars are cultured under low nutrient conditions. This would support the hypothesis that the tall cultivars are more capable of recruiting beneficial microbes, and when under stress they do this more readily than short cultivars. Furthermore, microbes have been shown to be capable of producing a range of phytohormones 65 . It will be interesting to determine if the plant breeding process has influenced the plant hormone status for other classes of phytohormones including abscisic acid, auxins, brassinosteroids, ethylene, jasmonates, salicyclic acid and strigolactones, as there is a strong interaction among them 66 and whether their levels influence microbiome composition.
In addition to plant hormonal levels, the observed differences in microbiome structure between short and tall cultivars could be attributed to root exudation profiles 23 . Iannucci et al. 67 showed that even though soil type dramatically changed the composition of root exudates, domestication and breeding also had major effects on root exudates in the rhizosphere of ten tetraploid wheat genotypes. The sensing and metabolism of root metabolites are likely important factors in determining microbiome assembly 68 . Effect of wheat breeding on key rhizosphere taxa and network complexity. We found that for both tall and semi-dwarf networks, zero network hubs were identified and that most nodes were categorised as peripherals (Fig. 3A,B). However, fifteen nodes were classified as module hubs in tall cultivars and they represented a variety of taxonomic groups: five belonged to Proteobacteria (genera Pseudomonas, Luteibacter, Sorangium and orders Myxococcales and Xanthomonadales), three to Verrucomicrobia, two to Acidobacteria and one each from Bacteroidetes, Chloroflexi, Cyanobacteria, Gemmatimonadetes and Latescibacteria. Ten module hubs were identified in semi-dwarf cultivars and three belonged to Proteobacteria (one from ß-Proteobacteria, one from the family Rhodospirillaceae and one from the family Nitrosomonadaceae), three from Bacteroidetes (genera Parafilimonas and Fluviicola and family Saprospiraceae), three from Verrucomicrobia and one from Latescibacteria. Connectors were also detected in both cases, however, semi-dwarf cultivars displayed only  www.nature.com/scientificreports www.nature.com/scientificreports/ four connectors, identified as belonging to genera Asticcacaulis and Pseudohongiella, both from phylum Proteobacteria, genus Luedemannella (Actinobacteria) and one Acidobacteria genus (RB41). On the other hand, tall cultivars showed fifteen connectors: two Acidobacteria (classes Subgroup 22 and Holophagae), one Actinobacteria (class MB-A2-108), one Armatimonadetes (family Fimbriimonadaceae), one Bacteroidetes (genus Pedobacter), one Chlorobi (order Ignavibacteriales), one Chloroflexi (class S085), one Planctomycetes (family Tepidisphaeraceae), two Verrucomicrobia (one identified as the genus Chthoniobacter and other belonging to the family Verrucomicrobiaceae) and five Proteobacteria (two belonging to the genus Haliangium, one betaproteobacterium from the family Rhodocyclaceae, one deltaproteobacterium from the family Sandaracinaceae and one α-Proteobacterium from the order Rhodospirillales (family DA111)). Members of the phylum Proteobacteria were found to be the most prominent putative keystone taxa, accounting for more than 30% of all module hubs and connectors in both tall and semi-dwarf cultivar networks. Members of Bacteroidetes and Verrucomicrobia were the second most prominent keystone taxa for semi-dwarf cultivars, each accounting for 21.43% of module hubs. On the other hand, members of Bacteroidetes accounted for less than 7% of module hubs and connectors and Verrucomicrobia members accounted for 16.67% in tall cultivars. Module hub and connector OTUs have been proposed as putative keystone taxa and ecological generalists, critical for maintaining community structure and function; peripherals on the other hand being considered isolated specialists 41,[69][70][71][72][73][74] . As such, Proteobacteria seem to be critical in maintaining the structure and functionality of bacterial communities in wheat rhizosphere for both tall and semi dwarf cultivars and these results are in accordance with data provided by Oberholster et al. 71 , who also found Proteobacteria to be important keystone taxa for sorghum and sunflower rhizospheres.
Phylogenetic molecular ecological networks (pMENs) were obtained for both tall and semi-dwarf cultivars and the overall structure of each pMEN is presented in Supplementary Table S1. Topological properties indicated that R 2 values for both networks (0.907 for tall cultivars and 0.899 for semi-dwarf cultivars) followed a power-law model, indicating that only a few nodes have many connections with other nodes 41 . Average path distance (GD) obtained for tall cultivars was lower than semi-dwarf cultivars, indicating that nodes in the network from tall cultivars are closer to each other. The average degree (avgK) obtained for tall cultivars indicated a more complex network structure, as a higher value was obtained. This higher complexity was also highlighted by networks calculated with sparCC represented in Fig. 3C,D. High values of modularity for both networks (tall = 0.749 and semi-dwarf = 0.890) indicate that they can be separated into modules using fast greedy modularity optimisation. Modularity values range from 0 to 1 and they measure the degree to which the network was organised into clearly delimited modules 74 and modular structure of complex networks display critical roles in their functionality 75 . The number of modules obtained for tall cultivars was ninety-four, whereas semi-dwarf cultivars resulted in one-hundred and fifteen modules. A module is defined as a group of OTUs which are highly connected, however these are poorly connected with OTUs outside a given module 41 . We found that modules comprised of 4 nodes or more numbered 38 for dwarf varieties compared to 21 for tall varieties, indicating that the OTUs are better connected within the tall cultivar network (Supplementary Fig. S5).
The number of edges obtained in the pMEN from bacterial communities from tall cultivars was higher than those from semi-dwarf cultivars and this was also highlighted by the connectedness index which was 0.498 as opposed to 0.293 for semi-dwarf cultivars. This index ranges from 0 to 1, with values close to 1 indicating a highly-connected graph 41 . The ratio of positive to negative co-occurrence patterns for tall cultivars was 3.094, whereas the ratio of positive to negative interactions for semi-dwarf cultivars was 5.331, indicating that bacterial communities from semi-dwarf cultivars displayed more positive co-occurrence relationships than tall cultivars. The composition and number of nodes per modules varied and the composition of modules is shown in Supplementary Fig. S6. Similarly to nodes identified as keystone taxa, the majority of nodes were classified as Proteobacteria, followed by Acidobacteria, Bacteroidetes, Actinobacteria and Verrucomicrobia. For tall cultivars, the biggest module (#6) comprised of 70 nodes, with more than 32% of nodes belonging to Proteobacteria. As for semi-dwarf cultivars, the largest module (#1) was composed of 71 nodes, with more than 35% belonging to Proteobacteria.
In the present study, tall cultivars showed bacterial lower diversity than semi-dwarf cultivars as well as a more complex bacterial network structure in the rhizosphere. Shi et al. 69 observed that microbial diversity decreased as network size and connectivity increased, which resulted in an increased community organisation. It has also been shown that more complex networks are able to cope with environmental changes 76 , increase crop productivity 73 or suppress soil borne pathogen infection on plants 77 . Indeed, semi-dwarf wheat cultivars carrying Rht alleles were previously shown to be more susceptible than wild cultivars to initial infection by some pathogens 78,79 . This might be linked to the severity of DELLA effect on plant stature, however, differences in differential abundance of specific bacterial keystone taxa might point out some microbial influence on disease susceptibility. For instance, in the network from tall cultivars, OTUs 1135 and 427 in modules 2 and 3, respectively, were classified as belonging to the genus Haliangium, whose species have been described as putative biocontrol agents 80 and another three-keystone species, from the same order, Myxococcales, were also observed in tall cultivars. Myxobacteria include Gram-negative gliding bacteria with predatory features, mostly due to the direct cell to cell contact, production of hydrolytic enzymes and secondary metabolites with antibiotic activity 81 . The third module in tall cultivars show one module hub (OTU 575 -classified as belonging to the genus Chryseobacterium) which has previously been described in the rhizosphere of wheat [14, 82;], and these bacteria have been shown to display plant growth-promotion 82 and are enriched in the presence of 2,4-DAPG-producing species, such as Pseudomonas fluorescens 83 . Interestingly, the other module hub identified in the same module as Chryseobacterium was OTU 171 (assigned to the genus Pseudomonas) and they positively co-occur, meaning that their abundance changed along the same trend but not necessarily meaning they directly interact with each other 74 . OTU 280, assigned to the genus Pedobacter and OTU 563 (Luteibacter) have also been identified as keystone taxa in tall cultivars. They have been found to inhibit the growth of root pathogens, such as the fungus Rhizoctonia solani 84,85 . When comparing both networks, only two OTUs (601 and 118) assigned to the genus RB41 (Blastocatellaceae family from the phylum Acidobacteria) were commonly found and were identified as keystone species. The breeding process seems to have reduced the amount of keystone taxa known for their potential to antagonise pathogens, possibly affecting the degree of resistance of semi-dwarf cultivars to specific diseases. Future work should investigate whether tall and semi-dwarf cultivars select different beneficial bacteria with PGP abilities, as well as antagonistic activities against economically important soil-borne pathogens.

Effect of wheat breeding on 16S rRNA gene-predicted functions. NSTI values obtained for 16S
rRNA gene-predicted functions on rhizosphere soil samples are on average 0.1679 ± 0.013 for all samples and this value is in accordance to what is observed for soils 44 and is lower when compared to other rhizosphere soils 86 . Rhizosphere samples collected from tall cultivars, on average have a NSTI value of 0.1534 ± 0.007 which is lower than NSTI value obtained from the rhizosphere of semi-dwarf cultivars (0.1766 ± 0.006). Functional prediction resulted in two hundred and twenty-one KEGG orthologs (KOs), of which one hundred and forty-one were differentially abundant (p < 0.05). Of these, 45.39% were significantly more abundant in tall cultivars and the remaining 54.61% were significantly more abundant in semi-dwarf cultivars. Seven 16S rRNA gene-predicted pathways were significantly enriched in rhizosphere bacterial communities of tall cultivars, as opposed to nine pathways enriched in the rhizosphere of semi-dwarf cultivars (Fig. 4A). Tall cultivars were predicted to be enriched in functions related to membrane transport, such as ABC transporters (Fig. 4B). Mahoney et al. 14 suggested that differences in 16S rRNA gene-predicted functions observed between wheat cultivars could be driven by differential root exudate chemistry. The breeding process may have affected the quality and quantity of root exudates, impacting predicted functions in the rhizosphere. The presence of bespoke solute binding protein-dependent ABC transporters may be required for the use of root exudates by bacteria, aiding their colonisation of this niche 87 . When analysing specific KOs, tryptophan metabolism was also predicted to be enriched in bacterial communities from the rhizosphere of tall cultivars (Fig. 4B) and this suggests higher production of Scientific RepoRtS | (2020) 10:1452 | https://doi.org/10.1038/s41598-020-58402-y www.nature.com/scientificreports www.nature.com/scientificreports/ indole-3-acetic acid (IAA) hormone, as tryptophan is the main precursor for IAA biosynthesis, acting as a signalling molecule between plants and microbes as well as for plant growth promotion 88 . It could also mean that tall wheat plants are secreting more tryptophan upon colonization by specific bacterial groups 89 . On the other hand, bacterial communities from semi-dwarf cultivars were predicted to be enriched in predicted functions related to cell motility, such as bacterial motility proteins, bacterial chemotaxis and flagellar assembly (Fig. 4B). These features would facilitate bacterial migration towards roots, followed by attachment and establishment 90 . Some root exudates are known to induce motility in bacteria, which is an advantageous characteristic for root colonisation, when compared to non-motile bacteria 91,92 . 16S rRNA gene-predicted functions of streptomycin, novobiocin, tetracycline and vancomycin biosynthesis were significantly enriched in the rhizosphere of semi-dwarf cultivars (Fig. 4B) and this might have some relation to the decrease in keystone taxa which resulted in a more hostile, less coordinated environment compared with tall cultivars, where microbes compete with one another in a less structured manner to gain access to rich nutrients in the rhizosphere. It should be stated that PICRUSt is a tool which provides a functional prediction of microbiome based on marker genes, but it is not an actual measurement of such functions 93 . Future work on shotgun metagenomics would enable to assess these predicted observations.

Conclusions, Limitations and Future Directions
In recent years, there has been a growing interest in modifying plant traits to change the ability of plants to interact with beneficial microbes 13,94 . Collectively, our results showed that wheat breeding from tall to semi-dwarf plants resulted in plants less able to select and sustain a complex rhizosphere. By evaluating differences in root traits of tall and semi-dwarf wheat cultivars and the structure, composition and differential recruitment of specific taxa and keystone species, we can start to understand the impacts of breeding process on root biology. Future work should assess the effect that the breeding process has had on plant hormonal status, exudation profile, plant performance and microbiome selection under low and high nutrient supplementation. This will allow the identification of keystone species and the development of synthetic communities of increased complexity, so the significance of their presence and absence in facilitating the development of a healthy microbiome can be tested.

Data availability
16S rRNA gene amplicon data are available at the NCBI Sequence Read Archive (SRA) under accession number: PRJNA601112.