Microbial communities and diazotrophic activity differ in the root‐zone of Alamo and Dacotah switchgrass feedstocks

Nitrogen (N) bioavailability is a primary limiting nutrient for crop and feedstock productivity. Associative nitrogen fixation (ANF) by diazotrophic bacteria in root‐zone soil microbial communities have been shown to provide significant amounts of N to some tropical grasses, but this potential in switchgrass, a warm‐season, temperate, US native, perennial tallgrass has not been widely studied. ‘Alamo’ and ‘Dacotah’ are cultivars of switchgrass, adapted to the southern and northern regions of the United States, respectively, and offer an opportunity to better describe this plant–bacterial association. The nitrogenase enzyme activity, microbial communities, and amino acid profiles in the root‐zones of the two ecotypes were studied at three different plant growth stages. Differences in the nitrogenase enzyme activity and free soluble amino acid profiles indicated the potential for greater nitrogen fixation in the high productivity Alamo compared with the lower productivity Dacotah. Changes in the amino acid profiles and microbial community structure (rRNA genes) of the root‐zone suggest different plant–bacterial interactions can help to explain differences in nitrogenase activity. PICRUSt analysis revealed functional differences, especially nitrogen metabolism, that supported ecotype differences in root‐zone nitrogenase enzyme activity. It is thought that the greater productivity of Alamo increased the belowground flow of carbon into roots and root‐zone habitats, which in turn support the high energy demands needed to support nitrogen fixation. Further research is thus needed to understand plant ecotype and cultivar trait differences that can be used to breed or genetically modify crop plants to support root‐zone associations with diazotrophs.


Introduction
Associative nitrogen fixation (ANF) by diazotrophic bacteria provides an alternative to chemical fertilizers in support of crop yields. ANF is thought to be a loose form of quid pro quo, whereby energy and carbon fixed by the plant are exchanged for nitrogen fixed by bacterial diazotrophs. ANF has generally been considered a relatively small player in the annual nitrogen economy of temperate grasses; however, in some cases, it may provide up to 35% of nitrogen to agriculturally important nonlegumes and forage grasses, such as Miscanthus, energy cane, and switchgrass (Weier, 1980;Chalk, 1991;Wewalwela, 2014). ANF offers a natural solution to an immediate need for plant-available nitrogen while also reducing the footprint (Smil, 1999a(Smil, ,b, 2002Galloway et al., 2003;Oenema et al., 2009;Sutton et al., 2011a,b) of land applied synthetic chemical fertilizers. Some gaps in knowledge, however, need to be filled before using ANF for sustainable production. The natural variation in ANF that occurs across plant cultivars and plant growth stages will help to identify plant-microbial traits and interactions that underlie potentially high rates of ANF, and can be used to eventually breed and develop cultivars with greater ANF capacities. Research in this regard carried out in a model grass will likely be transferrable to other major worldwide grain crops (e.g., maize, wheat).
Switchgrass (Panicum virgatum L.) is a perennial, US native grass. Its biomass is used for forage for livestock and for bioenergy. Switchgrass is useful in the prevention of soil erosion and provides wildlife habitat and carbon sequestration (Ma, 1999;Skinner, 2009;Follett et al., 2012). Due to the various environmental and monetary benefits of switchgrass, identifying microbes for its sustainable production has been an active area of research (Jesus et al., 2010(Jesus et al., , 2016Mao et al., 2011Mao et al., , 2013Mao et al., , 2014Chaudhary et al., 2012;Hargreaves et al., 2015). These studies have shown presence of bacteria that are capable of nitrogen fixation; however, studies rarely confirm the occurrence and extent of nitrogen fixation associated with switchgrass.
It is known that cultivars of plants often have differences in N requirements, abilities to support nitrogen fixation (Porter, 1966;Day et al., 1975;Boddey & Dobereiner, 1988;Ledgard & Steele, 1992), unique rhizosphere microbiomes (Miller et al., 1989;Li et al., 2014), and thus are important sources of variation in plant-microbial traits. A recent study, for example, identified diverse sets of nitrogen-fixing bacteria associated with roots and shoots of switchgrass (Bahulikar et al., 2014). Other studies focused on targeted approaches of using nitrogen-fixing bacterial endophytes to improve productivity across cultivars (Ker et al., 2012;Xia et al., 2013;Lowman et al., 2016). Plant bacterial associations and how they differ among plant species and cultivars can help to understand the interplay of traits important for defining the drivers of a globally significant plant-bacterial function like that of ANF.
'Alamo' and 'Dacotah' are tetraploid cultivars of switchgrass, adapted to the southern and northern US ecotypes, respectively. Compared with Dacotah, Alamo shows higher biomass productivity, taller shoots, drought tolerance, and disease resistance. Nitrogen is required for the growth and maintenance of plants, so we hypothesized that a higher productivity cultivar of switchgrass, Alamo, has more nitrogen demand than Dacotah and would be associated with a different root-zone bacterial community with greater nitrogenase activity than a lower productivity cultivar, Dacotah. Here, we investigated the nitrogenase activity, amino acid composition, structure, function, and interactions of microbial communities in the root-zones of two switchgrass cultivars across multiple growth stages. The research offers novel insights into the dynamics of nitrogen-fixing communities in the rootzone of switchgrass that can help to explain the variation associated with this globally important plant-bacterial interaction.

Experimental setup
Root-zone soil was collected from the Virginia Tech Agronomy Farm/Urban Horticulture Center field growing Alamo, Dacotah (~7 years), and their F1 progeny lines (~2 years) to serve as an appropriate inoculum of potential microbes associated with switchgrass. Alamo and Dacotah seeds were imbibed and allowed to germinate at 28°C. After approximately 7 days, the germinated seeds were potted in sieved (4.75 mm) and homogenized field soil and grown in the glasshouse. The plants were kept moist, well aerated, and supplemented with light to maintain optimal growth. Sampling was performed in replicates (Table S1) at three time points from imbibition, indicative of important growth stages: stages V0 (~2.5 weeks old) and V2 (~1.5 months old) from the vegetative phase, and stage E3 (~3.5 months old) from the elongation phase (Moore et al., 1991).

Nitrogenase enzyme activity
To estimate the relative rate of nitrogen fixation, nitrogenase activity was evaluated by the acetylene reduction (AR) assay, which measures the quantity of acetylene reduced to ethylene (Weaver & Danso, 1994). Briefly, the day before the nitrogenase measurement, the Alamo and Dacotah containing pots were wetted to saturation and allowed to drain gravimetrically to~0.03 MPa. Pots were randomly sampled to obtain replicates (n = 3 for V0, n = 4 for V2, and n = 5 for E3) (Table S1) that received acetylene ('treatment') or remained untreated ('controls'). For the vegetative growth stages, V0 and V2, the entire pot with soil + plant was placed into the Mason jars for the AR assay. At stage E3, we gently removed the plant from the pot and placed the plant + rhizosphere soil into the Mason jar for the AR assay. The data was normalized as per the weight of the 'system' (soil + plant for V0 and V2; plant + rhizosphere soil for E3). Each replicate sample was sealed within a 500-mL to 1-L Mason glass jar. Calcium carbide (CaC 2 ) and water were mixed to produce acetylene gas, which was injected to create a~10% (v/v) headspace. Mason glass jars with sample but without acetylene served as controls, to assess the natural (background) production of ethylene. At four time points (immediately and then approximately every 30 min), headspace gas was mixed using a 10-mL syringe and needle to extract gas from the container and analyzed for acetylene and ethylene concentrations using a gas chromatograph and column (GS-Carbonplot, Part # 113-3133, Agilent Technologies, Inc., Santa Clara, CA, USA, 30 m 9 0.32 mm 9 3 lm) in line with a flame ionization detector. The ethylene accumulation rate was used to determine nitrogenase enzyme activity.
Known ethylene standards were used to generate a linear calibration curve that was used to calculate ethylene concentrations for each sample and reported as the rate of ethylene production, lg g system À1 h À1 . Following the confirmation of normality, equality of variances, and determination of outliers using Grubb's test (Grubbs, 1969) in GRAPHPAD (http://graphpad.com/quickcalcs/Grubbs1.cfm), the rates of ethylene production between Alamo and Dacotah were compared (P-value <0.05) using two-sample t-test (V0 and E3) and Mann-Whitney U-test (V2). Only the 24 controls (3, 4, and 5 samples each of Alamo and Dacotah from V0, V2, and E3, respectively), samples were used for the following experiments.

Soluble amino acids
Shoots were removed, and the root and the rhizosphere soil were extracted in 10 mL of 0.9% NaCl in a 50-mL falcon tube. The roots were washed to remove any attached soil and used to measure fresh and dry weight. Falcon tubes were vortexed at 200 rpm on electric orbital shaker for 20 min, allowed to sit for 5 min. Eight milliliters of solution was transferred to a new 15-mL falcon tube and centrifuged at 3500 9 g for 15 min. The amount of supernatant (approximately 5 mL) transferred to new smaller size falcon tubes was recorded, and 8 lL of internal standard (2.5 mm a-aminobutyric acid; AABA) was added to each sample. The volume for each sample was made to 10 mL by adding 0.9% NaCl. After shaking for 2 min, 1.5 mL of sample was transferred into a 2-mL microfuge tube through 0.22-lm polyvinylidene fluoride (PVDF; Thermo Scientific TM (Thermo Fisher Scientific, Waltham, MA, USA) Target2 TM Syringe Filters, Cat# 03377155) membrane syringe filter; 500 lL of filtrate from the 2-mL tube was transferred to a new 1.5-mL microfuge tube and dried in a vacufuge (speedvac) for 2.5 h at 60°C (function 3). The dried pellet was derivatized using the AccQ FluorTM reagent kit (fluorescent 6-aminoquinolyl-N-hydroxysuccinimidyl carbamate derivatizing reagent; Waters Corporation, Milford, MA, USA; Cat# WAT052880) following the standard protocol (Bosch et al., 2006;Hou et al., 2009). Chromatographic separation on the HPLC 1260 Infinity system (Agilent Technologies, Inc.,) was carried out on a reversed-phase column (Waters X-Terra MS C18, 3.5 lm, 2.1 m 9 150 mm). The mobile phase consisted of A: an aqueous solution containing 140 mM of sodium acetate, 17 mM of triethylamine (TEA; Thermo Fisher Scientific, Cat# O4884100), and 0.1% (g L À1 , w/v) disodium dihydrogen ethylenediaminetetraacetate dihydrate (EDTA-2Na 2H2O; Sigma-Aldrich, St. Louis, MO, USA, CAS# 6381-92-6), pH 5.05, adjusted with phosphoric acid solution, and B: acetonitrile (ACN; HPLC grade, Thermo Fisher Scientific, Cat# A998-1): ultrapure water (60 : 40, v/v). The gradient conditions were 0-17 min 100-93% A, 17-21 min 93-90% A, 21-30 min 90-70% A, 30-35 min 70% A, 35-36 min 70-0% A, and then hold for 4 min before restoring to the initial composition at 40.5 min, with the final composition kept for 9 min. The column was thermostated at 50°C and operated at a flow rate of 0.35 mL min À1 . The sample injection volume was 5 lL. The analytes detection was carried out using a fluorescence detector (kex = 250 nm and kem = 395 nm) (Bosch et al., 2006;Hou et al., 2009). Soluble amino acids in the samples were qualified and quantified by comparison with amino acid standard solutions. Each amino acid standard solution contained 19 protein amino acids including alanine (Ala), arginine (Arg), aspartic acid (Asp), asparagine (Asn), cystine (Cys-Cys; more stable form in oxidative condition than monomer cysteine), glutamic acid (Glu), glutamine (Gln), glycine (Gly), histidine (His), isoleucine (Ile), leucine (Leu), lysine (Lys), methionine (Met), phenylalanine (Phe), serine (Ser), threonine (Thr), tyrosine (Tyr), tryptophan (Trp), and valine (Val) and 1 nonprotein amino acid ornithine (Orn). Note that Asn and Ser were co-eluted and Cys-Cys and Tyr were coeluted. A total of 20 peaks were quantified. The amino acid data was relativized per sample and used for multivariate data analyses using PC-ORD software version 6.0 (MjM Software, Gleneden Beach, OR, USA).

DNA extraction
Ten milliliters of 0.9% NaCl was added to the falcon tubes containing rhizosphere soil samples. The tubes were then turned upside down once and the soil solution was homogenized by vortexing 1 min each (at level 3). A 3 mL sample of the resuspended solution was taken immediately after vortexing in a 5-mL tube. After centrifugation of the 5-mL tubes at 10 000 g for 10 min at 20°C, the supernatant was removed and the tubes were turned upside down on a Kim wipe (with cap open) to dry for 10 min. For each tube, the leftover soil (pellet) was mixed and 0.25 g was used for the DNA extraction using MoBio PowerSoil â (MoBio Laboratories, Carlsbad, CA, USA) DNA Isolation kit as per the manufacturer's protocol. The quality and concentration of the DNA was checked using 0.8% (w/ v) agarose gel electrophoreses and NanoDrop 2000 spectrophotometry.
Real-time PCR to detect nitrogenase reductase gene copy numbers qPCR assay was performed on an ABI 7300 system (Thermo Fisher Scientific) to quantify the abundance of nitrogenase reductase gene (nifH). The 20-lL reaction volumes contained 10 lL of 2X PowerUp TM SYBR â Green Master Mix (Thermo Fisher Scientific), 2 lL template at 2 ng lL À1 , 4 lL each of PolF/PolR (10 lM) primers (Poly et al., 2001), and standard conditions as per master mix's protocol. Triplicates of nucleasefree water and no-template controls were included. Technical triplicates of each biological replicate were averaged to get mean C T value. DNA was extracted from Sinorhizobium meliloti 1021 (cultures provided by Dr. B. Scharf at Virginia Tech) using UltraClean â Microbial DNA Isolation Kit (Mo Bio). Dilutions of known concentrations of Sinorhizobium meliloti 1021 were used to generate the C T vs. log (N 0 ) standard curve. The standard curve was used to calculate the qPCR efficiency (E = 10 (À1/ slope) ) and nifH gene copy numbers in the samples (Brankatschk et al., 2012). For each cultivar at each growth stage, biological replicates were checked for outliers using Grubb's test (Grubbs, 1969) in GRAPHPAD (http://graphpad.com/quickcalcs/ Grubbs1.cfm) and boxplots were generated in BoxPlotR with Spear's criteria of whisker definition (Spitzer et al., 2014). The cultivars at each growth stage were tested for significant (Pvalue <0.05) differences in nifH gene copy numbers using Mann-Whitney U-test.

Sequencing and data analyses
As bacterial-fungal interactions are known to promote plant growth (Artursson et al., 2006;Gamalero et al., 2008;Frey-Klett et al., 2011), we investigated the bacterial and fungal communities to identify any possible interactions in root-zones of these switchgrass cultivars. The 16S rRNA and ITS gene amplification were performed using bacterial primer pair (S-D-Bact-0341-b-S-17 and S-D-Bact-0785-a-A-21) (Klindworth et al., 2013) and fungal primer pair (ITS1F and ITS2) (Smith & Peay, 2014), respectively. The library preparation was performed using Illumina 16S Metagenomic sequencing library preparation guide. DNA concentration of the pool was determined by fluorometric quantification using the Qubit â 2.0 platform with Qubit dsDNA HS Assay Kit (Life Technologies, Carlsbad, CA, USA). Multiplexed 250-base pair paired-end sequencing using Illumina MiSeq was performed at the Biocomplexity Institute (formerly Virginia Bioinformatics Institute) core facility at Virginia Tech. Data are submitted under Biosamples of SAMN04917354-SAMN04917359 (16S), SAMN04917374-SAMN04917379 (ITS), BioProject PRJNA320123 at NCBI (SRA).
Barcode adapters and primers were trimmed from each read using CUTADAPT v1.8.1 (Martin, 2011) with a quality-trimming threshold of 30 and minimum read length of 100. The pairedend reads were merged based on overlapping sequences into single reads using PANDASEQ v2.8 (Masella et al., 2012). The bacterial and fungal sequencing data were analyzed using QIIME v1.8.0 (Caporaso et al., 2010) as previously described (Rodrigues et al., 2015). Briefly, reads were clustered into OTUs based on 97% sequence similarity using UCLUST (Edgar, 2010) and USE-ARCH61 (Edgar, 2010) for bacteria and fungi, respectively, using an open-reference OTU-picking strategy. The representative sequence of an OTU was used to assign it a taxonomy, using UCLUST against the GREENGENES reference database version 13.8 (Desantis et al., 2006;McDonald et al., 2012) for bacteria, and RDP classifier (Wang et al., 2007) against the UNITE reference database version 12.11 (Abarenkov et al., 2010) for fungi.
Postprocessing included diversity and richness analyses, identifying and summarizing the most abundant taxons, and describing the taxons that are different and indicative of cultivars at different growth stages. Briefly, the alpha diversity was calculated on the OTU table for all samples using several different indices, including PD whole tree (only for bacteria), chao1, observed species, Good's coverage, Shannon, and Simpson. After using a sequence threshold, the beta diversity for bacteria and fungi was calculated using weighted UniFrac (Lozupone et al., 2007(Lozupone et al., , 2011 and Bray-Curtis (Beals, 1984) metrics, respectively, and were used for principal coordinate analysis (Gower, 2005) and visualization. The genus-level summary of communities was used for NMDS analysis using the Bray-Curtis dissimilarity in R (VEGAN). Using the collated alpha diversities from the rarefied OTU table, the PD whole tree (only for bacteria), chao1, and observed species indices were used to compare alpha diversity of groups (cultivar, stage, cultivar 9 stage) using a two-tail nonparametric t-test with FDR correction (q-value <0.05) using collated files. Multivariate data analysis methods of MRPP (Mielke, 1984), ADONIS (Anderson, 2001), and analysis of similarity (ANOSIM) (Clarke, 1993) were used to identify whether the cultivars and plant growth stage had an effect on the microbial communities.
Indicator species analysis (ISA) (Dufrene & Legendre, 1997) was performed in R (VEGAN) using the top 50% most abundant genera to identify taxa significantly (indicator value >50 and Pvalue <0.05) indicative of cultivars at the respective growth stage. PICRUST (Langille et al., 2013) analysis using the bacterial OTU abundance was performed before and after ISA analysis to identify whether functions differed between cultivars at the different growth stages. OTUs not part of the closed reference OTU picking were filtered out, and the metagenomes were collapsed into KEGG pathways. Using Statistical Analysis of Metagenomic Profiles (STAMP) (Parks et al. 2014), principal component analysis was used to visualize differences in KEGG pathways of bacteria between the cultivars at different growth stages. A two-sided t-test was performed to check whether Nmetabolism pathway was significantly different (P-value <0.05) between the cultivars.

Network analysis
The OTU table, at the class level for bacteria and order level for fungi, was summarized, relativized, and separated to create files as per samples belonging to cultivar irrespective of the plant growth stage, that is, Alamo (V2, E3) and Dacotah (V2, E3). Phylogenetic ecological networks were generated for each group (cultivar 9 growth stage) using MENAP (Deng et al., 2012) with an RMT threshold of 0.31 and module detection using leading eigenvector of the community matrix (Newman, 2006) and random walk methods (Pons & Latapy, 2005). A 'module' is a group of OTUs that were highly connected (positively or negatively correlated via abundance) among themselves, but had much fewer connections with OTUs outside the group. For each network, connections between OTUs between different modules were discarded to help focus on interactions among OTUs within a module. We refer these subnetworks as 'same-module networks' which show the connections between OTUs belonging to the same module. For each cultivar, the networks from leading eigenvector of the community matrix and random walk methods were compared to generate edges common between the two methods. These 'COMmon Edge between Two methods' (COMET) networks help to remove any potential biases due to the module detection algorithms and represent high confidence edges. Finally, the COMET networks were compared between cultivars.

Cultivable, free-living diazotrophs
For samples from stage V2, the NaCl solutions containing rhizosphere soil were (100x, 1000x, and 5000x) serial diluted, spread on yeast-mannitol (YM) agar (YMA) petri plates, and incubated at 28°C for 3-5 days to selectively grow diazotrophs. Bacterial colonies were differentiated based on morphology and color. Colonies were counted and will be reported as colony forming unit per gram of root-soil. Representative colonies were sampled and isolation plated on new YMA plates to obtain pure cultures. After 3-5 days of incubation, the isolation-plated samples were grown in yeast-mannitol broth and DNA was extracted using UltraClean â Microbial DNA Isolation Kit (Mo Bio). The quality of the DNA was estimated using 0.8% (w/v) agarose gel electrophoreses and NanoDrop 2000 spectrophotometry. The 16S rRNA genes were amplified using polymerase chain reaction using primers 27F (5 0 -A GAGTTTGATCMTGGCTCAG-3 0 ) and 1492R (5 0 -TACGGYT ACCTTGTTACGACTT-3 0 ) (Lane, 1991;Fredriksson et al., 2013) and KAPA2G Robust PCR kit (Kapa Biosystem, Wilmington, MA, USA), and cleaned using Agencourt AMPure XP PCR product cleanup kit (Beckman Coulter, Brea, CA, USA). Highquality, cleaned PCR products, as determined by 1.2% (w/v) agarose gel electrophorese and NanoDrop 2000 spectrophotometry, were sent to Biocomplexity Institute (Virginia Bioinformatics Institute) for Sanger sequencing. The 'ab1' files obtained from Sanger sequencing were converted to sequence files using 4Peaks and quality trimmed with a Phred score of 25. We used the taxonomy of the top hit from nucleotide BLAST (MEGABLAST) on the sequences to identify the bacteria. LIBSHUFF (Singleton et al., 2001;Schloss et al., 2004) analysis on the sequences and Pearson's chi-squared test on the genus-level CFU abundance were performed to identify significant differences (P-value <0.05) between the cultivable diazotrophic communities associated with Alamo and Dacotah.

Results
Alamo showed higher rates of biological nitrogen fixation (BNF) Quality control analysis using the Grubb's test (Grubbs, 1969) in GRAPHPAD identified an outlier data point of rate of ethylene production in Dacotah at stage V2 and hence was omitted from further statistics. The rate of ethylene production, which indicates the nitrogenase enzyme activity or BNF rates, was highest at stage V2 than those at stages V0 and E3 (Fig. 1). At stages V2 and E3, Alamo had significantly greater nitrogenase activity than Dacotah (Mann-Whitney U-test Pvalue = 0.017 at stage V2 and t-test with equal variance P-value = 0.030 at stage E3). Both Alamo and Dacotah showed (>22 lg ethylene per g system per h) reduction in BNF rates at E3 as compared to those at V2. We quantified the nifH gene in the rhizosphere samples of switchgrass cultivars (Fig. S1) using Sinorhizobium meliloti 1021 as the standard (C T = 22.064-4.0825 log(N 0 ); R 2 = 0.985; E std = 1.758), however, found no significant differences between cultivars.
The cultivars showed differences in the free soluble amino acid profiles Plant growth stages had a significant effect on the relative abundance of free soluble amino acids in the rootzone of the cultivars (adonis P-value = 0.001). Hence, further analysis was performed for each growth stage comparing the soluble amino acid profiles of Alamo and Dacotah. Adonis indicated that at stages V0 (P-value = 0.8) and E3 (P-value = 0.108), there were no differences in the relative content of root-zone soluble amino acids of cultivars. However, at stage V2 (Fig. S2), the relative abundance of soluble amino acids in the root-zones between Alamo and Dacotah was significantly different (adonis P-value = 0.04). The Mann-Whitney U-test showed that the cultivars had significantly (P-value <0.05) different levels of the following amino acids in their root-zones (Fig. 2): aspartic acid (Asp), glutamic acid (Glu), glutamine (Gln), glycine (Gly), and serine + asparagine (Ser + Asn) at V2; and arginine (Arg), histidine (His), and ornithine (Orn) at E3.

Sequencing data
Removal of barcodes, adapters, and primers using cutadapt and stitching with PANDAseq gave~1.9 million and~1.7 million high-quality 16S rRNA (bacterial) and ITS (fungal) gene sequence reads, respectively, for the 24 samples. The bacterial and fungal sequences, respectively, had average lengths of 333.2 and 254.3 bases with standard deviations of 57.7 and 48.7 bases. The raw data has been submitted to NCBI Sequence Read Archive (BioProject PRJNA320123) according to MIMS standard.

Alpha diversity
Using samples from all three growth stages and a sequence threshold (of 21 500 for bacteria and 16 000 for fungi), we compared the alpha diversity of groups (plant growth stage, cultivar, plant growth stage 9 cultivar) with a two-tail nonparametric t-test with FDR correction using the collated alpha diversity file. The bacterial alpha diversity at stage E3 was greater than that at stage V2 (PD whole tree q-value = 0.003, chao1 qvalue = 0.003). Furthermore, results from the interaction of plant growth stage and cultivar showed that the bacterial alpha diversity in Dacotah at stage E3 was greater than that at stage V2 (observed species q-value = 0.045, chao1 q-value = 0.045). The patterns of fungal alpha diversity resembled that of bacteria. The fungal alpha diversity at stage V2 was lesser than those at stages V0 (observed species qvalue = 0.0135) and E3 (chao1 q-value = 0.003, observed species q-value = 0.003). Also, results from the interaction of plant growth stage and cultivar showed that the fungal alpha diversity in Dacotah at stage E3 was greater than that at stage V2 (chao1 q-value = 0.03, observed species q-value = 0. 045). Certain patterns, however, were only present in the results from the interaction of plant growth stage and cultivar in the fungal data. The alpha diversity in Alamo at stage E3 was greater than those at stage V2 (chao1 q-value = 0.045, observed species qvalue = 0.0225) and stage V0 (chao1 q-value = 0.02).

Root-zone microbial communities structure (beta diversity) different between stage and cultivars
Multivariate data analyses in QIIME using MRPP, ADONIS, and ANOSIM on the weighted UniFrac distance between samples, using all OTUs, detected significant differences (P-value <0.05) in the bacterial communities as per stage and interaction of cultivar and stage (cultivar 9 stage), but not as per cultivar. These results were consistent when the analysis was performed on samples from (i) all three growth stages (V0, V2, E3), (ii) vegetative growth stages (V0 and V2), and (iii) vegetative (V2) and elongation (E3) stages (Table S3a). PCoA (Fig. S3a) and NMDS (Fig. S3b) analyses of bacterial communities using weighted UniFrac and Bray-Curtis distances, respectively, showed that samples from stages V2 and E3 clustered as per cultivars, however, samples from V0 did not. Hence, following analyses to compare cultivars were only performed at stages V2 and E3. Proteobacteria, Actinobacteria, and Acidobacteria were the most abundant phyla across samples (Fig. 3). Kruskal-Wallis test with Benjamini-Hochberg correction was performed on samples from stages V2 and E3 to identify significantly different (q-value <0.05) phyla. For each growth stage, the significant phyla, except 'Unassigned', were tested for differences (Pvalue <0.05) in relative abundance between cultivars using a bootstrap Mann-Whitney U-test in QIIME. For easier illustration, the 'Unassigned' phylum and phyla with total <1% abundance in samples across all growth stages were grouped as 'Other' taxa and phyla showing differential relative abundance between cultivars from the above test are shown by stars (Fig. 3). Stage V2 had more phyla with differential relative abundance between cultivars. Proteobacteria was significantly different only at stage V2; however, it consistently showed higher relative abundance in Alamo across all plant growth stages.
Similar to the bacterial data, multivariate data analyses in QIIME using MRPP, ADONIS, and ANOSIM on the Bray-Curtis distance between samples, using all OTUs, detected significant differences (P-value <0.05) in the fungal communities as per stage and interaction of cultivar and stage (cultivar 9 stage), but not as per cultivar (Table S3b). PCoA (Fig. S4) analyses of fungal communities using Bray-Curtis distances showed that samples from stage V2 clustered as per cultivars; however, samples from V0 and E3 did not. Sequences derived from the phylum Ascomycota were the most abundant across all growth stages (Fig. S5).

Root-zone bacterial communities showed functional differences between cultivars
Indicator species analysis (ISA) was performed for the bacterial and fungal data to identify statistically significant indicator taxa (Table S4). Functional profiles were predicted for bacteria using PICRUSt for the cultivars at the V2 and E3 stages before and after ISA. The predicted KEGG pathways in the metagenomes were different between cultivars at the corresponding growth stages for both before (Fig. S6) and after (Fig. S7) ISA. Before ISA, the KEGG pathway 'N-metabolism' showed significant differential abundance between cultivars at stage V2 (P-value = 0.039) but not at stage E3 (Fig. 4). However, after ISA, the pathway was observed to have significant differential abundance between cultivars at  . Two-sided t-test (P-value <0.05) was used to estimate differential abundance between cultivars.
stages V2 (p-value = 0.027) and E3 (P-value = 0.005) (Fig. 5). These indicate that the N-metabolism potential of the cultivar-specific root-zone-associated bacteria is different between Alamo and Dacotah.

Cultivars differed in the root-zone bacterial network correlation analysis
For the bacterial data, V2 and E3 samples showed differences as per cultivar and hence were used to generate bacterial COMET networks. In the bacterial COMET networks (Fig. 6), the number of classes (nodes) in Alamo and Dacotah was equal, whereas the total number of bacterial interactions in Alamo (72) and Dacotah (129) was different, with 24 interactions present in both networks (colored in black). Alamo had a higher number of modules (6) and percentage of positively correlated edges (56%) compared with the 3 and 46%, respectively, in Dacotah. For the fungal data, only samples from the V2 stage showed differences as per cultivar. As MENAP software needs at least eight samples per group, fungal COMET networks could not be generated.

Cultivable diazotrophs
After DNA extraction, 16S rRNA amplification, and Sanger sequencing of representative colonies cultivated on YM media, we obtained 45 and 42 high-quality reads for root-zone bacteria from Alamo and Dacotah, respectively. At stage V2, as per Libshuff analysis (Dacotah-Alamo P-value = 0.0054) and Pearson's chi-squared test (P-value <2.2e À16 ) (Fig. S8), we observed significant differences in the cultivable nitrogen-fixing bacteria (diazotrophs) between Alamo and Dacotah.

Discussion
The results of this study support the working quid pro quo model of plant-microbial interactions and the hypothesis that a higher productivity cultivar of switchgrass would be associated with a different root-zone bacterial community with greater nitrogenase activity than a lower productivity cultivar. The potential for greater carbon flow in exchange for bacterial fixed nitrogen provides a consistent framework to describe plant-bacterial interactions associated with switchgrass.    Fig. 6 The COMET networks of bacterial sequence abundance correlations from root-zones of (a) Alamo and (b) Dacotah. The nodes represent classes and they are colored as per the phyla. The red and green colored edges represent negative and positive interactions, respectively. A black colored edge represents that the interaction is present in both the Alamo and Dacotah networks.
Differences observed in the composition of amino acids, structure, cultivable diazotrophs, correlation networks, and the general functional capacity (e.g., nitrogen metabolism) of the bacterial communities associated with the two cultivars further support the above hypothesis. The results reinforce the idea that plant-diazotroph interactions and feedbacks (communications) are important descriptors of the greater nitrogen-fixing potential of the high productivity Alamo, compared with the Dacotah ecotype. Root-zone diazotroph communities, therefore, may help to meet the nitrogen demands of feedstocks and other grasses to support greater plant biomass production. The experimental outcomes support the idea that molecular breeding of switchgrass can be used to manipulate plant-bacterial interaction.

Plant growth-related N demand as a driver of switchgrass-diazotrophic interaction
Biological nitrogen fixation rates changed with growth stage and may reflect the temporal nitrogen demands of cultivars (Dong et al., 2001;Scagel et al., 2007). The increased BNF rates during the vegetative stages might be to meet the N demands of increasing biomass. On the other hand, the reduction in BNF rates during the elongation stage might be due to the nearing of maturation and the reduction in plant nitrogen demand, especially with its ability to translocate nitrogen to its roots for conservation during winter and reuse during regrowth (Monti, 2012). While it cannot be unequivocally concluded that the higher activities of nitrogenfixing bacteria are a result of Alamo having greater productivity than Dacotah, the results are consistent with this idea and track expected growth stage changes in plant nitrogen demand and carbon flow to below ground roots (Qian et al., 1997;Voisin et al., 2003).

Plant-microbial amino acids in the root-zone
Plant root and microbial exudates, especially amino acids, are used as signals of communication between plants and the root-zone microbiome (Morgan et al., 2005;Badri et al., 2009). The amino acids glycine, glutamate, alanine, aspartate, glutamine, and asparagine were found to be the dominant soluble amino acids in relatively mature switchgrass root-zone soil habitats (Figs 2 and S2); however, the contributions of these amino acids were altered with plant growth stage. The flux rates of these same molecules as exudates from plant roots have previously been shown to be relatively high (Lesuffleur et al., 2007;Carvalhais et al., 2011Carvalhais et al., , 2013 but can vary due to fertilization (e.g., N, Fe) and plant maturation. The relative concentration of glycine, in particular, was much higher in the root-zone of mature (E3) relative to young switchgrass plants (V0 and V2). There may thus be important changes in amino acid profiles that arise during plant maturation. The observations of amino acids in the root-zone of switchgrass represent a more complex habitat than studies that focus on axenic, microbialfree root exudation. It is thus difficult to come to firm conclusions regarding the biological role that amino acid dynamics might play in the root-zone. The variation across cultivar and growth stage, nevertheless, is consistent with the idea that root-zone amino acids reflect the contributions of two different but interactive plant-microbial systems.
Asparagine and glutamine tend to be relatively low in the rhizosphere, but can vary between <1 and~6% of the amino acid pool (Paynel et al., 2001;Lesuffleur et al., 2007). Intracellular concentrations within the roots of numerous plants (Paynel et al., 2001;Lodwig et al., 2003;Lesuffleur et al., 2007) of glutamine and asparagine tend to be relatively high, especially within phloem and xylem. Because of the importance of nitrogen-rich asparagine and the plant N shuttle glutamine, there is likely to be conservation and influx of these amino acids into organism biomass rather than efflux, particularly under N-limiting growth conditions. Hence, the relatively low concentrations of these two amino acids in the root-zone of switchgrass, relative to their constituent acids, aspartic and glutamic acids, may be the result of N-limiting growth conditions.
An amino acid feedback model has been proposed from the extensive study of legume (family Fabaceae) nodules (Dong et al., 2001). However, the application of this model to plant roots and associative diazotrophic bacteria in the rhizosphere remains to be determined. The rhizosphere (associative BNF system) is a different habitat than that of a nodule (symbiotic BNF system); however, application of this model provides a useful starting point for understanding interactions and exchanges between plants and bacteria in the root-zone. Glutamic and aspartic acids and the associated bases of glutamine and asparagine appear to provide a feedback signal that helps to control nitrogen fixation in nodules of legumes (Lodwig et al., 2003). These amino acids may thus be an important 'amino-stat' that helps to regulate the exchange of plant-bacterial carbon and nitrogen.
We used this nitrogen quid pro quo model to describe the dominant and dynamic amino acid species in the root-zone of the switchgrass ecotypes. We observed that aspartic acid was~threefold greater in Alamo than Dacotah, and conversely, glutamic acid was~1.4-fold greater in Dacotah than in Alamo. These differences could reflect different regulatory interactions between the plant and bacteria, perhaps related to differences in plant N demand and the availability of reduced nitrogen from diazotrophs (Lodwig et al., 2003). Assuming that the different amino acid concentrations reflect differences in the flows of these molecules between plant and bacterial biomass, the high levels of aspartic acid in the root-zone of Alamo could be interpreted as the result of higher rates of nitrogen fixation by diazotrophs. The greater flow of this bacterial-supplied N transporter could help meet the higher plant nitrogen demand of Alamo (Lodwig et al., 2003). In contrast, the relatively greater glutamic acid concentrations in the root-zone of Dacotah may reflect a buildup of this amino acid as a consequence of relatively low rates of nitrogen fixation. The results are also concordant with the observation that glutamic acid is secreted at lower levels by nitrogen-fixing bacteria under diazotrophic compared with adiazotrophic conditions (Gonzalez-Lopez et al., 1995). Reductions in glutamine and asparagine, and elevated aspartate-to-asparagine ratios, on the other hand, have been shown to be a response to low N availability that is translated as deficiency or stress within the plant (Amarante & Sodek, 2006) and thus help to describe plant N demand (Pate et al., 1980(Pate et al., , 1981(Pate et al., , 1984Atkins et al., 1983;Lea et al., 2007). Although we did not measure concentrations of amino acids within plant tissues, differences in these ratios in the root-zone may help to provide information explaining the greater nitrogen fixation potential associated with Alamo relative to Dacotah.
Changes in the structure, functions, and interactions of rhizosphere microbiome The broad differences in cultivable ( Fig. S8 and Libshuff results) and noncultivable bacterial (Figs 3, 6 and S3) and fungal (Figs S4 and S5) communities (Table S3) in the root-zone of Alamo and Dacotah are consistent with the above described functional indicators in nitrogen fixation potential (Fig. 1), root-zone exudates (e.g., amino acids, Fig. 2), KEGG analysis (Figs 4, 5, S6, and S7), and plant-microbial interactions (B€ urgmann et al., 2005;Mao et al., 2014). Different plant cultivars and growth stages have previously been shown to impact microbial community structure in plant root-zones (Berg & Smalla, 2009;Mao et al., 2011;Chaudhary et al., 2012;Chaparro et al., 2014). The research reported in this manuscript build upon these differences to show specific changes in microbial communities, potentially mediated through plant-microbial signals (amino acids), and concomitant functional changes in the root-zone.
Bacterial nitrogen fixation potential tracked changes in the structure and metagenomic functional potential of the root-zone microbial community. The increased N-metabolism potential in Alamo compared with Dacotah may reflect the overall differences in the nitrogen status of the ecotypes. The phylogenetic dominance of Proteobacteria, Actinobacteria, Acidobacteria, and Firmicutes observed in our results from culture-dependent and culture-independent methods across samples (Fig. 3) has been observed previously in roots of switchgrass (Jesus et al., 2010(Jesus et al., , 2016Mao et al., 2011;Plecha et al., 2013;Bahulikar et al., 2014). At stage V2, these phyla are among those that show differential abundances between ecotypes (Fig. 3), suggesting the potential importance of phylogenetic changes for driving functional level (Figs 4 and 5) shifts, such as that observed for the relative increase in N-metabolism potential in Alamo relative to Dacotah. The relative abundances of Proteobacteria and Actinobacteria, for example, are indicators of Alamo and Dacotah (Table S4), respectively. The observations from our study (Table S4, Figs 4 and 5) are consistent with the differential role that microbial communities and Ncycling bacteria may play in the growth of switchgrass ecotypes (Mao et al., 2011(Mao et al., , 2013. The ecological niches of the two switchgrass cultivars also provide a framework for understanding the potential for different types of root-zone bacterial interactions in relation to the greater nitrogen fixation rates. The differences in bacterial correlations of COMET networks between cultivars point out significantly different rootzone community dynamics among the ecotypes (Fig. 6). In addition, COMET networks summarize community differences using a different analytical framework than the multivariate pattern and statistical analysis, but nevertheless complement those findings (Tables S3 and S4, Figs S3, S4, S6, and S7). The common edges between Alamo and Dacotah COMET networks can be interpreted as the core set of bacterial associations (possible interactions) that are associated with switchgrass, and with the unique edges being cultivar specific. The lower number of modules but greater proportion of positive relative to negative correlations, if interpreted as associations, in Alamo relative to Dacotah, hints at the potential of bacteria-bacteria and plant-bacteria interactions that could serve to help the plant achieve higher productivity. Although it cannot be known whether these correlations reflect different interactions in the diazotroph root-zone communities of the ecotypes, the network model is useful for identifying differences between communities and the potential positive (or negative) feedbacks between microorganisms.
Positive and negative feedbacks between plants and microbes are often driven by environmental and nutrient limitations and thus provide a working hypothesis to explain how two ecotype root-zone habitats with different energy and nutrient conditions may develop different microbial communities and plant-microbial interactions (Morgan et al., 2005;Badri et al., 2009;Berg, 2009). The Proteobacteria, for example, have diazotrophic bacterial members and could help to explain their greater relative abundance and greater nitrogen fixation potential in Alamo relative to Dacotah. Similarly, the positive association of bacilli with many other taxa within the bacterial network, and its greater abundance in the cultivable community in Alamo compared with Dacotah is consistent with a growth supportive role played by these bacteria (Lopez-Bucio et al., 2007;Nautiyal et al., 2013;Xia et al., 2013;Gagne-Bourque et al., 2015). Presumably, the surrounding taxa that benefit would in return provide resources that support the growth of the bacilli. While the diversity of potential interactions in root-zone and soil habitats is staggering, and description of them is still in scientific infancy, the networks described herein are useful for understanding, predicting, and confirming possible positive or negative interactions among organisms.
In conclusion, although there is a scientific agreement regarding the importance of microbial interactions in the root-zone, the research presented is relatively unique in that it simultaneously studies the structure, potential for interaction, and functions of microbial communities in the root-zone of switchgrass. The research highlights possible linkages in plant nitrogen cycling (KEGG, amino acids) and demand associated with changes in the structure (COMET networks, PCoA, NMDS) and diazotroph activity (nitrogenase) of rootzone microbial communities. Understanding the interplay between the diazotroph communities with associated plants offers opportunities for genetic engineering and breeding of feedstock grasses and grain crops to select varieties that can better manage microbial communities to support growth and satisfy nitrogen demands of plants.

Supporting Information
Additional Supporting Information may be found online in the supporting information tab for this article: Figure S1. Boxplot of nifH copy numbers in cultivars at stages V2 and E3. Figure S2. Non-metric multi-dimensional scaling plots of the free soluble amino acid profiles in the root-zones of the cultivars at stages (a) V2 and (b) E3. Figure S3. (a) PCoA and (b) Non-metric multi-dimensional scaling showing differences in bacterial communities as per plant growth stage and cultivars. Figure S4. PCoA showing differences in fungal communities as per plant growth stage and cultivars. Figure S5. Phylum level taxonomic summaries of fungal communities in cultivars at the different growth stages. 'Unassigned' and less abundant taxa were grouped in 'Other'. Taxa in bars and legend from bottom to top are sorted as per decreasing abundance across all samples. Figure S6. Principal component analysis of the KEGG pathways predicted in the metagenomes of the cultivars before ISA at stages V2 and E3. Figure S7. Principal component analysis of the KEGG pathways predicted in the metagenomes of the cultivars after ISA at stages (a) V2 and (b) E3. Figure S8. CFU abundance of cultivable root-zone diazotrophs from Alamo and Dacotah at Stage V2. Table S1. Details about replicates used in this study. Table S2. Alpha diversity and counts per sample of (a) bacterial and (b) fungal data. Table S3. Multivariate data analyses of samples using all (a) bacterial and (b) fungal OTUs to test for significant effects (P-value <0.05) of stage, cultivar, and their interactions. Table S4. (a) Indicator bacterial taxa for Alamo and Dacotah at stages V2 and E3. Alamo-V2: 21 taxa; Dacotah-V2: 11 taxa; Alamo-E3: 14 taxa; Dacotah-E3: 38 taxa. (b) Indicator fungal taxa for Alamo at stage V2. Alamo-V2: 4 taxa; No significant taxa were found for Dacotah-V2.