Soil was collected from drained rice fields at the Italian Rice Research Institute (IRRI) located in the valley of the Po river (Vercelli, Italy). Slurry microcosms were set up by filling 125 mL bottles with 40 g dry soil and 35 mL of autoclaved water. A total of 0.5 g rice straw (Oryza sativa var. Roma) pieces (1–2 cm) was added as the main carbon source (e.g., Conrad and Klose, 1999, 2006; Wegner and Liesack, 2016). The carbon–nitrogen (C/N) ratio of rice straw is in the range of 50:1 to 60:1, with a carbon content of 30–40% and a nitrogen content of 0.6–0.7% (Maramba, 1978; Wang et al., 2014). The thoroughly mixed slurries were sealed with butyl rubber stoppers and repeatedly flushed with N₂ to establish anoxic conditions. Prior to the NaCl treatments, slurries were pre-incubated for 7 days in the dark at 30°C. Then, NaCl was dissolved in 5 mL autoclaved water to adjust slurries to a concentration of 0, 300, and 600 mM NaCl corresponding to 0, 17.5, and 35 g/L. The slurries were shaken vigorously and then incubated for 2 days without shaking. After the NaCl treatment, the unstirred 2-cm water layer and the rice straw were withdrawn and discarded. The slurry material beneath the top layer was collected for molecular analysis and immediately shock-frozen using liquid nitrogen and then kept at -80°C until further processing. Gas samples (0.2 mL) and liquid samples (0.5 mL) were taken from each slurry microcosm at three time points for process measurements: 0 and 7 days after flooding (pre-incubation period); and 2 days after the addition of NaCl. Molecular analyses were performed for control (0 mM NaCl) and the two salt treatments (300 and 600 mM NaCl) after a total incubation period of 9 days (7-day pre-incubation followed by 2-day stress exposure). Independent analysis of triplicate slurries for each experimental treatment included process measurements, qPCR and RT-qPCR, as well as metatranscriptomic analysis of both 16S rRNA and enriched mRNA. Taxonomic profiling of 16S rRNA was based on metatranscriptomic analysis of total RNA. Given the relatively low amounts of mRNA in total RNA (see below), stress-induced changes in the taxonomic and functional mRNA profiles, however, were assessed using enriched mRNA.

Methane production was determined by analysis of headspace samples using a gas chromatograph equipped with a Sephadex^® column and a methanizer (Conrad and Klose, 1999). The availability of acetate was monitored by ion exclusion HPLC (Conrad and Klose, 2006) applying an Aminex^® HPX87H organic acid analysis column (Bio-Rad, Munich, Germany). Obtained data were analyzed using PeakSimple (SRI Instruments, Bad Honnef, Germany) and Clarity (DataApex, Prague, Czechia).

Total RNA and enriched mRNA were extracted from slurries using a previously described method (Mettel et al., 2010). Six aliquots (6 × 0.5 g of soil) were extracted from each replicate slurry for analysis of total RNA. The six RNA extracts were pooled, treated with DNase I (Ambion, Austin, TX, USA), and purified using the RNA Clean and Concentrator kit (Zymo Research, Irvine, CA, USA) according to the manufacturer’s instructions. Enriched mRNA was obtained by subjecting total RNA (>2.5 μg) to subtractive hybridization using the Ribo-Zero^TM Magnetic Kit for Bacteria (Epicentre, Madison, WI, USA). Generally, 16 extracts of total RNA (each from 0.5 g of soil) were obtained from each replicate slurry and purified using the RNA Clean and Concentrator kit (Zymo Research) according to the manufacturer’s instructions. These 16 RNA extracts were pooled prior to quantitative removal of rRNA by subtractive hybridization. Quality and integrity of the purified RNA and enriched mRNA were checked by Bio-Rad Experion^TM and RNA HighSens Chips (Bio-Rad, Hercules, CA, USA). Genomic DNA was extracted using the FastDNA^® SPIN kit for soil (MP Biomedicals, Eschwege, Germany).

Bacterial and archaeal 16S rRNA gene and transcript numbers per gram of dry soil were determined by qPCR and RT-qPCR, respectively (Wegner and Liesack, 2016). In addition, numbers of genes and transcripts encoding a key enzyme of methanogens, methyl coenzyme reductase A (mcrA), were quantified. Total RNA was randomly reverse transcribed using the GoScript Reverse Transcription System (Promega, Mannheim, Germany). Standard curves were constructed using genomic DNA from Escherichia coli (bacterial 16S rRNA genes/transcripts, calibration range from 10 to 10⁸ copies) and PCR amplicons (archaeal 16S rRNA genes/transcripts, calibration range from 10 to 10⁸ copies; and mcrA genes/transcripts, calibration range from 10 to 10⁹ copies). The analyses were performed using SYBR Green-based assays (Stubner, 2004; Kemnitz et al., 2005; Angel and Conrad, 2013). qPCR was carried out using a CFX Connect Real-Time PCR detection system (Bio-Rad). All reactions were carried out in three technical replicates. The PCR efficiency was at least 85% (R² > 0.99). The absence of non-specific products was confirmed by melt curve analysis.

In total, nine samples each from total RNA and enriched mRNA were subjected to cDNA synthesis using the NEBNext^® Ultra^TM Directional RNA Library Prep Kit for Illumina^® (New England Biolabs, Ipswich, MA, USA) according to the manufacturer’s instructions. cDNA yield and quality were determined using DNA HighSens Chips (Bio-Rad). The 18 cDNA libraries were sequenced at the Max Planck Genome Centre Cologne using 250-bp Illumina MiSeq.

The raw Illumina paired-end reads were trimmed by removal of residual adaptor sequences using Cutadapt (Martin, 2011) and filtered using a minimum mean quality score of 20. 16S rRNA reads were extracted from total RNA datasets by SortMeRNA 2.0 (Kopylova et al., 2012). For enriched mRNA, reads mapping to rRNA and non-coding small RNA were filtered out by SortMeRNA 2.0, using SILVA and RFAM as reference databases. The remaining reads were considered to be mRNA.

Analysis of 16S rRNA (total RNA) was carried with a customized pipeline that applies USEARCH (Edgar, 2013) and QIIME (Caporaso et al., 2010) as overall framework. OTU clustering was done using the UCLUST algorithm integrated in USEARCH and an OTU threshold of 95%. Using a custom PYTHON script, the resulting OTU table was converted to match the particular specifications of QIIME. Taxonomic assignments to the OTU table were carried out using the RDP classifier method (Wang et al., 2007) implemented in QIIME with SILVA SSU (release: 119; Quast et al., 2013) as the reference database. A total of 5,993,643 quality-filtered reads were obtained from total RNA, of which 3,098,833 reads were derived from 16S rRNA (Supplementary Table S1).

Taxonomic assignment and functional annotation of enriched mRNA was carried using the UBLAST algorithm implemented in USEARCH applying an e-value cutoff of 1e-5 for database searches against NCBI’s non-redundant protein database (Pruitt et al., 2012). MEGAN5 (Huson et al., 2011) was used for parsing and downstream analysis of the UBLAST output. Gene information identifier mapping and tabular BLAST output files were loaded to MEGAN. Transcripts that had a hit in the NCBI’s nr protein database were functionally classified into SEED subsystems. A total of 1,607,909 reads had a taxonomic hit in this database (Supplementary Table S2).

Significant (P < 0.05) differences in process measurements and RT-qPCR data between control (0 mM NaCl) and the NaCl treatments (300 and 600 mM NaCl) were determined by analysis of variance (ANOVA). Statistical Analysis of Metagenomic Profiles (STAMP 2.0.8; Parks et al., 2014) was used to identify biologically relevant differences within the (i) 16S rRNA and (ii) mRNA datasets. Significant (P < 0.05) differences were determined by ANOVA and subsequently corrected for multiple tests using Benjamini–Hochberg false discovery rate (Q-value; Benjamini et al., 2001).

The sequence data have been submitted to GenBank under ID BioProject numbers PRJNA313183 (total RNA) and PRJNA302859 (enriched mRNA). This involves the accession numbers SRP077890 (total RNA) and SRP066471 (enriched mRNA).

Acetate started to accumulate at pretreatment day 1 and reached its transient peak concentration of 11.9 ± 2.3 mM at pretreatment day 7. In the control treatments, the time period between days 7 and 9 was defined by high CH₄ production (from 3.8 ± 1.0 to 16.3 ± 1.1 mM), in correspondence to a strong decline in the acetate concentration (from 11.9 ± 2.3 to 5.2 ± 0.6 mM) (Figure 1). This is in good agreement with previous research (Wegner and Liesack, 2016). Salt stress inhibited CH₄ production, thereby leading to an accumulation of acetate. There was a significant correspondence between the inhibition effect and severity of salt stress (P < 0.01; Figure 1). Salt stress reduced net consumption of acetate by 55.3% (300 mM NaCl) and 70.9% (600 mM NaCl). Correspondingly, methane production was reduced by 48.9% and 68.5%, respectively. The corresponding values for reduced CO₂ production were 56.5% and 82.4%.

The strong inhibition of CH₄ production agrees well with previous results about how salinization affects methanogenic communities in natural freshwater wetlands (Baldwin et al., 2006; Chambers et al., 2011; van Dijk et al., 2015) and paddy soils (van der Gon and Neue, 1995). Overall, these studies also observed a negative correlation between salinity and methane production. Similar to our results, the incubation of freshwater wetland soil with full-strength seawater salinity (35 g/L NaCl) led, relative to freshwater controls, to an inhibition of CH₄ production by approximately 70% within the first few days of incubation (Chambers et al., 2011).

Salt stress induced a significant decrease (P < 0.05) in the 16S rRNA gene and transcript numbers of both bacteria and archaea between control and 600 mM NaCl. The decrease was in the range of up to half an order of magnitude. Stress also induced a significant decrease in the mcrA gene and transcript numbers. In 600 mM NaCl, copy numbers were decreased (P < 0.05) by 34% (gene) and 59% (transcripts) relative to the control (Supplementary Figure S1).

The percentage proportion of mRNA in total RNA significantly decreased (P < 0.005) with increasing salt stress, from 4.22 ± 1.15% (control) to 2.15 ± 0.23% (300 mM NaCl) and 1.54 ± 0.11% (600 mM NaCl) (Supplementary Table S1). Relative to the control (100%), the mRNA proportion thus decreased by ∼50% and ∼65%, respectively.

The stress-induced decline in process activities corresponds well to the decrease in absolute transcript abundance (16S rRNA, mcrA) but, even more pronounced, to the decrease in the mRNA proportion in total RNA. Apparently, this proportion is an excellent proxy for changes in the community-wide expression level of functional genes. As such, it also seems to be a suitable molecular indicator for changes in process activities. There was a nearly exact correspondence between the stress-induced decrease in net acetate consumption, methane production, and the proportion of mRNA in total RNA [∼50% (300 mM NaCl) and 65–70% (600 mM NaCl)]. Obviously, mRNA turnover is more responsive to environmental change and adverse conditions than rRNA dynamics, in good correspondence to the mRNA half-lives that are known to be for bacteria in the order of minutes and for methanogenic archaea in the range of tens of minutes and a few hours (Evguenieva-Hackenberg and Klug, 2009; Peterson et al., 2016). The strong decrease in functional gene expression, however, was related to tremendous changes in the taxonomic and functional profiles of enriched mRNA.

Salt stress induced significant changes in the microbial community metatranscriptome of total RNA (16S rRNA) and enriched mRNA (Figures 2A,B). Except Chloroflexi, all phylum-level groups showed significant changes in their relative rRNA abundance. Acidobacteria and Planctomycetes were the only phylum-level groups whose relative mRNA abundance decreased significantly by 77.0% (P < 0.001) and 51.0% (P < 0.01), respectively (Supplementary Figure S2).

On family level, changes in relative taxon abundance were greater for 600 mM NaCl than for 300 mM NaCl stress exposure and clearly more pronounced in mRNA than rRNA analysis (Figure 3A and Supplementary Figure S3). The stress-induced response, however, greatly varied between the different family-level groups. Among bacteria, Clostridiaceae were most competitive and the only bacterial group that significantly increased in relative mRNA abundance, from 8.7% (control) to 19.8% (600 mM NaCl; P < 0.05). By contrast, their contribution to total 16S rRNA did not significantly change. Lachnospiraceae were able to combat salt stress at 300 mM NaCl. Their relative mRNA abundance increased significantly (P < 0.01), but declined when exposed to 600 mM NaCl (Figure 3A). The Peptococcaceae, Veillonellaceae, and Syntrophomonadaceae were unable to compete, even at 300 mM NaCl. The relative mRNA abundance of the Ruminococcaceae was fairly stable at 300 mM NaCl, but significantly decreased at 600 mM NaCl. Bacillaceae showed no significant abundance change, neither on rRNA nor on mRNA level (Supplementary Figure S3).

Among the methanogenic archaea, Methanosarcinaceae showed a highly competitive response to salt stress relative to other family-level groups. Their relative rRNA (15.9–17.7%) and mRNA (37.6–47.9%) abundances were by trend, but not significantly, increased in response to salt stress (Figure 3A). The Methanosarcinaceae are known to harbor halotolerant and halophilic members (Kendall and Boone, 2006) and to outcompete Methanosaeta spp. at high acetate concentrations. The ability to respond quickly with physiological activity and growth to increased acetate levels is linked to three rrn gene copies in all Methanosarcina spp. versus one or two in Methanosaeta spp. (Stoddard et al., 2015). The early stages of rice straw decomposition are characterized by a rapid but transient increase in acetate levels (Figure 1). Thus, the genetic and metabolic potential coupled with high substrate availability may be the major reason why the Methanosarcinaceae were able to respond to severe salt stress more effectively than most other family-level groups and to maintain their relative dominance in the community-wide mRNA pool (Figure 3A).

The most striking effect of salt stress was observed for the Methanocellaceae. While their contribution to total 16S rRNA only decreased from 1.1% (control) to 0.7% (600 mM NaCl), their relative mRNA abundance strongly declined from 7.4% (control) to 1.3% already when exposed to 300 mM NaCl (Figure 3A). Their unability to cope with salt stress agrees well with the low salt tolerance of taxonomically described Methanocella spp.: 17 mM NaCl (M. paludicola; Sakai et al., 2008); 86 mM NaCl (M. conradii; Lü and Lu, 2012); and 340 mM NaCl (M. avoryzae; Sakai et al., 2010). M. avoryzae would be able to tolerate 300 mM NaCl but has a temperature growth range of 37–55°C. Significant changes on both rRNA and mRNA level were also evident for the Methanosaetaceae and Methanobacteriaceae (Supplementary Figure S3).

The bacterial single-copy gene encoding RNA polymerase β subunit (rpoB) was shown to provide a phylogenetic resolution comparable to or even higher than the bacterial 16S rRNA gene (Case et al., 2007; Vos et al., 2012). We extracted the total number (2,247) of bacterial rpoB transcripts and archaeal homologs from the mRNA datasets obtained for control and NaCl treatments. The taxon-specific changes in rpoB abundance of the Clostridiaceae, Lachnospiraceae, Methansarcinaceae (at 600 mM NaCl), and Methanocellaceae were highly consistent with those observed for these families in total mRNA. The exactly same trend in taxon-specific abundance changes was also evident for mRNA encoding ribosomal proteins (Figure 3B).

The functional profiles of enriched and taxonomically assigned mRNA were analyzed according to SEED subsystems. The profiles clearly differed between control and NaCl treatments, even at SEED level 1 (Figure 2C). The number of categories whose transcript abundance was significantly affected by salt stress changed with SEED level and significance of the P and Q values (Supplementary Table S3). We first generated functional profiles across all SEED level 1 (Supplementary Figure S4) and 2 categories (Supplementary Figure S5). Based on that, 21 SEED level 2 categories were selected for taxon-specific analysis (Figure 4 and Supplementary Figure S6). These 21 categories were most significantly affected by salt stress and covered 37.1–49.7% of the enriched mRNA that could be functionally annotated (Supplementary Table S4). On family level, the stress-induced changes in the functional profiles corresponded well to the taxonomic shifts observed for these groups in total mRNA (Figures 3, 4 and Supplementary Figure S6). This was most evident for the Clostridiaceae, Lachnospiraceae, Methanosarcinaceae, and Methanocellaceae (Figure 4). Relative to the other family-level groups, clostridial mRNA was enriched in all the SEED level 2 categories (except for “ABC transporters”). Transcripts affiliated with Lachnospiraceae were enriched in various functional categories at 300 mM NaCl, but not at 600 mM NaCl. Among the Methanosarcinaceae, transcripts involved in “respiration,” “methanogenesis,” “glutamine, glutamate, aspartate, asparagine; ammonia assimilation,” “protein folding,” and “heat shock” were significantly enriched in response to salt stress. Greatest changes in mRNA abundance were related to the categories “respiration” and “methanogenesis.” Their relative transcript level increased from 5.0 ± 0.6% and 13.5 ± 1.1% (control) to 10.1 ± 1.2% and 21.0 ± 1.1% (600 mM NaCl), respectively. The relative transcript levels of all other functional categories either did not change or significantly declined. Transcripts affiliated with Methanocellaceae were underrepresented in all the 21 SEED level 2 categories (Figure 4). The mRNA abundance of the Lachnospiraceae and Methanocellaceae was insufficient for a more detailed gene expression analysis. Therefore, subsequent data presentation is only focused on the two most abundant and competitive families: Clostridiaceae and Methanosarcinaceae.

Using Illumina RNA-Seq, no pure culture study has yet been performed for members of the Clostridiaceae and Methanosarcinaceae to assess their transcriptome response to salt stress. Therefore, for these two family-level groups, we expanded our data analysis to the expression of those pathways, functions, and genes that were most responsive to the NaCl treatments (Tables 1, 2). This response was driven by the genera Clostridium and Methanosarcina. Nearly all the mRNA reads assigned to the Clostridiaceae and Methanosarcinaceae were affiliated with these two genera.