Seasonal dynamics alter taxonomical and functional microbial profiles in Pampa biome soils under natural grasslands

Study sites

The study was conducted during 2014 and 2015 in two sites under similar environmental conditions located in the Brazilian Pampa biome: Santa Maria municipality—SM (29°45′S, 53°45′W) and São Gabriel municipality—SG (30°20′S, 54°15′W), 83 km far from each other. The sites were chosen because of their similarities in soil features, vegetation, land use, and climate. Soil samples were collected under natural grassland, currently used for grazing of cattle, with no historic of fertilizers input (except for the manure added by animal activity). The sites exhibited similar soil physicochemical characteristics (Table S1) and were classified within the same soil taxonomy.

The sites chosen are exposed to large seasonal variations in soil moisture and temperature (Fig. 1). For this study, the seasons were characterized as warm season (October, November, December, January, February and March) and cold season (April, May, June, July, August and September), due to the similarities in temperature and rainfall between spring and summer as well as autumn and winter. The average air temperatures in the cold season were 16.2 °C ± 2.4 and 17.5 ± 2.7 °C in SM and SG, respectively. In the warm season the average air temperatures were 24.0 ± 2.8 °C in SM and 24.1 ± 1.4 °C in SG. The monthly average rainfall in the cold season was 107 mm in SM and 112 mm in SG. In the warm season, monthly average rainfall was 154 mm and 171 mm, in SM and SG, respectively (Fig. 1).

The plant species comprising 95% of the total biomass in SM were: Adropogon lateralis, Aristida laevis, Elephantopus mollis, Bacharis trimera, Paspalum plicatulum and Paspalum notatum, and in SG were: Andropogon lateralis, Axonopus affinis, Baccharis trimera, Eragrostis plana (invasive), Erianthus angustifolius, Paspalum notatum, Paspalum dilatatum, Paspalum plicatulum, Sporobolus indicus and Vernonia nudiflora.

Soil sampling, nucleic acids extraction and preprocessing

At each grassland (SM and SG), the sampling area was divided in three blocks, and five soil cores per block were randomly selected and pooled to make a single composite sample, resulting in three samples per site at each season (n = 24, that is, two sites (SM and SG), 4 months (May and August, 2014—cold season/November 2014 and February 2015—warm season) and three replicates per sampling). Samples were collected with a sterile spatula to a depth of 5 cm, kept in sterile 50 mL tubes in liquid nitrogen, and then stored at −80 °C until DNA and RNA co-extraction.

From each sample, 2 g of bulk soil (free of roots or rhizospheric soil) was used for simultaneous extraction of total RNA and DNA using the RNA PowerSoil kit and the PowerSoil^®DNA Elution Accessory Kit (MoBio laboratories, Inc., Carlsbad, CA, USA), following the manufacturer’s instructions. After RNA extraction, each sample was subjected to enzymatic digestion of DNA with TURBO DNA-free Kit (Applied Biosystems, Foster City, CA, USA). All RNA and DNA samples were submitted to quantification and quality-check using a Qubit RNA or DNA assay kit (Invitrogen, Eugene, OR, USA) and NanoVue™ spectrophotometer (GE healthcare, Chicago, CA USA), respectively, and further stored at −80 °C until library preparation.

Metatranscriptome library preparation and sequencing

One microgram of total RNA from each bulk soil sample were depleted by removing rRNA from total RNA with the MICROBExpress™ Bacterial mRNA Enrichment Kit (Thermo Fisher, Waltham, MA, USA) following the manufacturer’s instructions. The mRNA was further purified with the MEGAclear™ Transcription Clean-Up Kit (Thermo Fisher, Waltham, MA, USA) following the manufacturer’s instructions. The enriched mRNA was used to prepare the mRNA library with the Ion Total RNA-Seq Kit v2 and Ion Xpress™ RNA-Seq Barcode Kit (Thermo Fisher, Waltham, MA, USA). The libraries were amplified by emulsion PCR and sequenced on Ion 316™ Chip v2 using the Ion Torrent PGM system and the Ion PGM™ Sequencing 400 kit, according to the supplier’s instructions. After sequencing, all reads were filtered by the PGM software to remove low quality and polyclonal sequences.

16S library preparation and sequencing

Library preparation and sequencing followed the procedures described by Suleiman et al. (2016). For analyzing the composition of soil bacterial and archaeal communities the total microbial DNA was used as template for partial 16S rDNA amplification and sequencing. The V4 region of the 16S rRNA gene was amplified using the bacterial/archaeal primers 515F and 806R (Caporaso et al., 2012a), and sequenced using the PGM Ion Torrent (Thermo Fisher Scientific, Waltham, MA, USA). Multiple samples were PCR-amplified using barcoded primers. Each of the 25 µL of PCR mixture consisted of 2U of Platinum^® Taq DNA High Fidelity Polymerase (Invitrogen, Carlsbad, CA, USA), 4 µL 10× High Fidelity PCR Buffer, 2 mM MgSO4, 0.2 mM dNTP’s, 0.1 µM of both the 806R barcoded primer and the 515F primer, and approximately 100 ng of DNA template. The PCR conditions used were 95 °C for 5 min, 30 cycles of 94 °C per 45s denaturation; 56 °C per 45 s annealing and 72 °C per 1 min extension; followed by 72 °C per 10 min. The resulting PCR products were purified with the Agencourt^® AMPure^® XP Reagent (Beckman Coulter, Brea, CA, USA) and the final concentration of the PCR product was quantified by using the Qubit Fluorometer kit (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s recommendations. Finally, the reactions were combined in equimolar concentrations to create a mixture composed by 16S rRNA gene amplified fragments of each sample. This composite sample was used for library preparation with Ion OneTouch™ 2 System with the Ion PGM™ Template OT2 400 Kit Template (Thermo Fisher Scientific, Waltham, MA, USA). The sequencing was performed using Ion PGM™ Sequencing 400 on Ion PGM™ System using Ion 316™ Chip v2.

Bioinformatics analysis and statistics

The annotation of metatranscriptome sequences was performed with the Metagenomics Rapid Annotation (MG-RAST) pipeline version 3.3.6 (Meyer et al., 2008), using the standard parameters for sequence quality control. The data was compared to the SEED Subsystem using a maximum e-value of 1⁻⁵, a minimum identity of 60%, and a minimum alignment length of 15 measured in aa for protein and bp for RNA databases. The 16S rRNA raw sequences were analyzed following the recommendations of the Brazilian Microbiome Project (Pylro et al., 2014), using the BMP Operating System (BMPOS) (Pylro et al., 2016).

Analysis of Metagenomic Profiles v2 (STAMP) software package was used to determine differences in relative abundances of microbial functions and taxa (Parks et al., 2014). Statistical hypothesis tests were performed using the Welch’s t test while confidence intervals were calculated using the Welch’s inverted method and Bonferroni multiple test for p-value correction. The mRNA and the 16S rRNA datasets were rarefied to the same number of sequences per database (Lemos et al., 2011) and used to construct dissimilarity matrixes generated by Bray Curtis and Binary distances using the “phyloseq” package in R. The matrixes were ordinate by Principal Coordinate Analysis (PCoA) and adonis function was used to calculate the Permutational Multivariate Analysis of Variance (PERMANOVA) and verify the strength and statistical significance of groups among location, season and the combined effect of location and season with the vegan package (Oksanen et al., 2015). Sampling effort was estimated using Good’s coverage (Good, 1953).

Raw sequences obtained by metatranscriptome sequencing and associated metadata were submitted to MG-RAST server (http://metagenomics.anl.gov/) and are publicly available under the string mgp12046. All raw sequences obtained by amplicon sequencing were submitted to NCBI Sequence Read Archive (SRA) and are available under the experiment number SRX2779549 and run number SRR5499445.

Overall microbial differences among seasons and locations

After quality filtering the 16S rRNA reads, a total of 1,630,136 high-quality sequences longer than 200 bp were retained. The average Good’s coverage of 98% was calculated (Table S2), indicating the dataset was representative of the microbial communities analyzed. The soil metatranscriptome sequencing yielded a total of 20,584,533 reads for all 24 samples, and the average Good’s coverage was 77% (Table S2).

To evaluate the differences in microbial community structure and function between the two natural grasslands and along different seasons, the abundances of microbial taxa and mRNA encoding functions were used to compute a Bray–Curtis dissimilarity matrix coordinated by using PCoA. A Binary similarity matrix was also used to compute differences of either presence/absence of taxa or mRNA encoding functions. Contrasting patterns of microbial taxonomic structure and functions were observed. Although the sites were similar, these analyses revealed differences in microbial taxonomic structure between the two natural grasslands (SM and SG) and between seasons (Fig. 2). Four distinctive groups were observed taking into account the microbial taxonomic structure meaning the soil microbial phylotypes shifted along the year according to the season. When the abundance of mRNA encoding microbial functions was compared, two distinctive groups were observed, representing the cold and the warm seasons. Those results indicate that microbial functions are more resilient to environmental changes than the microbial taxonomic structure. No significant differences in microbial taxonomic structure or microbial functions were observed within the same season in either location. The location did not influence the main functions performed by the microbial communities (Fig. 2). Microbial community differences were further confirmed by PERMANOVA (Table 1). Season was the main driver of shifts of abundance of taxa and mRNA encoding functions of the soil microbiota. This factor contributed 22% and 23% for the variation in microbial taxa and mRNA encoding functions, respectively (Table 1). Altogether, these results indicated that soil microbial community structure was affected by location (e.g., two natural grasslands from the same biome, but 83 km apart) and by the natural environmental changes throughout the year (seasons). Also, location and season interaction was significant at 5% probability, meaning the impact of one factor depends on the level of the other factor or, in other words, that microbial community shifts emerge from the interaction location and seasons. On the other hand, the functions performed by soil microbes in those two sampling sites were similar, not being affected by location, but by seasons.

Identification of taxa and functions that shifted along seasons and locations

Once overall taxonomic and functional differences were identified, the functional differences as measured through RNAseq, and taxonomical differences as measured by 16S rRNA sequencing, were determined. The seasonal frequency of most abundant microbial phylotypes (relative frequency greater than 1%) in both locations was determined (Table 2). The data indicated the grassland ecosystem maintained a core microbiota along the year but the relative abundance of taxa varied in response to seasons. As the amount of rain was well distributed along the year (Fig. 1), soil moisture was not considered an important variable shaping microbial community assembly in our experiment. Temperature shifted during the year with an average minimum of 10 °C and an average maximum of 27 °C, representing an important source of environmental variation influencing microbial life in our experiment. During the warm season, the uncultured member of Verrucomicrobia designated DA101 was the most abundant phylotype in both locations. Mycobacterium and Rhodoplanes were also found in both grasslands with similarly high abundances, but better adapted to cold conditions (Table 2).

In agreement with the differential abundance analysis, alpha diversity measurements indicated similar microbial richness (number of observed OTUs) during the seasons, but greater microbial diversity (Simpson diversity index) during the cold season (Fig. 3). Cold temperatures allowed for greater evenness while warm temperatures decrease the diversity of taxa. The diversity of microbial functions presented similar trends (high diversity during the warm season and low diversity during the cold season), displaying a core set of microbial functions along the year in both locations (Fig. 3 and Table 3). Main functional differences between seasons were consistently detected within both grasslands. Despite the difficulties of defining the functions codified by complex mixtures of mRNA in metatranscriptome library, the metabolism of carbohydrates was the dominant set of functions performed during the cold season in both grasslands. An average of 10.5% of mRNA with known annotated function was committed to carbohydrate degradation and utilization in the model bacterium Thermotoga maritima, during the cold season in the grassland located in SG. Nevertheless, only 1.1% of the mRNA sequences were related with this function in the warm season (Table 3). The same pattern was observed in the grassland located in SM (7.1% of mRNA sequences committed to sugar utilization during the cold season and only 0.9% during the warm season). D-ribose utilization and deoxyribose and deoxynucleoside catabolism were also functions highly expressed during the cold season in both grasslands. During the warm season functions were more evenly distributed, without dominance of specific ones. The most abundant set of functions during the warm season were those dedicated to the protein metabolism, composing 5.7% of mRNA dataset in the grassland located in SG, and 3.5% in the SM grassland. The functions related to Cofactors, Vitamins, Prosthetic Groups and Pigments were also higher expressed during the warm season in both grasslands.