16S rRNA Gene Amplicon Based Metagenomic Signatures of Rhizobiome Community in Rice Field During Various Growth Stages.

Rice is a major staple food across the globe. Its growth and productivity is highly dependent on the rhizobiome where crosstalk takes place between plant and the microbial community. Such interactions lead to selective enrichment of plant beneficial microbes which ultimately defines the crop health and productivity. In this study, rhizobiome modulation is documented throughout the development of rice plant. Based on 16S rRNA gene affiliation at genus level, abundance, and diversity of plant growth promoting bacteria increased during the growth stages. The observed α diversity and rhizobiome complexity increased significantly (p < 0.05) during plantation. PCoA indicates that different geographical locations shared similar rhizobiome diversity but exerted differential enrichment (p < 0.001). Diversity of enriched genera represented a sigmoid curve and subsequently declined after harvest. A major proportion of dominant enriched genera (p < 0.05, abundance > 0.1%), based on 16S rRNA gene, were plant growth promoting bacteria that produces siderophore, indole-3-acetic acid, aminocyclopropane-1-carboxylic acid, and antimicrobials. Hydrogenotrophic methanogens dominated throughout cultivation. Type I methanotrophs (n = 12) had higher diversity than type II methanotrophs (n = 6). However, the later had significantly higher abundance (p = 0.003). Strong enrichment pattern was also observed in type I methanotrophs being enriched during water logged stages. Ammonia oxidizing Archaea were several folds more abundant than ammonia oxidizing bacteria. K-strategists Nitrosospira and Nitrospira dominated ammonia and nitrite oxidizing bacteria, respectively. The study clarifies the modulation of rhizobiome according to the rice developmental stages, thereby opening up the possibilities of bio-fertilizer treatment based on each cultivation stages.

Introduction

Soil comprises mineral particles of different shapes and sizes with distinct chemical characteristics and organic compounds in various stages of decomposition. Molecular studies suggest that soil harbor the highest biodiversity of organisms on Earth, approximately 1,000 Gbp of microbial genome sequences per gram of soil (Vogel et al., 2009). These microorganisms play important roles in maintenance of soil fertility, nutrient cycling, and carbon sequestration. The rhizosphere is the soil around living roots influenced by root activity (Hartmann et al., 2008); the chemical composition of the exudates, border cells and mucilage released by plant roots are important to signal and also establish the microbial community in this environment (Mendes et al., 2013). The rhizobiome can influence the plant directly and/or indirectly (Schnitzer et al., 2011), as seen by the PGPB (plant growth promoting bacteria), which acts synergistically on plant growth promotion and disease suppression (Liu et al., 2018). The microbial diversity present in rhizosphere of native and cultivated crops can be different due to the species-specific effects (Philippot et al., 2013). For example, the bulk microbiome is more complex than rhizobiome of Dendranthema grandiflora Tzvelev (chrysanthemum) (Duineveld et al., 2001), sagebrush (Gonzalez-Franco et al., 2009), and oak (Uroz et al., 2010). On the other hand, the opposite trend was observed for other plants, such as wheat (Donn et al., 2015), wild oat (Shi et al., 2016), and switch grass (He et al., 2017). Additionally, the root exudate composition changes during plant growth and influence microbial community assembly in the rhizosphere, as shown in Arabidopsis (Chaparro et al., 2014), soybean (Sugiyama et al., 2014), grapevine (Novello et al., 2017), and maize (Walters et al., 2018) rhizosphere. The microbial community assembly in rhizosphere is also determined by the abiotic and biotic factors influencing both natural and agricultural ecosystems (Philippot et al., 2013). Studies performed in greenhouse showed differences in the rice rhizobiome cultivated in soils with different fertilizer (Innes et al., 2004) and in rice-wheat rotation cultivation in the field (Wang et al., 2016).

Rice is the major stable food for more than half of the world’s population, and more than half of the supply comes from India and China (Muthayya et al., 2014). Due to the high consumption of this crop, rice paddies occupy around 11% of the arable land in the world. However, rice cultivation is responsible for 20% of total agricultural CH4 emission (Van Groenigen et al., 2011). High greenhouse-gas (GHG) emission occurs by the growth of methanogenic archaea present in the vicinity of rice roots (Neue, 1993). On the other hand, the increase of CO2 and temperature will elevate the CH4 emission, indicating a positive feedback loop. This suggests that higher photosynthetic rate results in an increase of substrates availability for the methanogenesis (Tokida et al., 2010). The release of CH4 (methane) into the atmosphere is further aided by the rice vascular system, which serves as a passage for a part of CH4 from soil to the atmosphere, bypassing the methanotrophic community (Hernández et al., 2015). The anaerobic conditions in paddy rice also preclude the nitrification, decreasing the NO3– production, and promoting the complete denitrification of NO3– to N2, resulting in low N2O emission (Verhoeven et al., 2018).

Although scientific understanding about the diversity and composition of soil microbiome can directly and/or indirectly influence a wide range of ecosystem processes, few efforts have been directed toward understanding the microbial biogeography (Fierer and Jackson, 2006). To address this question, we have undertaken a spatial and temporal characterization of rhizobiome and bulk soil from rice paddy using targeted amplicon-based (16S rRNA gene) metagenomic approach, and we posit the following predictions: the microbial community (i) is significantly different in bulk and rhizosphere soils, (ii) changes across the rice cultivation stages, and (iii) is distinct according to the geographical scale analyzed.

Materials and Methods

Sampling and Environmental DNA Isolation

The rhizosphere and bulk soil samples were collected during the month of December 2016 to April 2017 in seven different India states (Table 1 and Supplementary Figure S1). All of them are situated in rural areas away from any industrial set up in the vicinity. All the sites practiced rice monoculture for more than 10 years and during off-season are kept as fallow land, and cow manure was used after plowing. Ammonium, phosphate, and sulfate (Potash, FACTAMFOS) were used during seedling stage. The rice rhizosphere soil was sampled in triplicates by pulling out the plant to collect the soil loosely attached to the roots from three different geographical sites in all growth stages: vegetative, reproductive and ripening. Moreover, bulk soil was also collected in pre-plow, post-plow, and after harvest stages. To investigate the relationships among the bacterial community similarity and geographical distance, additional samples were taken from vegetative, and post-plow stages in paddy fields in different sites in India (Table 1). All the samples were collected from upper 5 – 15 cm in 5 replicate subsamples, which were maintained in individualized sterile plastic bags. The environmental DNA was extracted using HiPurATM Soil DNA purification kit (HiMedia) following the manufacturer’s protocol. In order to increase the DNA yield, the bead-beating step was extended to 15 min, and incubation during elution step was extended to 30 min at 37°C. Soil pH was measured, and ammonia and nitrite were analyzed using 2M potassium chloride extract as described elsewhere (Klotz and Stein, 2011).

TABLE 1

Description of study sites.

Geographic coordinates	Sampling stage (samples code)
Vadakara, Kerala(11°33′42.4″N 75°39′28.3″E)	Pre-plow (kvp1)Post-plow (kvp2)Vegetative (kvv)Reproductive (kvr1)Ripening (kvr2)Post-harvest (kvph)
Chemmattamvayal, Kerala(12°18′48.0″N 75°06′30.2″E)	Pre-plow (kcp1)Post-plow (kcp2)Vegetative (kcv)Reproductive (kcr1)Ripening (kcr2)Post-harvest (kcph)
Kheezmala, Kerala(12°16′50.1″N 75°12′42.7″E)	Pre-plow (kkp1)Post-plow (kkp2)Vegetative (kkv)Reproductive (kkr1)Ripening (kkr2)Post-harvest (kkph)
Dimapur, Nagaland(25°54′52.7″N 93°47′18.6″E)	Vegetative (nl)
Rasipuran, Tamil Nadu(11°26′40.5″N 78°10′26.6″E)	Vegetative (tn)
Ukhurl, Manipur(25°05′45.6″N 94°21′19.3″E)	Vegetative (mn)
Burari, Delhi(28°46′26.9″N 77°12′11.9″E)	Post-plow (dl)
Panvel, Maharashtra(19°01′01.4″N 73°07′32.1″E)	Post-plow (mh)
Mysore, Karnataka(12°30′39.4″N 76°54′05.7″E)	Post-plow (my)

16S rRNA Gene Amplification and Next Generation Sequencing

The prokaryotic hypervariable V3–V4 region from 16S rRNA gene was amplified using the primers set Pro341F (5′ – CCTACGGGNBGCASCAG – 3′) and Pro805R (5′- GACTACNVGGGTATCTAATCC- 3′) (Takahashi et al., 2014). Equimolar quantities of PCR amplicon in replicates were pooled and sequenced using Illumina MiSeq platform at Eurofins Scientific (Bangalore, India). The raw fastQ files were uploaded to the metagenome rapid annotation using subsystem technology (MG-RAST) server (Meyer et al., 2008) and annotated using default parameters. Briefly, the pipeline removes low quality sequences having phred scores below 15 (Cox et al., 2010). Artificial duplicate reads were removed using DRISEE (duplicate read inferred sequencing error estimation) (Gomez-Alvarez et al., 2009). The sequences were annotated using the RDP (Ribosomal Database Project) database having a minimum cutoff identity of 60% and e-value of 5. All samples are publicly available at MG-RAST under the IDs provided in Supplementary Table S1.

Statistical Methods

The alpha diversity was estimated by observed richness, Chao1, and Shannon diversity. The diversity was estimated using Shannon (H′) index [H′ = −Σni/n ln (ni/n), where ni is the number of individuals in the taxon i and n is the total number of individuals], which is influenced by both species richness and evenness. Difference in alpha diversity in various growth stages was calculated using Kruskal Wallis. Rarefaction curves were plotted to estimate whether the sequencing depth is enough to wholly capture the microbial diversity. Statistical difference in bacterial composition based on plantation, cultivation stages and geographical locations was calculated using PERMANOVA (permutational multivariate analysis of variance using distance matrices) and Bray Curtis distance. The p values were based on 999 permutations, and the adonis function of vegan package was used. Taxonomical network was constructed in R (R Core Team, 2014) with a jaccard maximum distance of 0.05 and 0.5 through make_network function from phyloseq v1.16.2 package. Network visualization was performed in Gephi-0.9.2 (Bastian et al., 2009). The samples were clustered in two groups, cropped (vegetative, reproductive, and ripening stages) and uncropped (pre-plow, post-plow, and post-harvest steps).

Results

Soil analysis for each growth stages showed a consistent pH of 4 – 6. Similarly, growth stages also had a little effect on nitrite concentration ranging from 0.8 to 1.2 μM having higher fluctuation between samples during the growth stages (Supplementary Figure S2). However, ammonia concentration was quiet distinguishable between uncropped (15 – 30 μM) and cropped (25 – 45 μM) stages showing almost 2× fold increase once the plant reached reproductive stage. High-throughput sequencing of the 16S rRNA gene amplicons generated more than 6.53 Gb, representing a total of 7.1 million reads. Over 207.825.194 bp of the total reads passed quality control (Supplementary Table S1). The high-quality reads were clustered using >97% sequence identity into 3856 microbial OTUs. The rarefaction curves reached a stable asymptote for all groups (Supplementary Figure S3). A significant portion of the microbiome fell under uncultured/unclassified taxa (Supplementary Table S2).

Soil Microbial Community Differs in Rhizosphere Across Rice Cultivation

Twenty-eight phyla were obtained from all groups (Supplementary Figure S4) and the most abundant were Proteobacteria (25.69 ± 13.92%), Firmicutes (20.82 ± 8.69%), Actinobacteria (16.68 ± 5.93%), and Acidobacteria (13.28 ± 10.47%). The rhizosphere had a higher number of genera (n = 844) than the bulk soil (n = 735), and both of them shared 79.36% of these genera (n = 700) (Figure 1A). Twenty-three genera showed more than 1% of occurrence from any sample (Supplementary Figure S5). The most abundant genus from rhizosphere group was a Candidatus Koribacter (7.5–18.6%) while the most frequently found genus in bulk soil samples was Ktedonbacter (13.4%). Measures of within sample diversity (alpha diversity) showed a gradient from bulk soil to rhizosphere samples; however, the difference in alpha-diversity could not be considered statistically significant (Table 2). PERMANOVA analysis showed that the community composition of the bulk soil was statistically different (Table 3) and higher complex network structure than rhizosphere (Figures 1C–E). Pairwise comparison of community dissimilarity variances showed that significant differences between post-plow to ripening and post-plow to post-harvest groups (Table 4). The community structure was also statistically different considering geographic location (Table 3), which was corroborated by weighted unifrac PCoA (principal coordinates analysis) (Figure 2B). On the other hand, unweighted PCoA analysis showed the separation of rhizosphere and bulk soil samples irrespective of geography, indicating similarities in the microbial composition instead of abundance (Figure 2A).

FIGURE 1

(A) Venn diagram of cropped and uncropped microbiome shows a high similarity of abundance. (B) Genera with more than 10 folds enrichment in rhizobiome. (C,D) Number of nodes and edges for 0.5 and 0.05-jaccard distance, respectively. (E) 0.05 jaccard distance network complexity map.

TABLE 2

Pairwise comparisons of alpha diversity between each stage.

	Pre-plow	Post-plow	Vegetative	Reproductive	Ripening	Uncropped
Post-plow	1.00
Vegetative	0.62	0.90
Reproductive	0.90	0.33	1.00
Ripening	0.90	0.33	0.90	1.00
Post harvest	1.00	1.00	1.00	1.00	1.00
Cropped						0.005

TABLE 3

Permutational multivariate analysis of variance using distance matrices (PERMANOVA) of bulk and rhizosphere microbiome compared to cropping, growth stages and location.

Group comparisons	Bray		Manhattan		Canberra		Jaccard
	R2	P value	R2	P value	R2	P value	R2	P value
Cropped vs. Uncropped	0.102	0.0002	0.102	0.0005	0.074	0.0004	0.083	0.0015
Cultivation stages	0.162	0.0193	0.162	0.0173	0.181	0.0174	0.162	0.0304
Geographical location	0.505	0.0003	0.505	0.0003	0.384	0.0006	0.478	0.0001

TABLE 4

Significant differences in microbial community among the groups using ANOSIM.

	Pre-plow	Post-plow	Vegetative	Reproductive	Ripening	Post-harvest
Pre-plow		0.22	0.08	0.26	0.33	–0.22
Post-plow	0.14		0.02	0.27	0.40	0.5
Vegetative	0.32	0.36		–0.22	–0.07	0.06
Reproductive	0.20	0.14	0.91		0.15	–0.04
Ripening	0.33	0.046	0.57	0.5		0.19
Post-harvest	1	0.027	0.36	0.6	0.1

FIGURE 2

PCoA rhizobiome analysis of different cultivation stages and location using (A) unweight unifrac shows clustering of cropped samples irrespective of location and cultivation stage. In contrast, (B) weight unifrac PCoA shows culturing majorly based on location implying that the diversity is relatively similar, but not abundance, among the samples during cropped stages.

Growth Stage-Dependent Variations in Microbial Community

All genera (n = 34) that were significantly (p < 0.05) altered in pre-plow stage were suppressed (Figure 3D). Similarly, several genera were suppressed in post-plow stage, except for Kibdelosporangium, Heliophilum, and Oceanimonas (Figure 3E). During vegetative stage, number of enriched genera increased further (n = 11) while the percentage of the depleted genera reduced (39%) (Figure 3F). The richness and abundance of enriched genera increased while the percentage of the down-regulated genera reduced with plant growth (Figures 3C–I). Interestingly, majority of the genera that experienced significant changes during plant growth were up-regulated while those genera with significant changes before plantation or after harvest were depleted. In order to catalog rhizobiome growth stage-dependent enriched genera in each location, all genera having above 1% abundance in any cultivation stage were screened to obtain the dominant genera and a heatmap was constructed (Supplementary Figure S6). In line to the weight PCoA, rhizobiome enrichment had more difference than similarities among the location. Several genera such as Ca. Koribacter, Ca. Solibacter, Clostridium, Pseudonocardia, Geoalkalibacter, Saccharopolyspora, and Acidobacterium were enriched commonly during the plantation in all locations. However, among those common genera, there was no consistent enrichment patterns with respect to growth stage between the locations. For instance, Ca. Koribacter was enriched during vegetative and ripening stage at site KK, while the same genus was enriched during reproductive stage and post-harvest stage for site KC and KV, respectively. Similarly, Ca. Solibacter was enriched during all the cultivation stages at site KK and KV while it was enriched majorly at reproductive stage in site KC. The highly abundant genus, Ktedonobacter, was depleted gradually as the cultivation progressed at site KK and KC, which was re-enriched only after post-harvest. However, it was enriched at site KV during post-harvest, as well as vegetative stage.

FIGURE 3

Growth stage-dependent selective enrichment (A) at genus level of the top 10 most abundant genera from each cultivation stage. (B) Enrichment and depletion on a larger time frame (cropped vs. uncropped). (C–I) Growth stage-dependent genus enrichment and depletion (p < 0.05). Blue color denotes enrichment, red color denotes depletion.

Methane Metabolism and Nitrification

Eighteen methanogen genera were found in all samples, which accounted for 0.96 – 1.35% of overall abundance in each stage. The abundance rose sharply post-plow probably due to the application of organic cow manure, and declined with plant age until post-harvest. Dominant genera were Methanosaeta, followed by Methanobacterium and Methanocella (Figure 4A and Supplementary Figure S7). Archaeal genera belonging to type I and type II methanotrophs were present throughout the cultivation. Abundance of Methanosarcina was constant until vegetative period, decreased at reproductive (p < 0.046), and increased at ripening stage (p < 0.046). Methermicoccus abundance increased at ripening stage and decreased after harvest. Predominant methanotroph genera belonging to first group were Methylococcus and Methylocaldum, while Methylocystis and Methylosinus were predominant in the second group (Figure 4B and Supplementary Figure S8). There was a significant (p = 0.003) difference between type I and type II methanotroph abundance. Methanotrophs belonging to the Verrucomicrobial phylum, previously considered as extremophiles, were also enriched during the cropped and post-harvest stages. Ammonia concentration increased significantly at reproductive stage, including ripening and post-harvest (26 – 48 μM) in contract to the initial growth stages (19 – 28 μM). However, no significant differences were observed between bulk soil and rhizosphere groups based on either ammonia or on nitrite concentration (Supplementary Figure S9). AOA (ammonia-oxidizing Archaea) dominated over AOB (ammonia-oxidizing bacteria) throughout the cultivation ranging from 3 to 30 folds (p = 0.003) (Figure 4C). Genera Candidatus Nitrosocaldus and Nitrosospira were the dominant members of AOA and AOB, respectively. Within the nitrite oxidizers, Nitrospira was the most dominant followed by Nitrobacter and Nitrococcus (Figure 4D). Type I methanotrophs, nitrite oxidizers, and AOA showed a similar pattern during the ripening stage having reduced abundance while AOB increased.

FIGURE 4

Rhizobiome composition and diversity of (A) methanogens, (B) methanotrophs, (C) ammonia oxidizers, and (D) nitrite oxidizers in a stage wise manner shows differential enrichment within respective niche throughout the cultivation.

Discussion

Alpha Diversity Increases Gradually in Cropped Rhizosphere During Plantation

Higher alpha diversity have been reported in rhizosphere, rather than bulk soil, of rice (Edwards et al., 2015), maize (Yang et al., 2017), and cotton (Qiao et al., 2017). In contrast, studies have also highlighted higher alpha diversity in bulk soil, rather than rhizosphere, in maize (Peiffer et al., 2013), soy plant (Wang P. et al., 2017), and annual grass Avena fatua (Shi et al., 2015). Cultivation practice such as rotation crop system of alfalfa-rice, wheat-rice (Lopes et al., 2014) and continuous monocropping of black pepper (Xiong et al., 2015) exhibited higher alpha diversity in uncropped soil. Vegetative stage of rice rhizosphere during flooded condition had higher alpha diversity compared to reproductive (matured) stage (Pittol et al., 2018). Such contrasting literature shows a missing gap between rhizosphere and bulk soil alpha diversity. In this study, gradual increase of alpha diversity was observed with plant age irrespective of location or genotype (Figures 5A–C). During the reproductive and ripening stage, the alpha diversity had the lowest standard deviation indicating a stable diversity. Alpha diversity between any two consecutive stages didn’t show significant difference (P > 0.05) (Table 2). However, comparison of alpha diversity between rhizosphere and uncropped soil showed statistically significant difference in observed alpha diversity (P < 0.05). This indicates that the enrichment in rhizosphere alpha diversity is a gradual process with no significant changes from one cropping stage to another, however, the changes are prominent when viewed at a larger timeframe i.e., cropped vs. uncropped.

FIGURE 5

Growth stage-dependent (A) richness, (B) Shannon, and (C) evenness represented in box plot.

Rhizosphere Microbial Community Is Diverse and Complex Yet Specialized

Denaturing gradient gel electrophoresis (DGGE) studies on Dendranthema grandiflora Tzvelev (chrysanthemum) have indicated rhizosphere as a subset of bulk soil (Duineveld et al., 2001). These observations are further backed up by NGS based metagenomics studies on soybean (Mendes et al., 2014) having lesser network complexity in the rhizosphere relative to bulk soil. However, rhizosphere of wheat and wild oat had microbial diversity increased with plant age (Donn et al., 2015) and markedly more complex than the surrounding bulk soil although the diversity was decreased (Shi et al., 2016). Cultivation of switch-grass also showed an increased microbial network structure in rhizosphere (He et al., 2017). In this study, community co-occurrence network was highly dependent on the threshold maximum distance between the nodes. The number of nodes (vertices), as well as the edges (links) of the co-occurrence network increased with plant growth with a strong correlation (r = 0.9, p < 0.03) between nodes and edges when the maximum jaccard distance between the nodes were set to 0.5. Although the number of nodes remained similar pre and post-plow (Figures 1C–E), the edges increased ∼1.5 folds post-plow indicating the increased association among the genera. During the cultivation period, the complexity (nodes and edges) increased in every succeeding stage, which continued until after-harvesting stage corroborating to the high alpha diversity. High correlation was found between the number of nodes and edges against alpha diversity (r = 0.58 – 0.95; p = 0.003 – 0.066). In contrast, network map having maximum distance of 0.05 had number of nodes as well as edges decreased significantly during the cropped stages (Figures 1D,E and Supplementary Figure S10). The contrasting observation in network map indicates that the microbial network in rhizosphere constitutes a higher number of complex relationships compared to the bulk soil; however, they share a lower number of strong linked relationships. This decline of strong linked relationships could be due to selective enrichment of beneficial microbes exerted by the root exudates which leads to partial depletion of non-beneficial microbes. This also indicates that the cropped and uncropped soil samples shared a common pool and the selective enrichment is highly regulated by the root exudates. In support of this, Venn diagram of the cropped and uncropped samples indicated that the number of genera unique to cropped stages increased; however, they constituted merely 0.626% of the overall abundance (Figure 1A). Meanwhile the common genera (n = 700) between cropped and uncropped stages accounted for over 99% of the abundance. A closer look at each genus showed high enrichment of selected genus (Figure 1B). This indicates that the rhizosphere microbiome, in addition to high diversity and complexity, is highly selective in rhizobiome recruitment.

Geography Influences Community Abundance

Physiological factors such as growth stages have shown relatively weak influence over determination of the rice rhizobiome and highlighted the importance of geography in shaping the rice rhizosphere microbial community (Edwards et al., 2015). However, in such studies, the geographical location distances were relatively small and hence a larger variation is required as a confirmatory study. Hence, metagenome from after-plow and vegetative stages were selected from each location in Kerala and compared to metagenomes of similar growth stages from locations separating thousands of miles apart in order to have a large variation in geography. India, being a sub-continent, has tremendous variation in geography, the climatic condition across such distance exhibits a considerable variation. Weighted UniFrac PCoA analysis showed clustering of samples based on geography (Figure 2B). However, an unweighted UniFrac PCoA plot clearly clustered the cropped stage irrespective of location (Figure 2A). The disparity between the PCoA plots could suggest that the microbial diversity in rice rhizosphere is influenced by the plant and remains relatively similar in spite of different location. However, their abundance is highly influenced by environmental factors due to the differential enrichment based on the environment. Similarly, unweighted UniFrac distances of maize rhizobiome had strong similarity deepened with plant age (Walters et al., 2018). In agreement to the PCoA, ADONIS analysis revealed location (R2 = 0.468 ± 0.049; p < 0.05) (Table 3) as the major determinant of the microbial community followed by plantation and growth stages. The strong influence of geography in shaping the community structure is in line with previous studies which state that the soil type and soil development stages had a major role compared to the plant effect (Mapelli et al., 2018). Further visualization of the dominant genera diversity and abundance (>1% abundance in any cultivation stage) in heatmap showed clear difference, albeit with some similarities, for each location, and cultivation stage (Supplementary Figure S6). Enriched dominant genera common to all location included Ca. Koribacter (nitrogen metabolism) (Júnior et al., 2019), Ca. Solibacter (nitrogen metabolism), Clostridium (nitrogen-fixing) (Doni et al., 2014), Pseudonocardia (siderophore production), Geoalkalibacter (syntrophic organic matter degradation) (Neveling et al., 2017), Saccharopolyspora (phosphate-mineralizing) (Franco-Correa and Chavarro-Anzola, 2016), and Acidobacterium (indole-3-acetic acid production) (Kielak et al., 2016b). Their prevalence could be because of the symbiotic association with the plants. For instance, in addition to IAA production, Acidobacterium spp. are able to reduce Fe3+ (ferric) to Fe2+ (ferrous) (Kielak et al., 2016b) which can be taken up by the plant.

Rhizobiome Enrichment Is Strongly Linked to Plant Development

Abundance of ten top most abundant genera from each cultivation stage decreased during the cultivated stages as compared to uncultivated stage (pre-plow and post-harvest) (Figure 3A). The combined dominance of the top 10 most abundant genera during cropped stages reduced providing room for higher microbial diversity and complexity. Ca. Koribacter (14.61 ± 5.57%), Ktedonobacter (10.02 ± 3.56%), Ca. Solibacter (5.75 ± 2.51%), and Peptoniphilus (4.57 ± 1.68%) were among top 4 genera which consistently dominated in all the stages. Sharp differences in enrichment and depletion of the genera were prominent when the dominance was compared from one growth stage to another. Ca. Koribacter and Ca. Solibacter were enriched during the cropped stages while Ktedonobacter and Peptoniphilus were depleted. Ca. ndidatus Koribacter and Ca. ndidatus Solibacter possess large genomes with several gene duplications and paralogs obtained through horizontal gene transfer (Challacombe et al., 2011) and adapted to nutrient limited conditions (Eichorst et al., 2011), as well as to acidic soil (Grza̧dziel and Gała̧zka, 2018). High abundance of Ca. Koribacter and Ca. Solibacter were also reported in several soil (Navarrete et al., 2013; Júnior et al., 2019). However, their ecological functions is poorly understood due to lack of pure representative culture (Kielak et al., 2016a). Nonetheless, genomic insights have indicated that Ca. Koribacter and Ca. Solibacter are oligotrophic K-strategist have a strong role in reduction of nitrite and nitrate which helps in Nitrogen cycle (Fierer et al., 2007; Ward et al., 2009). Ca. Solibacter have a versatile genome which is adapted to breakdown complex recalcitrant organic compounds and carbohydrates and provides a suitable environmental condition for other microbes (Pearce et al., 2012; Rawat et al., 2012; Rime et al., 2015; Li et al., 2018). Previous reports have also highlighted their abundance in rhizosphere of black pepper (Umadevi et al., 2018), Panax notoginseng (Tan et al., 2017) and rice rhizosphere (Lopes et al., 2014). In line with previous report, Ca. Solibacter was positively correlated to rice cropping. Ktedonobacteria belonging to the Chloroflexi phylum were also highly abundant. However, their abundance depleted during cropped stages. Ktedonobacter were linked to tobacco disease since Ktedonobacter are non-nitrogen fixers which may compete with the microbial community for nitrogen source (Niu et al., 2016). However, in the pinus rhizosphere, Ktedonobacter were postulated to play a critical role in maintaining the microbial community structure (Wang et al., 2018). The genus Peptoniphilus under Firmicutes phylum was also found to decrease during cropped stages. Peptoniphilus has not been linked to plant rhizosphere to the best of our knowledge; however, there are several clinical reports of Peptoniphilus on human intestinal occlusion (Cobo et al., 2017), bloodstream infections (Brown et al., 2014), and various cancers (Shilnikova and Dmitrieva, 2015).

In order to study the patterns of genera having significant (p < 0.05) changes during cultivation stages, the rhizobiome enrichment in a stage wise manner was analyzed by identifying the entire genera which differed significantly (p < 0.05). A trend was observed in both enrichment as well as depleted microbes in a stage wise fashion (Figures 3C–I). Gradual enrichment was observed in the initial stages of cultivation followed by exponential increase between the vegetative and reproductive stage. This was followed by the formation of a plateau, which declined sharply after harvesting. This is expected since root exudates of rice as well as other plants have been documented to change according to plant maturity (Aulakh et al., 2001). Hence, the recruitment of beneficial microbes in the rhizosphere would vary as the plant matures. In addition, secondary metabolites from the rhizobiome aids in crosstalk with the microbes which further stimulates to secrete favorable exudates (Lareen et al., 2016). Such mutualistic relationship is necessary to avail several essential elements for plant development. For instance, bioavailability of phosphate in soil is critical for plant growth due to its role in morphological and physiological functions such as energy storage and transfer, cell division, root elongation, and photosynthesis, etc. (Yulianti and Rakhmawati, 2017). However, available forms of phosphate even in fertile soil is generally low at sub-micromolar levels (Arif et al., 2017). Phosphate solubilizing bacteria plays a critical role in solubilization of various phosphate compounds making it available to the plants (Arif et al., 2017), counteract soil calcification (Adnan et al., 2017), and enrich microbial diversity (Wang J. et al., 2017). Plant hormones are one of the major determinants of rooting and shooting (Habib et al., 2016), enhances drought resistance (Jung et al., 2015), up-regulates nitrogen fixation (Defez et al., 2017), and as bio-fertilizers (Shahzad et al., 2017). Such reprograming and cross talk between rhizobiome and plant leads to selective enrichment of beneficial microbes. Consistent with this, majority of the enriched dominant genera (p < 0.05, relative abundance > 0.1%) during the growth stages were PGPB involved in sulfur and nitrogen cycle as well as production of metabolites that aids in plant growth such as siderophore, IAA, ACC, phosphate solubilization, and antimicrobials. Similarly, further analysis between all the cropped stages vs. uncropped stages yield 143 genera (p < 0.05) (Supplementary Table S3) out of which 25 genera had relative abundance higher than 0.1% (Figure 3B). All the enriched genera (n = 19) related to IAA, ACC, siderophore, phosphate solubilization, etc. were represented during cropped stage while the suppressed genera (n = 6) were dominant in uncropped stage. Similar observations were made on soybean rhizosphere having higher abundance of plant growth promoting rhizobacteria in a stage dependent manner (Sugiyama et al., 2014). Enrichment of the rhizobiome was pronounced during the cropped stages and weaker during the off cultivation stages (Figure 3C). A strong correlation (r = 0.99; p = 0.02) was also observed between the number of enriched and suppressed genera during the growth stages. Previous studies have hypothesized that the high diversity and abundance of PGPB in the rhizosphere could decrease the non-PGPB due to the intense competition for resource and interference in the microbiome network.

Aceticlastic Methanogens Dominate Rice Rhizosphere

Methane is a candidate for renewable energy source having high yield per mass unit (55.7 kJ/g) and burns cleaner than the fossil fuels (Richards et al., 2016). Hence, it holds a promising position in the energy industry. However, methane is 25 folds more potent greenhouse gas than CO2 (Yvon-Durocher et al., 2014). It is also the second most abundant greenhouse gas. Rice paddy contributes a significant amount of CH4 emission into the atmosphere (Smith et al., 2007) due to the ideal anoxic condition brought about by the prolonged waterlogged irrigation system. Methanogenesis is exclusively carried out by methanogens belonging to Archeal domain mainly by using H2/CO2 (cytochrome-lacking hydrogenotrophic methanogens) or acetate/methylated compounds as substrates (cytochrome-containing methylotrophic methanogens) (Vaksmaa et al., 2017). Methane emission in rice paddy is directly correlated to abundance of methanogen (Singh et al., 2017). Hence, in this study, microbial community of methanogens were investigated. Abundance of methanogens consistently declined until harvest (Figure 4A). Soil ammonium concentration analysis indicated an opposite trend to methanogens; ammonium concentration as well as AOA/AOB relative abundance increased with plant age (Figure 4C and Supplementary Figure S2). Nitrification produces NOx– compounds which are partly taken up by the roots and the rest diffuses into the surroundings having adverse effect on methanogens (Banihani et al., 2009). In addition, methanogens require anaerobic condition for methanogensis. Hence, they are generally found with increased depth of 30–200 mm (Lee et al., 2015). We observed that the rice paddy methanogens mainly consisted of the Methanosaeta, Methanocella, and, Methanobacterium throughout the cultivation. Similar abundance has also been found in the rice paddy from South Korea, Philippine, Italy, and China (Lee et al., 2014; Breidenbach and Conrad, 2015; Hernández et al., 2015; Vaksmaa et al., 2017; Wang et al., 2018; Yuan et al., 2018). Methanosaeta was the most dominant methanogen in this study. It is capable of acetoclastic pathway for CH4 production and exhibits high affinity to acetate even at low concentration (Rotaru et al., 2014). Previous study from Indian rice paddy during cultivation has also reported the dominance of methanogenesis through acetoclastic pathway (Singh et al., 2012; Bhattacharyya et al., 2016). The second and third most abundance methanogen Methanocella and Methanobacterium, respectively, belongs to hydrogenotrophic methanogens. Methanocella belongs to the order Methanocellales which were formerly regarded as rice cluster I (RC-I) (Lü and Lu, 2012). They were found to produce CH4 from rice photosynthates (Lu and Conrad, 2005) and have adaption in low H2 environments (Liu et al., 2014).

Type II Methanotrophs Are Less Diverse Yet More Abundant Than Type I

Methanotrophs resides in the rhizosphere owing to the sub-oxic zones created by the rice exudates. Methanotrophs are responsible for oxidizing at least half of the methane produced by methanogens before it is released to the atmosphere (Thauer, 2010). Hence they have been considered as a CH4 biological sink under aerobic as well as anaerobic condition (Leng et al., 2015). A strong pattern was observed in type I Methanotrophs being enriched with plant age until the drainage of logged water (Figure 4B). Similar to previous report from Indian rice paddy (Pandit et al., 2016b), a high diversity of type I methanotrophs were observed (Rahalkar et al., 2016; Pandit et al., 2018). Unlike type II Methanotrophs, type I methanotrophs have atmospheric nitrogen fixing capability and derive nitrogen from atmospheric N in the absence of other N sources (Noll et al., 2008). Hence, selective enrichment by root exudates could establish a symbiotic relationship due to the nitrogen fixing properties of type I Methanotrophs. However, genus Methylocystis (type II methanotrophs) was the most abundant among the methanotrophs (Figure 4B). Various other studies have observed similar dominance of Methylocystis in India (Pandit et al., 2016a), China (Liu et al., 2017), Taiwan (Shiau et al., 2018), and South Korea (Lee et al., 2014). This could be due to its survival physiology in dry environment during the off season by forming desiccation-resistant cysts (Pandit et al., 2016a). It can also be attributed to the dual nature of Methylocystis to obtain carbon from methane as well as acetate (Leng et al., 2015). On a broader view, similar to previous studies (Lee et al., 2014), type II-Methanotrophs didn’t show any significant enrichment with plant age. In vitro studies by have shown that type II Methanotrophs were inhibited in paddy soils by ammonium treatment possibly due to the substrate competition between methane monooxygenase and ammonia monooxygenase (Reay and Nedwell, 2004) and also by the toxicity of ammonia oxidation byproducts such as hydroxylamine and nitrite (Saari et al., 2004). However, the negative effect of ammonia could also be counteracted to some extent by the end product of ammonia oxidation, nitrate, which stimulates the growth of both type I, and type II Methanotrophs (Hu and Lu, 2015). Genus Methylacidiphilum under phylum Verrucomicrobia was found throughout the cultivation (Figure 4B). Methylacidiphilum have also been recently reported in Italian paddy soil (Vaksmaa et al., 2017). This genus is of special interest because they are the most acidophilic methanotrophs which was originally discovered from a geothermal ecosystem (Op den Camp et al., 2009). They can thrive at pH 2 (Islam et al., 2008) and found mostly in extreme conditions such as geothermal environment and volcanic soil (van Teeseling et al., 2014). Its presence in paddy soil hints toward the need to understand its role in the rice paddy methane oxidation (Hou et al., 2008). Overall, type I methanotrophs showed a higher diversity (n = 12) than type II Methanotrophs (n = 6) yet type II methanotrophs had a significantly higher abundance (p = 0.003).

AOA Abundance Dominates Over AOB in Flooded Rice Paddy

AOA was represented several folds higher than AOB in all the cultivation stages (Figure 4C). In line to our observation, meta-analysis of East Asian paddy soils have found similar dominance of AOA over AOB in acidic as well as alkaline paddy soils (Mukhtar et al., 2017). AOA have been shown to be highly adaptive in microaerophilic condition (Wang et al., 2015) and greatly influence by the rice root exudates (Chen et al., 2008). We found that within AOA community, Nitrosocaldus cluster dominated the paddy soil. Nitrosophaera cluster was the 2nd most dominant group in our study. However, Nitrosophaera cluster was the most dominated within AOA community in several studies from China (Chen et al., 2008; Ke et al., 2013; Wang et al., 2015; Lu et al., 2018) and East Asian paddy soils meta-analysis (Mukhtar et al., 2017). The disparity between our study and literature could be attributed to the major difference in PCR primer. Previous studies have implemented the widely used Archaeal amoA gene primers which has several mismatch base pairs to Nitrosocaldus cluster (Ke et al., 2013). Hence, comparative analysis of amoA based literature should be taken with precaution. As for AOB community, Nitrosospira was the most dominant. Similar observation were found in alkaline as well as neutral paddy soil from China (Chen et al., 2008; Lu et al., 2018), Japan (Bowatte et al., 2006), East Asian paddy soils meta-analysis (Mukhtar et al., 2017). AOA has been found to be dominating in acidic soils (pH < 5.5) (Liu et al., 2013). They are also highly adaptive to low ammonia concentration having highest affinity toward ammonia among all ammonia oxidizers and tolerate sub-oxic condition created in the rhizosphere through the oxygen leakage from root surface and diffusion of atmospheric oxygen (Wang et al., 2015). In addition, AOA has also been reported to thrive in wet environment (Hu et al., 2015). Corroborating to this, AOA abundance increased sharply after-plow which declined after drainage during pre-harvest stage indicating that the flooded water irrigation system favors AOA. Candidatus Nitrosocaldus genus (AOA) was the most dominant among all ammonia oxidizers. AOB abundance increased gradually with plant growth till pre-harvest stage (Figure 4C). Hence, AOB were less affected by the change in the flooded condition. AOB are known to thrive in high ammonium concentration (Jia and Conrad, 2009). Within AOB, Nitrosospira was the most dominant and follows a K strategy having high affinity to low ammonia concentration (Zhang et al., 2016). Overall the diversity of AOA (n = 3) as well as AOB (n = 4) was low. Previous study (Han et al., 2018) has noted that the diversity of ammonia oxidizers is minimum in long-term water input. Since the sampled geographical location experiences wet tropical climate, the low diversity of ammonia oxidizers is expected. Within the nitrite oxidizing bacteria, Nitrospira genus was found to be most abundant (Figure 4D). Nitrospira play a major role in terrestrial ecosystems and withstand high concentration of nitrite (Han et al., 2018). Nitrospira are fast-growing K-strategists having high affinity to substrate even at low nitrite and oxygen concentrations having resistance to ampicillin and Acriflavine. It has been postulated that the high efficiency of Candidatus Nitrospira defluvii could be due to its unique presence (found only in Nitrospira, Nitrospina, and Anammox bacteria) of nitrite oxidoreductase, a key enzyme in NO2– oxidation, in periplasmic region rather than cytoplasm. This orientation favors the generation of more proton-motive force per nitrate oxidized (Rodríguez et al., 2015).

Conclusion

16S rRNA gene amplicon based metagenomic analysis of rice bulk soil and rhizobiome was carried out at different cultivation stages showing sharp distinction in the alpha diversity and enrichment as the plant matured. The rich network complexity during the plantation points out the enormous diversity of microbes, which can be harnessed not only for improved plant health as bio-fertilizer but as a reservoir of novel enzymes and metabolites yet to be explored. Looking at an ecological point of view, it is counterintuitive to find that the methanotrophs declined with plant growth in the rhizobiome. However, methanogensis at the higher depth rhizosphere would be interesting to investigate further in Indian paddy rice irrigation system. In addition, the higher abundance of AOA over AOB further adds to the evidence that AOA dominates the ammonia oxidation in a majority of ecosystems.

Although the analysis of microbial community through 16S rRNA gene amplicon based metagenomic is a powerful tool and provides a significantly larger view as compared to DGGE or traditional culture based techniques. However, it comes with certain limitation. Firstly, such technique relies on the relative abundance rather than absolute abundance of the microbial community. Recent reports as per Naylor et al. (2017) and Fitzpatrick et al. (2018) have indicated that the relative abundance from such high throughput data could reflect the absolute abundance fluctuation. Secondly, the functional potential of the microbial community analyzed based on genera could have multiple outcomes based on the species and strains. Thirdly, 16S rRNA gene amplicon based metagenomic relies on the use of primers for amplification of the 16S rRNA gene, which could introduce biases. However, Ibarbalz et al. (2014) reported that despite the 16S rRNA gene amplification bias, it provides a strong quantitative and qualitative information of the microbial community yet the archeal community such as methanotrophs, methanogens, and ammonia oxidizers could be under represented (Raymann et al., 2017).

Data Availability

The raw fastQ files for the Kerala, India samples were uploaded in MG-RAST server and publicly available from the MG-RAST server under IDs mgm4758099.3 to mgm4758122.3.

Author Contributions

RK and MI conceived and designed the experiments. RK, MI, DB, and AV analyzed the data. RK contributed reagents, materials, and analysis tools. RK, MI, AV, AG-N, NK, OP, DB, PG, and VA wrote the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.