Rice Plant–Soil Microbiome Interactions Driven by Root and Shoot Biomass

: Plant–soil microbe interactions are complex and affected by many factors including soil type, edaphic conditions, plant genotype and phenotype, and developmental stage. The rice rhizosphere microbial community composition of nine recombinant inbred lines (RILs) and their parents, Francis and Rondo, segregating for root and shoot biomass, was determined using metagenomic sequencing as a means to examine how biomass phenotype inﬂuences the rhizosphere community. Two plant developmental stages were studied, heading and physiological maturity, based on root and shoot biomass growth patterns across the selected genotypes. We used partial least squares (PLS) regression analysis to examine plant trait-driven microbial populations and identify microbial species, functions, and genes corresponding to root and shoot biomass as well as developmental stage patterns. Species identiﬁed correlated with increases in either root or shoot biomass were widely present in soil and included species involved in nitrogen cycling ( Anaeromyxobacter spp.) and methane production ( Methanocella avoryzae ), as well as known endophytes ( Bradyrhizobium spp.). Additionally, PLS of species, functions, and genes. Many of the community functions and genes observed during the heading stage were representative of cell growth (e.g., carbohydrate and nitrogen metabolism), while functions correlated with physiological maturity were indicative of cell decay. These results are consistent with the hypothesis that microbial communities exist whose metabolic and gene functions correspond to plant biomass traits.


Introduction
Soil microbial communities can increase nutrient availability to plants and influence plant growth and overall health [1]. In turn, rhizosphere soil microorganisms rely on root exudates, such as carbon metabolites and other nutrients, as growth substrates. Plant species directly influence soil microbial communities through these root exudates which change as the plant matures [2,3]. Even within a plant species, the rhizosphere soil microbial community can be altered by plant genotypic differences [4,5]. Plant breeding efforts have the potential to make use of beneficial plant-soil microbiome interactions to increase the health and productivity of a crop [6]. However, more needs to be learned regarding how plant developmental stages and their physiological traits influence soil microbial communities before this plant breeding potential can be realized. Rice is a staple crop for half the world's population; thus understanding rice-soil microbial community interactions is important [7][8][9][10][11][12][13][14][15].
Previous studies exploring plant-soil microbiome interactions report plant developmental stage, soil type, and genotype influence the rhizosphere soil microbial community [2,16,17]. It is established that plant developmental stage influences rhizosphere microbial community structure [4,5,8,18] in part, as root exudation patterns change with plant age [3]. Plant developmental impacts on the rhizosphere soil microbiome have been studied for rice [8,13,14,19] and other crops including potato [4,5,20] and maize [3]. Developmental stages from vegetative to maturity result in changes in numerous plant physiological traits. For example, roots and young leaves are major carbon and nutrient sinks during the vegetative stage whereas, after anthesis (i.e., heading) and during the subsequent grain fill stages, fruits and seeds become the dominant metabolic sinks as plants reallocate carbon to these reproductive tissues [21]. The transition that occurs between heading and physiological maturity is, therefore, a critical period to study in order to capture the interactions between root and shoot biomass growth and soil microbial community composition.
Likewise, plant genotypic variation has been shown to modulate the soil microbial community structure in rice [7][8][9]22] and other crops such as maize [23][24][25]. Many studies examined the plant genotype effect on soil microbe interactions using one or a few genetically diverse cultivars, whereas other studies used plants possessing a single gene mutation that modifies a trait of interest. A study by Zhang et al. (2019) examined the effect of mutating NRT1.1B, a nitrate transporter in rice, on the soil microbial community composition compared to that of the wild type. They found the single gene mutation in nrt1.1b decreased the relative abundance of nitrogen-cycling microbial populations compared to the wild-type [26].
Although it is well known that plants can shape the rhizosphere soil microbial community structure through qualitative traits such as root exudates, it remains largely elusive whether quantifiable traits, such as shoot biomass, play a similar role. However, a study by van der Heijden (1990) showed that both root and shoot biomass were positively correlated with mycorrhizal fungal community diversity in grassland species [27]. More recent studies have examined the connection between aboveground plant quantitative traits and belowground soil microbial communities [28,29].
We performed this study with the globally important rice crop, to understand how changes in growth patterns of plant traits, such as root and shoot biomass, during the developmental transition from heading to ripening stages impact plant trait-driven shifts in microbial community structure. This study used nine rice recombinant inbred lines (RILs), and their parents, differing for root and shoot biomass to examine the impact of these traits, measured at heading and physiological maturity, on the rice rhizosphere microbial community structure. RILs represent a genetic composition of the genomes from two parents generally having contrasting phenotypes. They are valuable for studying the impact of differences in plant quantitative traits on the soil microbial community because of their shared genetic background.
Many soil microbiome studies use exploratory methods, such as principal coordinate analyses (PCoA) among others, to characterize the observed shifts in the overall microbial communities between treatments [30]. However, these methods may obscure specific taxa or genes that significantly affect growth but are present as a small proportion of the overall microbial community. For example, it has been observed that relatively high abundance taxa are found across all rice growth stages but low abundance taxa can be characteristic of specific developmental stages [19] and this principle may apply to other plant traits. Partial least squares (PLS) analysis is a statistical tool that can discriminate between treatment groups by a trait of interest. Several soil microbiome studies have used PLS methods to differentiate and group microbial taxa by treatment [31][32][33]. Here we used nine rice recombinant inbred lines, and their parents, segregating for root and shoot biomass at two important developmental stages, i.e., heading and ripening, to understand how changes in microbial taxa, functions, or genes are related to plant biomass traits and developmental changes using both PCoA and PLS regression analysis.

Rice Genotype Selection and Experimental Design
All studies were conducted in fields or greenhouses at the Dale Bumpers National Rice Research Center (DBNRRC) in Stuttgart, AR in 2017. From a Francis (FRCS) and Rondo (ROND) (FR) mapping population of 217 recombinant inbred lines (RILs) grown to maturity in the field, 62 RILs were chosen based on similarity in developmental and physiological stages, i.e., heading date similarity (within 10 days) and plant height (within 30 cm) for a preliminary greenhouse study. The 62 selected RILs were grown for six weeks to evaluate their shoot and root biomass traits ( Figure S1). Seeds were planted (roughly 5 seeds per pot) in small pots (7.6 cm square and 10.2 cm high). Pots used for transplanting were filled with soil collected from a field with a history of rice cultivation at the DBNRRC in Stuttgart, AR. The soil is characterized as a Dewitt silt loam soil (fine, smectite, thermic, Typic Albaqualf), slightly acidic (pH of 5.6), with total C and N contents of 0.66% and 0.085%, respectively. Prior to soil collection for the greenhouse study, basal fertilizer of P (29 kg P ha −1 ) and K (84 kg K ha −1 ) was applied to the field and incorporated. Flooded conditions were maintained in the pots beginning at transplanting and continued for 6-weeks with a water depth of 5-7 cm above the soil surface.
Based on this preliminary study, nine RILs and their parents (Francis and Rondo) were selected which displayed a range in shoot and root biomass. Three replicates of the nine selected FR-RILs (except for RIL 9 which only had 2 replicates), and the two parents were grown in a greenhouse (plant date: 4 April 2017) using a completely randomized design, as described above, and after emergence, seedlings were thinned to one per pot. Four weeks after sowing, seedlings were transplanted into larger pots (27 l volume) with one seedling per pot and pots were placed at a uniform distance from each other. Additionally, two pots, containing soil but no plant, were used as replicated soil controls. Flooded conditions were maintained in the pots beginning at transplanting and continuing to maturity with a water depth of 5-7 cm above the soil surface. Each pot received nitrogen fertilizer in the form of urea (80 kg N ha −1 ) applied in a three-way split during the two months following transplanting. Total photosynthetic active radiation (PAR) fluxes in the greenhouse during the study were 2398.3 mmol m −2 and the average air temperature in the greenhouse was 24.9 ± 0.14 • C.

Plant and Soil Sampling and DNA Extraction
Plants were destructively sampled for biomass traits at the reproductive stage (heading, R4-R5), and ripening stage (maturity, R8-9) [34]. Plant samples and their associated rhizosphere soil samples were taken based on plant developmental stages as the selected FR-RILs matured at different rates. The entire plant was removed from the pot, and soil samples were taken from the center of the soil-root aggregate at a depth between 13-18 cm below the crown, where roots were most abundant. The soil-root aggregate was shaken to remove loose soil, leaving roughly a 1-2 cm soil layer on the roots. Approximately 5 g of soil covered roots were placed in a sterile flask with 50 mL of sterilized phosphate-buffered saline (PBS) solution and shaken gently to wash soil from the roots. Sterile forceps were used to remove roots from the soil slurry solution which was subsequently used for DNA extraction (see below). From the remaining soil-root aggregate, the soil was washed off. Plant tissue was divided into root biomass (RB) (below the crown) and shoot biomass (SB) (crown and above). Plant tissues were weighed after drying at 60 • C for 4-5 days to determine root and shoot dry weight.
Soil DNA was extracted from samples according to Kepler et al. (2018). The DNeasy PowerSoil HTP 96 Kit (Qiagen, Germantown, MD, USA) was used to extract gDNA from 800 µl of soil slurry according to manufacturer's instructions [35]. DNA was quantified using a Qubit dsDNA high sensitivity assay (Invitrogen, Waltham, MA, USA). All samples were diluted to 2.5 ng µL −1 using deionized, autoclaved water for sequencing preparation. A subsample of soil was taken for total C and N analysis using a TruSpec CN analyzer (LECO Corp., Saint Joseph, MI, USA).

Shotgun Metagenomic Library Construction and Illumina Sequencing
Shotgun metagenomic sequencing was conducted in order to capture both the taxonomy and functions of the microbial community associated with SB, RB, and developmental stages heading and maturity. Diluted DNA samples were used for the shotgun metagenomic library construction using the Nextera XT DNA Library Preparation Kit (96 samples) and the Nextera XT Index Kit v2 Set A (96 indexes) (Illumina Inc., www.illumina.com, accessed on 14 March 2021) following the Illumina protocol. Briefly, DNA was fragmented and tagged with adapter sequences. Tagmented DNA was then amplified to add indexed adapters. After amplification, libraries were cleaned up using AMPure XP beads and quality checked on the Agilent Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Libraries were normalized to 0.5 nM using 10 nM TRIS-HCl and 0.1% Tween 20 (pH 8.5) and pooled. Rhizosphere soil samples were sequenced on the Illumina NextSeq 500 platform at the USDA Agricultural Research Service, Beltsville Area Research Center, Beltsville, MD, USA.

Sequence Processing
Reads were trimmed of adaptors and contaminant reads were removed using BBMap version 37.66 (https://jgi.doe.gov/data-and-tools/bbtools/bb-tools-user-guide/bbdukguide/, accessed on 1 February 2018). BBMap was then used to filter reads with kmers matching to the reference rice genome. DIAMOND (v0.9.17.118) [36] was used to perform searches of the reads against the NCBI NR database retrieved (February 2018) in protein space. Taxonomic and functional gene assignments to the NCBI taxonomy database (November 2018) and SEED DB, a functional gene database, (May 2015) [37,38], respectively, were made using the Megan 6 (version 6.15.2) command line tool "daa-meganizer" with the weighted lowest common ancestor (LCA) algorithm using a top percent setting of 3.0. Read counts from the analysis were exported from Megan Community Edition [39]. The data produced by sequencing is compositional count data [40]. For that reason, samples were closed, multiplicative replacement was run, and a centered log-ratio (CLR) transformation was performed using Skbio 0.5.5 (https://scikit-bio.org, accessed on 1 November 2018).

Community and Multivariate Statistical Analyses (PCoA and PLS)
Principal coordinate analysis (PCoA), using Bray-Curtis distance, was used to examine the soil microbial community structure as a whole using Megan Community Edition. Alpha diversity was measured by the Shannon index and calculations were performed using R (version 3.5.1) and Welch's t-tests were conducted to test significant differences at p = 0.05 between developmental stages by genotype [41]. Partial least squares (PLS) regression analysis was performed on the whole soil microbial community to differentiate species, metabolic pathways, and genes as a function of SB or RB dry weight. All PLS regression analysis steps were carried out in JMP v. 14.2.0 (2018 SAS Institute Inc., Cary, NC, USA).
As both PCoA and PLS regression analyses were performed on the entire microbial community, specific species, functions, or gene responses may have been masked by more abundant taxa, functions, or genes. For this reason, we performed a linear least squares fit analysis for each, individual dependent variable (species, function, or gene) as a function of SB or RB. Candidate lists were created for further analyses based on analysis of variance (ANOVA) significance tests (p < 0.05) of the linear least squares fit. These p-values were used only for creating candidate lists, and not used for multiple comparisons between species.
PLS analysis was performed on the selected candidate list of variables (species, function, or gene level), to differentiate their patterns with traits, i.e., SB, RB, and developmental stages (heading and physiological maturity), and to further identify which variable(s) significantly contribute to trait differences. PLS regression analyses can be used to discriminate between treatment groups by a trait of interest and are particularly useful when the number of independent variables (e.g., number of microbial taxa or genes) are significantly greater than the number of response variables (e.g., traits of interest). For the PLS analysis, SB and RB were used as continuous variables and developmental stage (heading or physiological maturity) was used as pseudo-continuous variable by coding samples as either members (1) or non-members (0) for the specific developmental stage [42]. Regression beta coefficients were calculated and used to present the correlation patterns of traits with the selected variables. Developmental stages coefficient values are presented alongside SB and RB coefficients, to visualize their responses, however because they are compared to each other in the model (i.e., members vs. non-members), the coefficient values display equal and opposite patterns. Variable importance in projection (VIP) values were used to approximate trait variation. VIP scores greater than 0.8 are generally considered significant [43,44] and this cut-off was used in our study. Hierarchical cluster analyses (using the Ward method) based on beta coefficient data were used only for ease of visualization to identify groups of variables showing similar trends with traits of interest. These clusters were then used to create coefficient plots which grouped variables by their similarity in response to our traits of interest.
Following selection of metabolic functions related to SB and RB through PLS analyses, a regression analysis of selected dependent variables as a function of plant biomass traits (i.e., RB or SB) meeting a p < 0.05 were examined using CLR transformed count data. This was done to confirm PLS results and examine the correlation between microbial functions and SB and/or RB.

Selection of Recombinant Inbred Lines Segregating for Root and Shoot Biomass
From a Francis and Rondo (FR) mapping population of 217 recombinant inbred lines (RILs), 62 were chosen based on similar developmental and physiological stages, i.e., heading date similarity (within 10 days, Figure S1A) and plant height (within 30 cm, Figure S1B). From this candidate list, nine Francis (FRCS) and Rondo (ROND) recombinant inbred lines (FR-RILs), along with their two parents, were selected for the soil microbiome study. The selected genotypes resulted in a range of low, intermediate, and high shoot biomass (SB) and root biomass (RB) phenotypes (range of 29.3-188.5 g dry weight and 4.6-46.3 g dry weight for shoot and root biomass, respectively, Figure S2) at heading and maturity. In general, root biomass decreased between heading and harvest maturity (harvest maturity was considered to be R7-9) while shoot biomass continued to increase at this stage ( Figure S2).

Whole Microbial Community Structure and Impact of Shoot and Root Biomass Traits
Principal coordinate analysis (PCoA) plots demonstrated that microbial community composition varies by plant genotype and developmental stage. Plant genotype influences microbial community structure at both heading ( Figure 1A) and at maturity ( Figure 1B). Francis and Rondo have distinct microbial communities at heading and maturity differentiated in PC1, with the exception of one Rondo replicate that groups close to Francis at maturity. Taxa identified in both developmental stages include the genera Sphingomonas, Bradyrhizobium, and Anaeromyxobacter and the phyla Acidobacteria and Chloroflexi. Both shoot biomass and root biomass influenced the community structure during heading and at maturity. In contrast, soil percent total C and N did not influence the community structure strongly likely because it did not vary much between samples (mean total C and N 0.73% ± 0.02% and 0.083% ± 0.002%, respectively). Alpha diversity, measured by the Shannon index, showed sample variance by genotype and developmental stage. The mean alpha diversity of each genotype tended to decrease from heading to maturity for RIL1, RIL5, RIL6, and Francis, increase for RIL2, RIL4, and RIL7 but did not change for RIL3, RIL8, RIL9, Rondo, and the soil control samples ( Figure S3). These differences by developmental stage were significant for Francis and RIL4 (p = 0.128 and p = 0.0322, respectively). Additionally, the top 1% of phyla ( Figure S4) and top 2% of species ( Figure  S5) by genotype were plotted to examine their overlap with those species influencing the microbial community structure in the PCoA. The relative abundance plots revealed that all taxa influencing the community structure were highly abundant at both heading and maturity stages across all genotypes.

Shoot and Root Biomass Driven Species Level Analysis
We used PLS to identify taxa, functions and genes that are related to changes in SB and RB at heading and physiological maturity. To focus our analysis on those species, functions, or genes most affected by SB and RB, we used a linear least squares regression analysis to identify those variables from within the larger microbial population. The species, functions, or genes identified as correlated to either biomass trait and passing an ANOVA significance test of p < 0.05 were then used in the PLS regression analyses to further distinguish which displayed differential responses by SB and/or RB as well as their patterns by developmental stages. Regression beta coefficients and variable importance in projection (VIP) scores were plotted to visualize microbial species, function, or gene patterns with SB, RB, and/or the two developmental stages. Hierarchical clustering was used to visualize groups of species, functions, or genes based on their response patterns for SB, RB and the two developmental stages.
Approximately 100 candidate species were identified through ANOVA significance testing as a function of SB and RB (Table 1). PLS beta coefficient values were clustered for visualization, based on their SB, RB, and two developmental stage patterns, resulting in 7 species clusters (SC) ordered by increasing magnitude from SC1-SC7 (Table 1 and Figure 2). The greater the difference between coefficient values of a dependent variable, i.e., species, the greater the contribution of that variable to the model differentiating traits of interest. Additionally, variable importance in projection (VIP) values are presented and are used to approximate trait variation. VIP scores greater than 0.8 were considered significant [43,44].
Species in SC1 through SC4 had a higher abundance at heading when SB was lower as compared to maturity when SB was greater. In contrast, taxa in SC5 to SC7 showed the opposite pattern and were more abundant at maturity. The response due to changes in RB were the same direction as SB but were comparatively small (Table 1 and Figure 2). Species clusters SC4 and SC7 contain the species most significantly differentiated by biomass traits and developmental stages as demonstrated by the comparatively large coefficient values and the VIP scores greater than 1. Species which increase as both shoot and root biomass increase, include several bacterial species from the Planctomycetes, Proteobactera (Alphaproteobacteria), and Verrucomicrobia phyla, and two archaea from the phyla Euryarchaeota and Thaumarchaeota (SC7 in Figure 2). In SC4, species displaying a negative correlation with both SB and RB include mostly those in the Actinobacteria, Chloroflexi, and Proteobacteria (Alpha-and Gamma-proteobacteria) phyla. Table 1. Species clusters (SC1-SC7) significantly correlated with shoot biomass (SB) or root biomass (RB) across two developmental stages. Species names, as well as numerical species ID used for coefficient plots, are shown. Beta coefficient values are reported for SB and RB. Variable importance in projection (VIP) scores are shown, VIP scores greater than 0.8 are considered significantly different by responses, i.e., biomass traits and developmental stages, in the model. All species in clusters 4 and 7 have VIP scores greater than 1.

Species Cluster
Species Name

Rhizosphere Soil Microbial Community Function Analysis
At the functional level, three distinct clusters were associated with SB and RB at the two developmental stages ( Figure 3A). To confirm observations made with PLS analyses, CLR transformed read counts for one to two representative microbial functions from each of the three clusters are shown that were significantly correlated to SB, RB, or both ( Figure 3B-F). All other metabolic functions were plotted to confirm results and are presented in Supplementary Materials ( Figure S6). All microbial functions from each of the three clusters exhibited the same pattern as the representative function shown in Figure 3B-F. Microbial functional cluster 1 (FC1) was divided into those functions (depicted by hashed bars in Figure 3A) displaying significant negative correlations with shoot biomass and significant positive correlations with root biomass (p < 0.05) ( Figure 3B), and those functions with a negative correlation to shoot biomass (p < 0.05) and no significant relationship to root biomass (p > 0.05) ( Figure 3C). CLR transformed functional count data was plotted as a function of SB and RB to verify the correlations found by PLS analyses for fatty acid synthesis ( Figure 3B), metabolism of aromatic compounds, amino acids and derivatives, metabolite damage and repair, cofactor, vitamins, prosthetics groups, and pigment ( Figure S6A). The remaining four functions identified in FC1, had VIP scores greater than 0.8 and showed significant negative correlations with shoot biomass, but weak or no correlation to root biomass including carbohydrate metabolism ( Figure 3C and Figure S6B). FC2 contains metabolic functions considered of lesser importance to the model (VIP < 0.8 for all) compared to all functions in FC1 and three functions in FC3. Among these FC2 functions are stress response and cell division and cell cycle functions ( Figure 3D). In FC3, dormancy and sporulation are significantly positively correlated with shoot biomass but negatively correlated with root biomass ( Figure 3E). All other functions in FC3 display the same significant positive correlation with SB. However, this was not the case for RB as the negative correlation pattern predicted from the PLS model was not present when CLR-transformed count data was plotted as a function of RB, indicating no correlation with RB. These functions include virulence, disease, and defense ( Figure 3F) and motility and chemotaxis ( Figure S6C).  Interestingly, microbial community functions correlated with root biomass were also related to heading stage, while shoot biomass and maturity stage displayed similar functional trends. However, functions correlated with RB and heading showed opposite trends to those correlated with SB and maturity. This opposing trend of RB and SB associated genes was observed for all functions in the PLS-selected candidate list.

Metagenomic Gene Level Analysis With PLS
The PLS analysis used in this study was particularly powerful when examining the more than 7000 microbial gene annotations, to discover genes differentially correlated with plant biomass traits and developmental stage. From the thousands of annotated genes, PLS analysis revealed 249 genes that were significantly associated with both SB and RB ( Figure 4 and Table S1). Genes displaying particularly contrasting patterns in SB and RB, and VIP scores greater than 0.8, are summarized in Table 2. Genes grouped into 11 clusters that displayed similar gene coefficient patterns with RB or SB (Figure 4).  , and virulence, defense, and disease (F). Contour lines in-dicate quantile density at 5% intervals. This means that approximately 5% of points generated from the estimated non-parametric distribution are below the lowest contour, 10% are below the next contour, etc. The highest contour has about 95% of the points below it.
Genes significantly positively correlated with RB and heading (and negatively correlated with SB and maturity) are largely found in gene clusters 1 through 5 (GC1-GC5). The most frequent genes identified in GC1 and GC2 were involved in amino acid metabolism, carbohydrate metabolism of mannose and maltose, fatty acid synthesis, and several stress response genes. Additionally, genes involved in iron acquisition by siderophore biosynthesis, polysaccharide synthesis, nitrogen metabolism genes, cell signaling genes, and metabolism of aromatic compounds, were all found to be negatively correlated with SB/maturity and positively correlated with RB/heading.    GC6 to GC11 represent genes positively correlated with SB and maturity and negatively correlated with RB and heading. Clusters GC8, GC10 and GC11 include virulence related genes involved in conjugative signaling and transfer and genes involved in potassium uptake, photosynthesis, and nucleotide synthesis. Additionally, several genes related to spore formulation are strongly related to increasing SB and decreasing RB. Many genes related to protein metabolism, specifically large subunit ribosomal proteins, and amino acid synthesis and degradation functions (including Acetyl-CoA acetyltransferase involved in lysine degradation) were also identified. Interestingly, the gene level patterns observed in GC8, GC10, and GC11 mirror the patterns observed in the functional level cluster 3 (FC3). PLS analyses were likewise carried out for genes specifically associated with SB or RB, individually and can be found in Supplementary Materials (Figures S7 and S8, respectively).

Discussion
A number of studies have shown that microbial community structure changes with plant genotype and developmental stage and this trend was observed in our study in which microbial communities differed across the FR-RILs and shifted between heading and maturity stages within genotypes [3,8,17]. Although we did observe genotypic differences, they were not based on the observed SB and RB trend that we selected for (i.e., two parents, Francis and Rondo, displaying contrasting SB and RB traits and their nine FR-RILs, encompassing a range of SB and RB) when the soil microbial community was examined as a whole.

Soil Microbial Populations Associated With Rice Shoot and Root Biomass
To identify specific taxa, microbial functions, and genes related to biomass traits of interest, we used two multivariate analysis tools: PCoA and PLS regression analysis. The PCoA was conducted first to examine the whole microbial community shifts in the heading and physiological maturity stages. RB and SB explained some of the rhizosphere community structure in both stages, with SB being dominant over RB in both stages. Several taxa explained the microbial community structure in both stages and were also found to be highly abundant (Figures S5 and S6) including the genera Bradyrhizobium and Anaeromyxobacter. However, PCoA did not allow us to interrogate which specific taxa, functions, and genes were correlated with SB and RB, to do this we used PLS regression analysis.
PLS regression analyses were performed on the entire microbial community to identify microbial species correlated with SB and RB. However, when linear regressions of species CLR-transformed count data were plotted against SB or RB to confirm the model, many of the coefficient patterns identified through the PLS model showed a very weak, nonsignificant correlation (p > 0.05). We speculated that the presence of abundant microbial taxa that do not necessarily respond to SB and RB may have obscured changes in community structure of SB or RB associated species. For this reason, we used linear least squares regression analyses to identify species, functions, and genes that are associated with SB or RB for use in the PLS analysis. In this discussion we consider only species showing the most differentiation in PLS results, i.e., species in SC4, SC6, and SC7, with VIP scores greater than 0.8 as they are considered significant effectors in the PLS models [43,44]. Microbial taxa identified through the PLS analysis are all widely present in soil, including species from the Actinobacteria, Chloroflexi, Planctomycetes, Proteobacteria, and Verrucomicrobia phyla. Within the Proteobacteria, several Alphaproteobacteria species were identified as correlated with shoot and root biomass. Within this phylum, a Phenylobacterium sp. (ID = 66, VIP = 1.192), aerobic or facultatively anaerobic non-spore forming bacteria [45,46], was negatively correlated with RB, SB, and maturity. Perhaps unsurprisingly, several Bradyrhizobium spp., which are known rice root endophytes capable of nitrogen-fixation [47], were identified as positively associated with root biomass, but, interestingly, also with shoot biomass and physiological maturity (Table 1 IDs 69, 70, 71 and VIP scores of 1.306, 1.209, 0.983, respectively). A previous study by Hirano et al. (2001), found that rice grown without nitrogen fertilizer, in soils with a high abundance of nitrogen-fixing bacteria in late growth stages, had similar yields to rice grown under conventional fertilization practices [48,49]. This underscores the importance of rice endophytes, including Bradyrhizobium spp., in rice cultivation. The fact that the abundance of these bacteria was found to correlate with both root and shoot biomass patterns with 11 genotypes suggests these traits could be used for breeding to select for nitrogen-fixing microbial populations.
Several species showed a significant positive correlation with rice shoot biomass and this pattern was well-presented at maturity. Especially interesting, Anaeromyxobacter spp. (Table 1 IDs 118, 119, 120 and VIP scores of 1.008, 1.182, 0.913, respectively) are Gram negative facultative anaerobes that are known to perform dissimilatory iron reduction [50] and have been shown to use nitrate [50] as an electron acceptor making it an important species for biogeochemical cycling of nutrients. A methanogen, Methanocella avoryzae (ID 279, VIP = 1.427), was also positively correlated with increased shoot biomass and physiological maturity. Peak methane emissions are typically observed during heading stage [51], however, higher shoot biomass likely creates a larger conduit for methane to diffuse from the soil to the atmosphere through aerenchyma tissue [9,52]. Additionally, Candidatus Nitrosocosmicus oleophilus (ID 282, VIP = 1.276), an ammonia oxidizing archaea (AOA) was found to correlate positively with biomass traits and physiological maturity.
AOA are autotrophs that fix inorganic carbon and are thought to be the dominant microbial population performing ammonia oxidation in soil [53,54]. AOA and ammonia oxidizing bacteria (AOB) are involved in nitrate formation on rice root surfaces which have been shown to be an important nitrogen source for rice grown in submerged conditions [55,56]. They are key players in the nitrification pathway, which oxidize ammonia to nitrite, a nitrogen form less amenable for assimilation by rice. Interestingly, this archaea is negatively correlated during the heading stage, when rice is transitioning from its reproductive phase to the ripening phase and requires more nitrogen in the easily assimilated form of ammonia [22].
There were a number of overlapping taxa between the predominant taxa present in the soil samples ( Figures S4 and S5) and the taxa identified by PCoA ( Figure 1) and PLS (Figure 2) analysis. Methanocella avoryzae and several Anaeromyxobacter species were among those identified as having relatively high abundance ( Figure S5) and positively correlated with shoot and root biomass ( Figure 2 and Table 1). However, other species, including Bradyrhizobium spp. and Candidatus Nitrosocosmicus oleophilus, were not in the top 2% of taxa but were significantly positively correlated with higher biomass. Species within Bradyrhizobium and Anaeromyxobacter influenced the community structure at both heading and maturity (Figure 1) with the latter species having relatively high abundance (in the top 2%) across genotypes and developmental stages ( Figure S5). Other high abundance taxa including Sphingomonas spp., and unclassified taxa within the Actinobacteria, Chloroflexi, Planctomycetes, Proteobacteria, and Verrucomicrobia phyla, likewise influenced the PCoA structure (Figure 1 and Figure S4). Thus, while the PCoA did identify some of the same species as PLS, these results underscore the importance of the PLS analysis which specifically identified those species positively correlated with biomass traits regardless of abundance.
PLS analysis was able to differentiate microbial species based on root biomass; however, coefficient values for root biomass were one to two orders of magnitude smaller than those of shoot biomass, indicating RB had less impact in the species level model than SB. The nine selected RILs and their parents exhibited a lower and wider range in root biomass (roughly 10-fold difference) as compared to shoot biomass (roughly 6.5-fold difference). This likely allowed for stronger trends to arise from PLS analysis with SB data, which surpasses any relationships seen with RB. Additionally, processing soil samples for root quantification is difficult and time-consuming which may lead to underestimation of plant root biomass [57]. Despite the potential for underestimation in RB quantification, PLS analyses were still able to identify functions related to RB. This demonstrates that PLS analyses, used in conjunction with ANOVA significance selection from the microbial community by traits of interest, have the potential to be used to identify trait-associated microbial populations and to understand the relationship between traits of interest.
Overall, through PLS analysis, we were able to identify species that were abundant with lines that had a lower range of root and shoot biomass at heading, and other species that were more abundant with lines that had a higher range of root and shoot biomass at maturity. Our results also indicate that microbial populations are influenced by RB and SB in the same manner. This suggests that shoot biomass has potential to be used as a proxy for root biomass in breeding to select for beneficial soil microbial species.

Microbial Community Functions Correlated to Biomass Traits
Microbial community functions identified in our analysis were distinct between the two developmental stages, heading and physiological maturity. Heading stage represents a dynamic phase of plant growth as rice transitions from reproductive to ripening stages. This growth phase was reflected in the microbial functions related to cell growth (synthesis of amino acids and fatty acids) and metabolism (carbohydrate and aromatic compound metabolism) in our study. In contrast, many of the microbial functions indicative of cell growth plateauing or decaying, such as sporulation, virulence defense, and motility and chemotaxis, were represented at maturity. Chemotaxis functions may indicate bacterial populations are searching for nutrients as root exudates decrease with plant maturation and they are required to search further out to find growth substrates [3,58].
Shoot and root biomass related functions followed the same trends as developmental stages maturity and heading, respectively. Root biomass related functions showed an opposing pattern to SB related functions. Although differences in coefficient values for RB were smaller than those of SB at the species level, they were greater at the functional level, indicating the PLS model was able to discern microbial functions differentially correlated with RB. This suggests that RB may more strongly influence overall community function than it does the abundance of specific taxa.

Gene Trends Related to Shoot and Root Biomass
Many of the genes correlated with RB or SB belonged to the pathways identified at the functional level. However, in contrast to the functional level analysis, RB showed a stronger pattern with gene abundance than SB. This suggests that PLS analyses can make use of relatively small variation in observed traits to differentiate changes in abundance among thousands of genes and identify those related to traits of interest. Even though only a few RB-associated species and functions were identified with PLS as compared to SB, it appears RB may be a more dominant driver of soil microbial community gene abundance than SB. This implies there may be functional redundancy across species; thus, although individual species or functions do not correlate with the RB phenotype, specific genes do.
As expected, the genes most frequently identified were involved in amino acid, carbohydrate, and protein metabolism, all processes central to microbial life. Genes identified as specifically positively associated with RB and heading include those involved in amino acid degradation including an aminotransferase. Additionally, two genes for mannose and maltose metabolism as well as a 2-ketogluconate kinase, part of the pentose phosphate metabolic pathway, were identified as positively associated with RB. A gene related to nitrogen metabolism by nitrite transport showed a similar positive correlation with RB and heading which indicates bacterial nitrogen cycling is prominent during the heading stage. Fatty acid synthesis genes and several stress response genes were likewise positively correlated with RB. Genes involved in iron acquisition by siderophore biosynthesis, polysaccharide synthesis, cell signaling, and metabolism of the phenolic compound vanillic acid, found in plant root exudates [59], are all representative of cell growth, and were found to be positively correlated with RB/heading and negatively correlated with SB/maturity.
Genes that were identified as being positively correlated with SB and maturity and strongly negatively correlated with RB and heading belong to many of the cell decaying functions observed at the functional level. Virulence related genes, for example, were identified as uniquely related to RB and maturity, particularly genes predicted to be involved in conjugative signaling and transfer, indicate cells are undergoing attack. Several genes involved in spore formulation are strongly related to increasing SB and decreasing RB, suggesting microbial populations are preparing for dormancy as plants mature from heading to maturity and RB begins to decrease. Two photosynthesis-related genes were identified as positively correlated with SB and maturity and negatively correlated with RB and heading. As roots decay and root exudates decrease with age, a decrease in available carbon metabolites that bacteria rely on may cause the soil microbial community to favor autotrophs capable of photosynthesis resulting in a negative correlation between photosynthesis-related genes and root biomass.

Relationship of Developmental Stage to Microbial Community Structure and Functions
In general, beta coefficient patterns during heading stage followed those of SB while coefficient patterns for maturity were more closely associated with RB patterns in our study. During heading stage, plants experience rapid growth whereas during maturation stages rice is remobilizing resources to grain development, and as a consequence, belowground carbon allocation decreases [60] as do overall root exudates [61]. As plants age, root exudates change and consequently modulate the microbial community structure and function. In our study, FR-RIL SB increased from heading through the grain fill stages while RB decreased between these stages. Thus, up until maturity, SB is in a growth phase, influencing, and likewise being influenced by, the microbial community. At maturity, FR-RIL root biomass decreased relative to heading. Since RB is in direct contact with the soil microbial community, it follows that during RB growth (i.e., heading stage), when root exudates are high, it is driving the microbial community structure and function.
It is difficult to disentangle the effect of biomass traits from developmental stage, thus focusing on two developmental stages, heading and physiological maturity, that show contrasting biomass traits was a key component of this study. However, further studies validating the results herein are needed. For example, inoculation studies using PLS-identified species or species harboring PLS-identified genes could be used to verify if they positively impact RB and SB growth during the critical heading and maturing stages. A study by Mayer et al. (2019), did just this, using identified endophyte isolates to inoculate several plant species and evaluating their impact on root and shoot growth [62]. Future work in this vein could prove useful for validation, following identification of specific taxa or genes related to RB and SB growth.

Conclusions
Overall, these results suggest species, functions, and gene abundance patterns vary by shoot and root biomass, and across developmental stages of heading and maturity. This study indicates changes in root and shoot biomass may play an equally important role in driving the microbial community as does developmental stage. However, our results suggest that SB may be used as a surrogate for root biomass in future plant traitbased breeding, as we found that rhizosphere microbial species vary with SB and RB in the same manner. Using SB as a breeding trait would be advantageous over RB, as RB quantification is time-consuming and prone to error, while SB quantification is more efficient and accurate. This study underscores the potential of exploiting rice phenotypic variation in plant breeding to promote beneficial plant-soil microbiome interactions.
Supplementary Materials: The following are available online at https://www.mdpi.com/1424-2 818/13/3/125/s1, Figure S1: Rice genotype selection based on shoot and root biomass, Figure S2: Biomass trends of nine rice genotypes, Figure S3: Shannon alpha diversity, Figure S4: Relative abundance of the top 1% of phyla, Figure S5: Relative abundance of the top 2% of species, Figure S6: Linear regressions of biomass traits and metabolic function abundance, Figure S7: Partial least squares results for genes related to shoot biomass, Figure S8: Partial least squares results for genes related to root biomass, Table S1: Genes correlated to biomass traits.