The study was conducted under furrow irrigation at the Texas A&M AgriLife Research and Extension Center at Vernon. The soil is classified as a Miles fine sandy loam (Fine-loamy, mixed, super active, thermic Typic Paleustalfs). Peanuts (Span17) were planted on 14 May 2021 at 16 seed m⁻¹ using a four-row vacuum planter (1-m row spacing). The preceding year’s crop was cotton. Seed was treated with a commercial inoculum (Exceed^® Peat for peanut/cowpea/lespedeza/mung bean; Visjon Biologics) prior to planting. Plots were arranged in a randomized complete block design and replicated four times. Treatments included no cover crop and various cover crop treatments. For this study, peanut plants were collected from the no cover crop control plots only.

Nodules of two independent plant roots were collected from randomized plots in three developmental stages (R2, R4 and R7; Prasad, 1999). Whole peanut plants with roots and nodules were transported overnight on ice from the field site to the laboratory. Roots and nodules were washed thoroughly with tap water and wiped with a paper towel before collecting two types of nodules (big and small) from roots for surface sterilization. Ten big nodules (~ 1.5–2.0 mm), five from each plant, and 20 small nodules (<1.5 mm), 10 from each plant, were selected and immersed in 70% (v/v) ethanol for 2 min, washed with sterile H₂O three times, and then soaked in 10% (v/v) commercial bleach for 5 min. After bleach treatment, nodules were carefully washed 5–6 times with sterile H₂O and rinsed of any trace elements of bleach. A total of 12 big and 12 small nodule samples (3 growth stages x 4 plots) were independently collected for total DNA isolation.

Total DNA was extracted from the nodules using DNeasy Plant Pro Kit (Qiagen), according to the manufacturer’s instructions. DNA quality and concentration were determined using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Waltham, United States).

To confirm DNA amplifiability, the 16S rRNA gene was amplified using primers Eub338 and Eub518 (Fierer et al., 2005). To check for the presence of symbiotic marker genes (nifH and nodC) in the isolated nodules DNA, primers for nifH (Poly et al., 2001) and nodC genes (Sarita et al., 2005) were used for PCR amplification. The same PCR conditions were used for all three amplicons (16S, nifH, and nodC) which is as follows: 1 cycle of 95°C for 5 min (preheating), 30 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 1 min, and a final extension period at 72°C for 5 min. The size of the PCR product was confirmed by electrophoresis on a 2% agarose gel. For PCR reactions, nuclease-free PCR grade water, 2x GoTaq Master mix, and polymerase were purchased from Promega.

DNA sequencing was performed at the Texas A&M University Institute for Genome Sciences and Society (TIGSS) using Illumina MiSeq. Amplicon of V3-V4 regions of 16S rRNA gene amplified from 10 ng of DNA and 250 bp paired-end library cleaning, indexing, library quantification, normalization, pooling, and loading for MiSeq sequencing was carried out following the standard library preparation guide (Illumina, United States).

Nodule microbiome sequencing results were analyzed with the QIIME 2 software package version 2021.8¹ (Bolyen et al., 2019) and plugins associated with this version. Raw paired-end sequencing reads were demultiplexed using the demux plugin² with the demux emp-paired command. Quality control, filtering chimeric sequences, and feature table construction were done using the q2-dada2 plugin (Callahan et al., 2016) with trimming parameters (--p-trim-left-f 10 --p-trim-left-r 10 --p-trunc-len-f 250 --p-trunc-len-r 250) based on the demux visualization (QIIME 2 View). Taxonomy classification alignment with 99% similarity was done against the SILVA 138 database (Quast et al., 2012; Yilmaz et al., 2014) using a pretrained naive Bayes classifier and the “feature-classifier” plugin (Bokulich et al., 2018) with the “classify-sklearn.”

Based on “feature table” visualization, an additional filtering step was used to filter out samples with a total feature count less than 3,281. Furthermore, “low abundance features” (i.e., ASV’s with low frequency) were filtered from “feature-table” using “feature-table filter-features” command with “-p-min-frequency 10.” In addition, taxonomy-based filtering was conducted to remove mitochondria, chloroplast, Archaea, and Eukaryota using “taxa-filter-table and taxa-filter-seqs” commands from both taxonomic and rep-seqs tables.

After analyzing raw sequence data, we found more than 5 million sequence reads for both forward and reverse fragments from the 24 nodule microbiome samples with an average count of 111,646 reads per sample generated (Supplementary Table 1). The demultiplexed sequence data resulted an average of 59,560 non-chimeric sequence reads (52.66%; Supplementary Table 2) and inferred a total of 950 amplified sequence variants (ASVs) with an average number of features 40 and mean frequency 59,812 per sample. Since one of the samples (05-PNMB-104R4; Supplementary Table 2) generated a smaller number of non-chimeric sequence reads (6.71%), we further filtered out ASVs based on sampling depth to those that were less than 3,281 feature count/samples, and the minimum feature frequency was set to 10, leaving behind 492 ASVs from 23 samples (lowest feature count was retained to 35,401), which was analyzed further for community diversity and taxonomic classification.

A taxonomic bar-plot was generated using the filtered table and assigned taxonomy file and visualized through the QIIME 2 view website. The resulting filtered sequences from above were used to construct a phylogenetic tree using “align-to-tree-mafft-fasttree” pipeline from the q2-phylogeny plugin.³ Subsequently, phylogenetic diversity metrics were constructed with the resulting output “rooted phylogenetic tree.” Alpha and beta diversities were calculated through the q2-diversity plugin⁴ using the “core-metrics-phylogenetic” method with a sampling depth (rarefaction) of 1,290. Associations between categorical metadata columns (two types of nodule samples) and alpha diversity data were calculated through “qiime diversity alpha-group-significance” with Kruskal-Wallis pairwise test (Kruskal and Wallis, 1952) for Shannon’s diversity index. In addition, we also used PERMANOVA (default) to evaluate if the distances between the two nodule groups (sample compositions) were similar using the “qiime diversity beta-group-significance” plugin (Anderson, 2001) for weighted UniFrac distance. Further, we also conducted Mann–Whitney/Wilcoxon significance tests for alpha and beta diversity in two nodule groups and Kruskal-Wallis pairwise test for beta diversity analysis in three developmental growth stages of small nodules using the “rstatix” package in R. To find similarity within samples or dissimilarity between sample groups, a principal component analysis (PCA) was conducted using the centered log ratio (CLR) transformed microbial compositions. The PERMANOVA (n = 999) was done on Aitchison distance matrices to test significant effect of microbial compositions between two nodule types using the vegan package in R. A compositional heatmap of CLR-transformed microbial abundance (prevalence = 5) across samples was made with R package “microViz.”

To identify differentially abundant core microbial communities from two types of nodule samples, we used a statistical power analysis called “ANCOM” (analysis of composition of microbiomes) using ANCOM plugin (Mandal et al., 2015) with “qiime composition ancom” command. For ANCOM analysis, the ASV table that filtered out for low features, frequency, and for chloroplast, mitochondria, Archaea, and Eukaryota was used. Initially, taxa were collapsed using “qiime taxa collapse” command with - -p-level 6 to create a genus level table and later, add-pseudocount 1 for log0 (if any), and finally, “qiime composition ancom” command was used to tabulate significant (W-score) groups of bacteria at genus level that were differentially abundant (Estaki et al., 2020).

To predict functional profiles of pathways from 16S rRNA gene sequences of big and small nodule microbiome data (i.e., ASVs and abundance), we analyzed functional metabolic pathways for KO, EC, TIGRFAM, COG, CAZymes, and Pfam using a recently published metagenome analysis tool called “MicFunPred” (Mongad et al., 2021).

All the bioinformatics and biostatistics data were analyzed using Qiime-2 (Bolyen et al., 2019), Phyloseq (McMurdie and Holmes, 2013), and R packages in R Studio platforms (RStudio Team, 2020), and all the figures were generated using R version 4. 1. 2 (R Core Team, 2013).

Raw sequencing data were submitted to the NCBI Sequence Read Archive (SRA) under Bio Project accession PRJNA865795 and will be available upon publication. Data files and scripts that were used to generate figures are available at github.com/shak71/PNM.

Despite being inoculated with host-specific rhizobia, the field-grown peanuts formed both regular (referred to as big) sized (Figure 1A) and small nodules (Figure 1B; Please see M&M section for nodule size). These small nodules were present at all observed growth stages throughout the plant life cycle. While the larger nodules had a red interior color suggesting active N₂ fixation, the small nodules showed variable color (white, light pink, and gray or green) suggesting limited N₂ fixation activity (data not shown). However, N₂ fixation genes were found in both types of nodules. The nifH gene was found in all big and small nodules, except one small nodule while the nodC gene was not amplified from few small nodule samples (Supplementary Table 3).

To assess and investigate the bacterial community richness and diversity in two nodule types, the alpha diversity measures for observed ASVs and Shannon Diversity Index were calculated in big and small nodules of peanut roots. The number of observed ASVs differed significantly (p = <0.0001) between the two nodule types (Figure 2A). Observed ASVs in big nodule samples ranged from 5 to 10 while in small nodule samples ranged from 11 to 235. A Mann–Whitney/Wilcoxon significance test analysis for Shannon Diversity index (Figure 2B) was used to find the community richness which indicated a significant difference and variation between the big and small nodule microbiomes (p = 0.00029) with the bacterial communities in the small nodules being more diverse than those in the big nodules.

To analyze the bacterial community compositions in two types of nodules as well as growth stages for small nodules, the beta diversity index was calculated based on Weighted UniFrac distances. A Mann–Whitney/Wilcoxon significance test analysis was used among all samples grouped as big and small for Weighted UniFrac distance for the significance test and the results clearly showed that bacterial compositions differ between big and small nodule types (p = < 0.0001; Figure 2C). Since bacterial compositions differ within the small nodule category, we further explored the possibility within three developmental growth stages of small nodules and indeed, we found that bacterial compositions differ significantly (Kruskal-Wallis test, p = <0.0001) among all three growth stages of small nodules (Figure 2D). Further compositional aspects of centered log ratio (CLR) transformed PCA were analyzed for the bacterial communities in the big and small nodules as well as three different growth stages (R2, R4, and R7) of the small nodules. Based on the PERMANOVA test, the PCA results indicated that there were significant differences in bacterial communities between big and small nodules (p = 0.001); however, the bacterial communities among three growth stages of small nodules were also slightly different (p = 0.061; Figure 3A). Within small nodules, bacterial communities tended to group separately in these three growth stages except one or two data points in R4 and R7 stages that were closer to other growth stages. In contrast, bacterial communities in the big nodules seemed to cluster together and vary less across growth stages. Notably, the small nodule communities were initially (R2) somewhat similar to the big nodules but became more distinct at later growth stages (R4 & R7). At the genus level, based on top taxa loading on the PCA plot as well as direct visualization of sample compositions on the iris plot (Figures 3A,B), big nodule samples contained mostly Bradyrhizobium, while other taxa were often more abundant in small nodule samples. Collectively, the data clearly indicated that the small nodules had more diverse and temporally variable bacterial communities than the big nodules.

All the filtered ASV’s from the big and small nodules were classified as bacteria and included 15 phyla, 25 classes, 57 orders, 87 families, and 151 genera (Figure 4; Supplementary Table 4). Both nodule types had large relative populations of Proteobacteria (Figure 4A), with the phylum representing ~99.96% of the bacterial community in the big nodules and ~ 68.67% in the small nodules combinedly from three growth stages. Thus, the big nodules harbored a tiny fraction of native community microbes while the small nodules were inhabited by a large portion (~31.32%) of diverse microbes. Further analysis revealed that only Actinobacteriota (~0.03%) co-inhabited with Proteobacteria in the big nodules; however, 15 different phyla including Actinobacteriota were detected in the small nodules (Figure 4B).

Interestingly, the microbiome data showed that these diverse bacterial populations in the small nodules dominated over time, becoming more abundant at later growth stages – the most noticeable being members of the Actinobacteriota, Bacteroidota, Chloroflexi, and Patescibacteria (Figure 4A). Surprisingly, the phyla Patescibacteria and Chloroflexi were not enriched at all in the early R2 growth stage. In addition, we found that Proteobacteria abundance and occupancy in the big nodules was steady throughout the growth stages. However, Proteobacteria were comparatively lower while the diverse microbial populations appeared higher in the small nodules irrespective of growth stages. At the phylum level, the microbial abundance and distribution frequency was analyzed since the most dominant and abundant phyla were Proteobacteria and Actinobacteriota compared to other taxa. The data revealed that microbial abundance and frequency varied from low to high abundance in the big and small nodules for these two phyla with more distinct bacterial populations present in small nodules. To further explore the variation of microbiota within these three growth stages of the small nodules, 118 taxa at the genus level were extracted based on relative abundance without any uncultured taxa (Figure 4C; Supplementary Table 5). The Venn diagram showed the taxa (relative abundance >0) distribution in R2, R4 and R7 growth stages of small nodules, and indicated there was diverse microbial dominance over the growth cycle of the small nodules. Interestingly, we also found some unique as well as overlapped microbiota between these three growth stages with the greatest overlap between growth stages R4 and R7.

To identify the dominant and abundant taxa in the big and small nodules, we first removed all the unclassified species-level data, the rare and uncultured taxa, the taxa with a minimum prevalence of 5, and then a CLR-transformed compositional heatmap was made at the genus level (Figure 5A). The small nodule samples from the three growth stages had a diverse bacterial community while the big nodule samples were dominated by a single genus Bradyrhizobium spp. The core bacterial genera in the big and small nodules were identified using ANCOM differential test analysis at the genus level (Figure 5B). Based on ‘W’ score > 100, 7 genera were identified: Bradyrhizobium, Streptomyces, Niastella, Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium, Sphingomonas, Novosphingobium, and Chitinophaga. These seven core genera were also found in the compositional heatmap, suggesting that these are the dominant and highly abundant core genera in big and small nodules. Of these seven abundant and dominant core genera, based on W score, only Bradyrhizobium was highly abundant in the big nodules, while at least 6 genera including Streptomyces, Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium, Niastella, Novosphingobium, Sphingomonas, and Chitinophaga were primarily abundant in the small nodules.

Predicted functional profiles of pathway genes (based on pathway abundance) of bacterial communities from big and small nodules of peanut roots were generated using ‘MicFunPred’ (Mongad et al., 2021). Utilizing 16S rRNA gene sequences and normalized ASV’s abundance from three developmental stages of big and small nodules, metagenomes in terms of KEGG Orthology (KO) functional categories (based on log10 transformed pathways abundance) were used to investigate the enriched function for pathway genes related to nitrogen (N) and phosphorus (P) cycling processes (Figures 6, 7; Supplementary Table 6). Based on search criteria as “nitrogen metabolism” at C-level hierarchy (Supplementary Table 6), a total of 32 genes from the big nodules and 36 genes from the small nodules related to the N cycling process were enriched. A heatmap was made based on all the 32 genes that were common in both big and small nodule developmental stages (Figure 6). The relative abundance of most biological N₂ fixation (nifH, nifD, and nifK), assimilatory (nirA, nasA, and nasB) and dissimilatory (nirB, napA, and napB) nitrate reduction, and denitrification (norB and norC)-related genes were higher in all three growth stages of big nodules compared to small nodules at all growth stages (Figure 6). However, the predicted pathway data showed that the relative abundance of some genes or different subunits related to the N cycling process was higher in at least two growth stages (R4 and R7) of the small nodules. For example, nirD—nitrite reductase (NADH) small subunit, nitrilase, gdhA—glutamate dehydrogenase (NADP+), nasB—assimilatory nitrate reductase electron transfer subunit, GLUD1_2, gdhA—glutamate dehydrogenase (NAD(P)+), arc—carbamate kinase, nitrate reductase / nitrite oxidoreductase alpha subunit (narG, narZ, nxrA), nitrate reductase / nitrite oxidoreductase beta subunit (narH, narY, and nxrB), and nitrate reductase gamma subunit (narI, narV; Figure 6). Interestingly, a few genes related to N cycling processes that were uniquely enriched only in the small nodules were observed. These unique nitrogen-related pathway genes were Glutamate synthase (ferredoxin; GLU, gltS), ferredoxin-nitrate reductase (narB), nitrous-oxide reductase (nosZ), and hydroxylamine reductase (hcp; Supplementary Table 6).

Like the N cycling processes, KO pathway genes related to the P cycling processes were analyzed (Figure 7). Based on search criteria “phosphate metabolism/pathway” at C-level hierarchy (Supplementary Table 6), a total of 48 genes from three growth stages of big nodules and 58 genes from small nodules related to P cycling processes were enriched. A heatmap was made based on all the 48 genes that were common in both big and small nodule developmental stages (Figure 7). Among 48 phosphate-related genes, almost 50% (26) genes showed higher relative abundance in big nodules growth stages compared to small nodules. However, the relative abundance of 22 genes related to P cycling processes was also higher in all three stages of small nodules category, particularly 3-phytase and 4-phytase/acid phosphatase, myo-inositol 2-dehydrogenase (iolG), fructose-1,6-bisphosphatase II (glpX), deoxyribose-phosphate aldolase (deoC), ribose 5-phosphate isomerase B (rpiB), myo-inositol-1-phosphate synthase (INO1), 6-phosphofructokinase 1 (pfkA), transaldolase (talA, talB), and glucose-6-phosphate isomerase (GPI; Figure 7). The unique phosphate-related genes enriched in small nodules were inositol oxygenase (MIOX), 1-phosphatidylinositol phosphodiesterase (plc), 2-keto-myo-inositol isomerase (ioll), diphosphate-dependent phosphofructokinase (pfp), 6-phospho-3-hexuloisomerase (hxlB), 3-hexulose-6-phosphate synthase (hxlA), dehydrogluconokinase (kguK), glucose-6-phosphate isomerase-archaeal (pgi1), and CDP-diacylglycerol--inositol 3-phosphatidyltransferase (CDIPT), and fructose-1,6-bisphosphatase III (fbp3). Besides N & P cycling genes, other pathway genes were analyzed that were unique only in the small nodule growth stages. All the duplicates for KO identifications from both big and small nodules were removed, and then sorted out for unique pathway genes in small nodule category only, of which 383 genes were observed that were unique and enriched (log10 abundance >4.0 in all three stages; Supplementary Table 7). Among the observed 383 pathway genes, a heatmap was made for 35 genes that were common among three stages of pre-selected top 50 pathway genes (Figure 8). All 35 pathway genes were highly enriched in growth stage R7 compared to R2 and R4, suggesting that these genes might play some role in the late stage of the small nodules. The most noticeable and highly enriched genes were IMP dehydrogenase (IMPDH, guaB), NADH-quinone oxidoreductase subunit J, K, N (nuoJ, nuoK, and nuoN), and succinate dehydrogenase / fumarate reductase (sdhC, frdC), and phosphate regulon sensor histidine kinase PhoR (phoR). It is quite interesting to observe that some key metabolic categories other than N and P were specifically increased in the small nodules (Supplementary Table 7) and this warrants future studies for further characterization.