Chronic Physical Disturbance Substantially Alters the Response of Biological Soil Crusts to a Wetting Pulse, as Characterized by Metatranscriptomic Sequencing.

Biological soil crusts (biocrusts) are microbial communities that are a feature of arid surface soils worldwide. In drylands where precipitation is pulsed and ephemeral, the ability of biocrust microbiota to rapidly initiate metabolic activity is critical to their survival. Community gene expression was compared after a short duration (1 h) wetting pulse in both intact and soils disturbed by chronic foot trampling. Across the metatranscriptomes the majority of transcripts were cyanobacterial in origin, suggesting that cyanobacteria accounted for the bulk of the transcriptionally active cells. Chronic trampling substantially altered the functional profile of the metatranscriptomes, specifically resulting in a significant decrease in transcripts for nitrogen fixation. Soil depth (biocrust and below crust) was a relatively small factor in differentiating the metatranscriptomes, suggesting that the metabolically active bacteria were similar between shallow soil horizons. The dry samples were consistently enriched for hydrogenase genes, indicating that molecular hydrogen may serve as an energy source for the desiccated soil communities. The water pulse was associated with a restructuring of the metatranscriptome, particularly for the biocrusts. Biocrusts increased transcripts for photosynthesis and carbon fixation, suggesting a rapid resuscitation upon wetting. In contrast, the trampled surface soils showed a much smaller response to wetting, indicating that trampling altered the metabolic response of the community. Finally, several biogeochemical cycling genes in carbon and nitrogen cycling were assessed for their change in abundance due to wetting in the biocrusts. Different transcripts encoding the same gene product did not show a consensus response, with some more abundant in dry or wet biocrusts, highlighting the challenges in relating transcript abundance to biogeochemical cycling rates. These observations demonstrate that metatranscriptome sequencing was able to distinguish alterations in the function of arid soil microbial communities at two varying temporal scales, a long-term ecosystems disturbance through foot trampling, and a short term wetting pulse. Thus, community metatranscriptomes have the potential to inform studies on the response and resilience of biocrusts to various environmental perturbations.

Introduction

Drylands are defined by their lack of a vital resource, precipitation. Thus, water availability is a significant force in shaping the structure and function of drylands or arid ecosystems. Biological soil crusts (biocrusts) are a community of cyanobacteria, algae, associated heterotrophic bacteria, fungi and mosses that colonize arid soils worldwide (Belnap, 2003; Pointing and Belnap, 2012; Belnap et al., 2016). These communities inhabit plant interspaces and are important players in fixing carbon and nitrogen in dryland soils (Evans and Lange, 2001; Barger et al., 2016). Because of the unpredictable periodicity and amount of precipitation, biocrusts respond quickly to wetting events to exploit the available moisture to repair cell damage and synthesize new biomass. For example, when precipitation wets surface soils, biocrust cyanobacteria display a property of “hydrotaxis” and move toward the available water (Garcia-Pichel and Pringault, 2001; Pringault and Garcia-Pichel, 2004). Pigments, such as chlorophyll a, are rapidly produced, resulting in a “greening” of the surface soils (Abed et al., 2014). Concurrently, consumption of oxygen and evolution of CO2 by the biocrusts can be detected, indicating vigorous metabolic activity (Rajeev et al., 2013; Tucker et al., 2017). Wetting pulses can also stimulate “blooms” of specific bacteria, highlighting that specific taxa have differing thresholds and responses to wetting (Karaoz et al., 2018). Despite the multitude of metabolic and physiological shifts associated with wetting in biocrusts, the specific genes and biochemistry underlying these processes are largely unknown.

The dominant cyanobacterium of biocrusts of the Colorado Plateau is Microcoleus vaginatus. RNA microarrays specific to M. vaginatus, have investigated the transcriptional response of this organism to wetting cycles in laboratory microcosms. These experiments show that M. vaginatus rapidly resuscitates after the addition of water with an upregulation of genes for photosynthesis, anabolic metabolism, and DNA repair (Rajeev et al., 2013). Similarly, the desert moss Syntrichia caninervis increases gene expression for photosynthesis, membrane integrity, and protein synthesis upon wetting (Gao et al., 2014). However, these results consider biocrust species in isolation and not the community as a whole. Thus, they cannot be scaled to address ecosystem interactions, such as trophic exchanges among biocrust members.

In this study, metatranscriptome sequencing was employed to characterize community gene expression. Biocrusts and below-crust soils were collected from undisturbed sites for comparison to soils that have experienced annual physical trampling for 17 years. This study was designed to address three hypotheses. First, the activity of surface biocrusts would differ from below-crust soils, driven by gradients in factors such as light and moisture that vary largely at small vertical scales. Second, a short duration wetting pulse of 1 h will be sufficient to capture the initial response of the arid soil communities to a water pulse through activities such as the initiation of photosynthesis and carbon fixation. Finally, a severe ecosystem disturbance in the form of chronic trampling will significantly alter the functional profile of the soil community and possibly influence the community response to wetting.

Materials and Methods

Field Site

Samples were collected from the Island of the Sky district of Canyonlands National Park, near Moab, UT, United States (WGS84, 38°27′27.21″N, 109°32′27.93″W). The study site was located in Coleogyne ramosissima (blackbrush) shrubland and harbored biocrusts that differed in development and composition of lichens and cyanobacteria. This site has been employed in a chronic trampling experiment and has been described in detail previously (Kuske et al., 2012; Steven et al., 2015).

The sampling date was October 26th, 2012, the most recent previous precipitation event was 2 days earlier, registering as a trace event (too small to measure) and the year to date precipitation was 6.6 cm, about 35% of the normal year to date average; so the biocrusts were dry at the time of sampling. Samples were collected from the plant interspaces from areas visibly devoid of mosses and lichen to focus on cyanobacteria-dominated biocrusts. The biocrusts were defined as the layer of soil bound together by the cyanobacterial filaments (generally ca. 1 cm depth) and the below-crust soils were the unbound soils directly below at about 2 cm depth.

Wetting Treatment

For wetting, three soil cores, approximately 13 mm in diameter, of ca. 2 cm depth were collected and placed in a Petri dish and 25 ml of sterile distilled water was added from the top, simulating a rainfall. This volume of water is approximately equivalent to a 5 mm precipitation event (based on the volume of the Petri dish). This is within normal precipitation amounts for the region (Hereford and Webb, 1992). The purpose of placing the sample in the Petri dish was to limit the lateral diffusion of water to maintain a consistent wetting treatment between samples. While the Petri dish would not allow for drainage, given the volume of water compared to the volume of biocrust (5 mm water depth compared to 2 cm soil depth) it is unlikely that a significant amount of drainage would have occurred, particularly in the short duration (1 h) of the wetting pulse.

Sample Collection

Dry samples (ca. 10 g) were collected and immediately placed separately into a 50 ml solution of RNAlater (Sigma-Aldrich) and directly stored on dry ice. For the wetted samples, soil cores were collected and placed in Petri dishes, wetted, and after 1 h of incubation in the field, the biocrust and below-crust samples were separated and individually placed in RNAlater and stored on dry ice. After collection, samples were transported frozen to Los Alamos National Laboratory, Los Alamos, NM, United States, for processing.

RNA Extraction and Library Preparation

Prior to RNA extraction, each sample was thawed and 4 g of material was removed. RNA was purified using the RNeasy purification kit (Qiagen), and eluted in a final volume of 50 μl. RNA concentration was determined using a Nanodrop spectrophotometer (Thermo Scientific). DNA was removed using Turbo DNase I enzymatic digestion (Life Technologies). A quantitative real-time PCR (qPCR)-based assay was used to insure that residual DNA was not present in the purified RNA preparations.

Ribosomal RNA (rRNA) was removed from the total RNA via hybridization to biotinylated probes provided in Ribo-Zero kits (EpiCentre). Metabacteria and eukaryote (human/mouse/rat) probes were combined, hybridized to total RNA, and the complexes were removed with the magnetic beads. Samples were processed according to the manufacturer’s protocols.

RNA samples were quantified using a Life Technologies RNA quantification kit and a Qubit 2.0 instrument (Life Technologies). The size distribution of the RNA molecules in each sample was determined using a BioAnalyzer 2100 and RNA pico chips and reagents (Agilent). A maximum of 50 ng of RNA from each sample was converted to an Illumina sequencing library using ScriptSeq kits (EpiCentre). Each library was generated using a unique 6 bp multiplex barcode provided in the ScriptSeq kit. Following library generation, the 12 samples were pooled in equimolar concentrations and run on two lanes of the Illumina HiSeq2000 sequencing platform at the Los Alamos National Laboratory sequencing center, with 2x150 sequence chemistry. Raw sequences were uploaded to the NCBI Sequence Read Archive (SRA) under the Submission ID: SUB3320760 and BioProject ID: PRJNA421954.

Sequence Quality Trimming and Assembly

Raw sequence reads were quality screened and contaminating sequences were removed with the BBDuk package (v. 37.17) within the BBTools suite[1]. Briefly, sequence adaptors were removed using the provided adaptors.fa reference file with a minimum kmer match of 11 bp and a 1 bp mismatch allowed for removal. The sequences were then quality trimmed to Q10 using the Phred algorithm with a kmer size of 23 (including the “tbe” and “tbo” flags). Small and large subunit ribosomal RNA sequences that passed the magnetic bead pulldown were bioinformatically removed using BBDuk and the ribokmers.fa reference file (resulting in a removal of 5 to 22% of the sequence reads). Finally, common sequence contaminants were removed with the hg19_main_mask_ribo_animal_allplant_allfungus.fa reference file. Both sequence screenings were performed with a kmer size of 31.

The quality-filtered sequences were pooled and assembled with MEGAHIT using default metagenomic parameters for paired-end reads (Li et al., 2016). Only contigs greater than 500 base pairs were retained. The number of contigs in the assembly was 209,032 ranging in length from 500 to 35,855 bp The N50 of the assembly was 835 just below the predicted average length of a bacterial mRNA of 924 bp (Xu et al., 2006). These data indicate that many potential full-length sequences were recovered. Additionally, contigs > 10,000 base pairs suggest that even operon length transcripts were assembled. The assembled FASTA file will be made available by the authors upon request.

Estimation of Transcript Abundance and Taxonomic/Functional Annotation

The quality filtered raw reads were pseudo-aligned to the assembly to generate transcript abundances in the samples using the kallisto software package (v0.43.1; Bray et al., 2016). First, the assembly was used to generate an index file, using the default parameters (31 kmer length) and the quality trimmed raw read files were mapped to the index file. One hundred bootstrap analyses were performed for pseudo replication and import into the Sleuth software package to generate contig count tables (Pimentel et al., 2017).

For sequence annotation contigs were aligned to the NCBI non-redundant (NR) protein database (downloaded July 9th, 2017) using DIAMOND with default parameters (Buchfink et al., 2014). Alignments were uploaded to the community version of MEGAN 6 (Huson et al., 2016) and contigs were assigned to SEED subsystems (Overbeek et al., 2014) within MEGAN using the long-read enabled sequence assignment and the May-2015 accession to SEED mapping file. The LCA classification of reads was performed with the prot.accession2taxid mapping file. Merging the contig taxonomic and hierarchical functional annotations with the estimated counts generated a table with LCA classifications and hierarchical SEED annotations for each contig in the assembly and was the raw data on which statistical comparisons were made.

Statistical Analyses

Non-metric multidimensional scaling (NMDS) analyses were performed with the phyloseq R package (McMurdie and Holmes, 2013). Inter-sample distances were calculated with the Bray–Curtis dissimilarity metric, and statistical differences in composition were assessed using PERMANOVA with the “adonis” function in the vegan package for R (Dixon and Palmer, 2003).

To identify significant differences in sequences assigned to the same SEED categories (annotation bins), the annotated count tables were uploaded to the STAMP software package (Parks et al., 2014) along with a relevant metadata table that encoded the information about each sample. For comparisons between multiple groups, differentially abundant bins were identified with an ANOVA test to determine if there was a significant difference in means between the groups, with a p-value cutoff of ≤0.05 and an effect size filter of 0.5. All p-values were corrected for multiple testing using the Bonferroni method. For two group comparisons (i.e., between control/trampled and surface/subsurface) significant differences were identified by Welsh’s t-test with Storey’s FDR method for multiple test correction. We additionally tested for an interaction between wetting treatments and trampling, or between biocrust and below-crust soils and did not identify any significant interactions for any of the differentially abundant bins.

To identify the specific contigs that were differentially abundant due to wetting, unnormalized estimated count data was analyzed with the phyloseq implementation of DESeq2 (Love et al., 2014). Initially, contigs with a summed total count of less than 10 across the samples were removed to avoid the analysis of extremely rare contigs. Differentially abundant contigs were identified with the Wald test on data with a parametric fit. Contigs were considered significantly different in estimated counts if the Benjamin Hochberg multiple test corrected p-values were ≤0.05 with a false detection rate of 0.1. Only contigs showing at least a twofold change in mean estimated counts were considered biologically meaningful.

Results

Clustering of Sequence Datasets

Differences between the biocrusts and trampled soils were readily apparent at a visual scale (Figure ). To investigate the relationships between the metatranscriptomes from the different samples, the transcript composition and abundance were employed to generate a NMDS plot (Figure ). The control and trampled samples clustered distinctly (PERMANOVA p = 0.001), indicating that chronic trampling significantly altered transcript composition (Figure , large ellipses). The wetting treatment also induced an alteration in transcript composition, as evidenced by the clustering of the wet samples from the dry samples (PERMANOVA p = 0.001 for biocrusts, p = 0.004 for trampled soils), indicating that gene expression was significantly altered within 1 h of a water pulse (Figure , inset ellipses).

(A) Photograph of field plots. The control plots show well-developed biocrusts directly adjacent to a trampled plot (separated by a nylon rope). A pen is included for scale. (B) NMDS plot showing the relationship of metatranscriptome datasets. Intergroup distances were calculated with the Bray–Curtis metric. The colored ellipses are included as a visual aid only and were not part of the analysis. The groupings within each ellipse are indicated by the bolded text. The stress value of the ordination is indicated in the top left corner.

Taxonomic Bins in Metatranscriptomes

Each contig in the assembly was assigned to a taxonomic rank using a lowest common ancestor (LCA) approach, and the estimated counts were used to identify the most common taxa in the datasets. As the LCA method identifies the common taxonomic rank that is shared by all homology matches, sequences are often assigned to broad taxonomic ranks. For example, among the most abundant taxonomic bins in the contigs, were groups such as Cyanobacteria, Bacteria, and cellular organisms (Figure ).

Taxonomic bins in the metatranscriptomes. Contigs were assigned a lowest common ancestor (LCA) and binned based on these assignments. Estimated counts of the contigs were then employed to assign a percent of total estimated counts for each LCA bin, of which the 20 most abundant are displayed. Each point represents the mean abundance from three replicate samples.

Three taxonomic bins accounted for a substantial proportion of the estimated counts, namely the Oscillatoriales an order encompassing many filamentous Cyanobacteria, the phylum Cyanobacteria, and the kingdom Bacteria (Figure ). The next most abundant bin was Microcoleus vaginatus the specific cyanobacterium that is considered to be a keystone species of biocrusts. Thus, these data suggest that the bulk of RNA in biocrust soils is derived from Cyanobacteria, including the category Bacteria, which could also include sequences derived from Cyanobacteria.

Functional Annotations in the Metatranscriptomes

Contigs were functionally annotated to the SEED hierarchy, and the estimated counts were used to investigate the abundance of different functions in the datasets. At the broadest annotations levels the categories Protein Metabolism and Carbohydrates were dominant (Figure ). However, for the biocrusts the order was reversed with Carbohydrates-related reads being more abundant than for Protein Metabolism (Figure ).

SEED functional bins in the metatranscriptomes. Contigs were annotated within the SEED hierarchy and binned based to the broadest categorical assignments (level 1). Estimated counts of the contigs were then employed to assign a percent of total estimated counts for each bin, of which the 20 most abundant are displayed. Each point represents the mean abundance from three replicate samples.

To further investigate the most abundant metabolic pathways in the metatranscriptomes, the sequences in the categories Protein Metabolism and Carbohydrates were annotated to the second level in the SEED hierarchy (Figure ). Within the category Protein Metabolism the most abundant subsystems were largely involved in structural components of the ribosome or protein synthesis. However, the subsystems also included ATP dependent proteolysis and Bacterial proteasomes, indicating a capacity for protein recycling. Only a single subsystem (ribosomal protein S12p) was found to have a significant difference in abundance among the samples.

Subsystems within the category Protein Metabolism. The 20 most abundant level 2 subsystems within the category are shown. Sample abbreviations: CD, crust dry; CW, crust wet; BCD, below-crust dry; BCW, below-crust wet; TSD, trampled surface dry; TSW, trampled surface wet; TBD, trampled below dry; TBW, trampled below wet. The size of the points represent the proportion of estimated counts as a percent of the total reads in the parent subsystem (i.e., percent of Protein Metabolism reads). Each point represents the mean of three replicate samples. Points colored red were significantly different by an ANOVA difference of means test.

The most abundant subsystems in the category Carbohydrates were related to carbon fixation (Figure ; Calvin-Benson cycle, and Alpha/beta carboxysomes) supporting that these soil communities primarily rely on autotrophic carbon fixation. Two central metabolic pathways, the Entner–Doudoroff pathway and TCA cycle were also represented among the most abundant subsystems. Within this category were also subsystems that point to the metabolism of specific carbon compounds, namely maltose, glycerol, and glycogen (Figure ). Among the Carbohydrate subsystems, three were identified as significantly different among the metatranscriptomes. Beta-carboxysomes were most abundant in the trampled subsurface soils. Beta-type carboxysomes tend to be more common in heterotrophic bacteria (Kerfeld and Melnicki, 2016), suggesting potential taxonomic shifts in the organisms performing carbon fixation in the trampled subsurface soils. In comparison, Maltose utilization and Photorespiration were both elevated in the wet biocrusts, suggesting increased expression under wetting. Yet, the most abundant subsystems in the datasets tend to be present in similar proportions across the datasets.

Subsystems within the category Carbohydrates. The 20 most abundant level 2 subsystems within the category are shown. Sample abbreviations: CD, crust dry; CW, crust wet; BCD, below-crust dry; BCW, below-crust wet; TSD, trampled surface dry; TSW, trampled surface wet; TBD, trampled below dry; TBW, trampled below wet. The size of the points represent the proportion of estimated counts as a percent of the total reads in the parent subsystem (i.e., percent of Carbohydrate reads). Each point represents the mean of three replicate samples. Points colored red were significantly different by an ANOVA difference of means test.

Differentially Abundant Subsystems due to Trampling

Chronic physical trampling induced a significant shift in composition of the metatranscriptomes as evidenced by the NMDS analysis (Figure ). To identify functions underlying the differences, contig annotations were binned to the second level in the SEED hierarchy and differentially abundant functional bins were identified (Figure ). A total of 15 functional bins were identified as significantly different. The largest proportional difference was in nitrogen fixation, which was significantly reduced in the trampled soils (Figure ). Concurrently, sequence reads to deal with nitrosative stress, the ability to detoxify reactive nitrogen species, was also reduced in the trampled soils, supporting a shift of nitrogen metabolism. Other subsystems point to alterations in amino acid metabolism, specifically methionine and lysine. Two subsystems were also related to the transport and maintenance of metal levels in cells, with a decrease in genes for copper homeostasis as well as nickel and cobalt transport (Figure ). Thus, these observations support that chronic trampling alters the metabolism and transport of nitrogen and metals by the soil community.

Differentially abundant subsystems due to trampling. Boxplot of differentially abundant subsystems between biocrust and trampled soils. For this analysis the wet and dry, as well as the surface and subsurface samples were considered together. Thus, the boxes summarize the data for 12 replicate samples. For each box the mean is represented by the bolded bar and boxes represent the interquartile range. Outliers are shown as points outside of the box and whiskers.

Differentially Abundant Subsystems due to Soil Depth

The datasets were separated into surface and subsurface samples, irrespective of trampling, to investigate if specific subsystems differentiated the metatranscriptomes based on soil horizon. A total of five significantly different subsystems were identified between surface and subsurface soils (Figure ). Of the five subsystems, three were related to the biosynthesis or metabolism of amino acids, and one for the biosynthesis of a carbohydrate polymer, namely arabinogalactan, which consists of arabinose and galactose monosaccharides that can complex with peptidoglycan in the envelope of some bacteria (Crick et al., 2001). Finally, the subsystem for transport of the metals nickel and cobalt was significantly reduced in the subsurface samples, matching an observation for the trampled soils (Figure ). Overall, a relatively small number of subsystems discriminated the metatranscriptomes of surface versus subsurface soil communities.

Differentially abundant subsystems due to soil depth. Boxplot of differentially abundant subsystems between surface and subsurface soils. For this analysis the wet and dry, as well as the biocrust and trampled samples were considered together. Thus, the boxes summarize the data for 12 samples. For each box the mean is represented by the bolded bar and boxes represent the interquartile range. Outliers are shown as points outside of the box and whiskers.

Differentially Abundant Contigs due to Wetting

Each contig in the assembly was plotted based on its mean abundance and log-fold change between dry and wet samples (Figure ). These data clearly demonstrate that the number of contigs that were significantly altered due to wetting was larger in the biocrusts than for any of the other sample types (Figure , red points). In comparison, no contigs were identified as significantly enriched in the wet trampled surface soils, whereas 27 were significantly enriched in the dry conditions (Figure ). The number of differentially abundant contigs was more synonymous between the below-crust and trampled below soils (Figures versus ).

Bland–Altman plots displaying shifts in contig abundance due to wetting. The average contig abundance (log10 scale) is displayed on the x-axis. As the means are for both wet and dry samples, each point is the average of six samples. The fold-change (log2 scale) is displayed on the y-axis, such that contigs enriched in dry samples are shown below the 0 line and contigs enriched in wet are shown above. Points having a significant difference with False Discovery Rates (FDRs) of less than 0.1 and at least a twofold difference in abundance (indicated by the light gray dashed lines) are colored red. (A) Differentially abundant contigs in biocrusts. (B) Differentially abundant contigs in trampled surface soils. (C) Differentially abundant contigs in below-crust soils. (D) Differentially abundant contigs in trampled subsurface soils.

To investigate the functional role of the differentially abundant contigs, they were annotated to the second level of the SEED hierarchy and the results for the 20 most abundant bins in the biocrusts are summarized in Figure . The functional assignments of the contigs split the annotations into three groups: (1) Functions specifically enriched in dry crusts, which included 5-FCL-like proteins [a protein family implicated in thiamin metabolism (Pribat et al., 2011)] and hydrogenases (Figure , top two yellow bars). (2) Functions specifically enriched in wet biocrusts, Photosystems I and II (Figure , bottom two green bars). (3) Functions that included contigs enriched in both the dry and wet biocrusts, which encompassed the remaining functional bins (Figure and Supplementary Table ). Results for the surface trampled soils are not displayed as no contigs were significantly enriched in the wetted samples, and of the 27 significantly enriched in the dry samples, only four were functionally annotated, all to the category hydrogenases (Supplementary Table ).

Differentially abundant subsystems due to wetting in surface soil samples. The level 2 SEED annotations were used to functionally bin the contigs identified as differentially abundant (Figure ), of which the 20 most abundant (sorted by the total number of contigs in each bin) are displayed. A full list of differentially abundant contigs in the biocrusts is listed in Supplementary Table (only annotated contigs are presented). For the trampled surface samples, no contigs were identified as significantly enriched in the wetted samples (Figure ) and of the 27 significantly enriched in the dry samples only four were annotated to SEED subsystems, so the data is not shown. All four subsystems were in the category Hydrogenases. A full list of differentially abundant contigs in the trampled surface soils is listed in Supplementary Table .

The below-crust and trampled below-crust samples showed a similar number and pattern of differentially abundant contigs (Figures ). Similarities included hydrogenases, protein chaperones, pathogenicity islands, GroEL/GroES, Fe-S-cluster assembly, ATP dependent proteolysis, and biofilm formation all of which were enriched in the dry soils (Figure ). In sum, these observations suggest that while there were substantial differences in the response of the surface soils to the wetting pulse, within the subsurface the effects of wetting were more congruent.

Differentially abundant subsystems due to wetting in subsurface soil samples. The level 2 SEED annotations were used to functionally bin the contigs identified as differentially abundant (Figure ). (A) Below-crust soils, (B) below trampled soils. A full list of differentially abundant contigs in the below-crust soils is listed in Supplementary Tables for the below trampled soils.

Biogeochemical Cycling Genes in the Metatranscriptomes

The contig annotations were examined for genes in carbon and nitrogen cycling, of which a reference set were selected to investigate the effect of wetting on biogeochemical cycling (Table ). As can be seen from Table , multiple contigs were mapped to the same annotation, suggesting multiple sequence variants exist of each of the genes. The multiple transcripts for the same gene product were easily explained by multiple organisms encoding the same genes, multiple gene copies, or incomplete assembly of full-length transcripts.

The biogeochemical cycling genes were assessed for their abundance in the dry and wet biocrusts. The biocrusts were chosen as they showed the most pronounced response to the wetting pulse (Figure ). Nitrogen cycling genes in the biocrusts did not show a consensus response to wetting. The marker genes for each of the nitrogen cycling processes included contigs that were more abundant in wet and dry conditions, with the exception of nitrogen fixation (Figure ). Similarly, carbon cycling genes did not move in concert in response to wetting. In this regard, these data highlight the complexity of the defining a common response of the biocrust community to wetting and linking the changes in transcript abundance to the potential for biogeochemical cycling.

Biogeochemical cycling genes in the biocrust metatranscriptomes. The biogeochemical cycling genes listed in Table were mapped onto the biocrust datasets. Each point represents the average normalized abundance of each contig in the dataset (n = 6). Biogeochemical cycling genes are colored based on the biogeochemical pathway to which they belong.

Discussion

One observation from this study was that a relatively high proportion (48%) of the assembled transcript sequences had no significant homology to any sequences in the non-redundant protein database. Furthermore, for many of the transcripts the best taxonomic assignments were to broad groups such as Bacteria, Cyanobacteria, or Actinobacteria (Figure ). This limits the predictions we can make from homology mapping of transcripts. Several studies have found that combining metagenomic and metatranscriptomic reads can increase the assembly length and annotation of the resulting sequences (Franzosa et al., 2014; Ye and Tang, 2016). The improvements of multi-omic sequence assemblies arise from longer contigs giving more accurate annotations or larger contigs linking sequence reads to phylogenetically anchored genes (Narayanasamy et al., 2016). In this regard, matched metagenomes may have increased our ability to assign and interpret the transcriptome sequences.

Metabolic Activity in Arid Soils as Revealed by Metatranscriptome Sequencing

The most abundant functional bins across the metagenomes indicate that central pathways in protein and carbohydrate metabolism account for a relatively large proportion of the metatranscriptomes across the samples (Figure ). Broadly similar categories have been described by metagenomic sequencing of biocrust soils, suggesting a correspondence between the most abundant genes and transcripts (Steven et al., 2012, 2014, 2015).

The pathways represented among the subsystem Carbohydrates point to specific carbon metabolic processes that support these arid soil communities (Figure ). Pathways for the fixation of CO2 into biomass are dominant, highlighting the scarcity of organic carbon in these soils and the dependence on autotrophic carbon fixation (Klemmedson, 1989). Beyond carbon fixation, pathways for glycerol, glycogen and maltose metabolism were all abundant in the metatranscriptomes (Figure ). Cyanobacteria accumulate intracellular stores of glycogen as a carbon and energy source (De et al., 2017) and maltose can be employed by mixotrophic filamentous cyanobacteria as a heterotrophic carbon source (White and Shilo, 1975). A common strategy employed by microorganisms to survive desiccation is the intracellular accumulation of compatible solutes. These solutes replace water in the cell and aid in maintaining cellular and macromolecule structure (Esbelin et al., 2018). Compatible solutes can take the form of small sugars, glycocides, or peptides (da Costa et al., 1998). A common compatible solute is the sugar trehalose (Hershkovitz et al., 1991; Page-Sharp et al., 1999). M. vaginatus encodes genes to synthesize trehalose from both glycogen and maltose (Starkenburg et al., 2011). Thus, glycogen and maltose play a central role in protecting biocrust community members from starvation and desiccation. Thus, the metatranscriptomes point to a carbon-processing pathway that may be important in supporting the function of biocrusts.

Metatranscriptomic Characterization of an Ecosystem Disturbance

Biocrust soils subjected to chronic physical trampling were visually distinguishable from healthy biocrusts by the loss of mature developed biocrusts (Figure ). The metatranscriptomes of the trampled soils clustered distinctly by NMDS (Figure ), demonstrated several differentially abundant subsystems (Figure ), and displayed a notably altered response to a wetting pulse (Figure ). Taken together, these data clearly indicate that chronic physical trampling alters the structure of the biocrust soils and community gene expression.

Trampling of the biocrusts compromises the integrity of the biocrust structure, but some cyanobacteria survive the treatment and can be detected in the soils following the trampling (Kuske et al., 2012; Steven et al., 2015). In this respect, the trampled soils may be comparable to early successional stage biocrusts. In the initial stages of biocrust formation, the founder species M. vaginatus is dominant, forming the physical structure of the biocrust. Accordingly, sequences within the bin Oscillatoriales, the order to which M. vaginatus belongs, tended to be highest in the trampled datasets (Figure ). The dominance of filamentous cyanobacteria in early successional stage crusts comes at a cost. These cyanobacteria are not capable of nitrogen fixation. A lack of diazotrophic cyanobacteria in early and disturbed crusts is supported by reductions of nitrogen fixation potential as measured by acetylene reduction by up to 250% (Evans and Belnap, 1999; Zhou et al., 2016). However, the presence of non-cyanobacterial diazotrophs, predominantly within the phyla Firmicutes and Proteobacteria in early successional stage crusts has been presented as evidence that the nitrogen fixation potential of early successional biocrust may be under-estimated (Pepe-Ranney et al., 2016). The metatranscriptome data presented here suggest that nitrogen fixation potential is significantly reduced in trampled soils as assessed by the abundance of nitrogen fixation transcripts in the trampled soils (Figure ). However, these data only consider dry soils and a brief wetting pulse; perhaps the wetting threshold for activity of diazotrophic organisms in trampled soils is larger and thus requires a longer or larger wetting pulse. Yet, these data could also go some way to explain why biocrusts were generally unresponsive to inorganic nitrogen additions (Mueller et al., 2015). Perhaps a more robust response would be observed with nitrogen addition to disturbed or early-successional stage biocrusts, where nitrogen fixation potential may be limited.

In addition to reduced nitrogen fixation, the trampled surface soils also displayed a dramatic difference in the number of contigs that were altered in abundance due to wetting in comparison to the biocrusts (Figures versus 8B). No contigs were significantly enriched in the wet trampled surface soils in comparison to more than 4,000 for the biocrusts. Biocrusts produce an exopolysaccharide matrix and dark pigments that play a strong role in retaining moisture (Mazor et al., 1996; Rossi et al., 2012). Well-developed biocrusts also alter the soil surface properties, increasing roughness and forming pinnacles, properties that can additionally foster moisture retention and infiltration (Belnap, 2006). These properties are visually absent in the trampled soils (Figure ), presumably also affecting the ability of the community to take advantage of a wetting pulse. In this respect, these data suggest that not only are the trampled soil communities altered in their metabolic repertoire, but are less able to take advantage of the pulsed wetting events that control the productivity of these arid soils.

Differences Between Surface and Subsurface Metatranscriptomes

One of the unexpected findings of this study was the relatively small difference in the metatranscriptomes between the surface and subsurface samples. Because of gradients in resources such as light and moisture, biocrust communities are vertically stratified, often at the millimeter scale (Garcia-Pichel et al., 2003; Steven et al., 2013). Yet, only a small number of functional bins were identified as significantly different between the surface and subsurface samples (Figure ). The most parsimonious explanation for these results is that the cyanobacteria in both soil layers accounted for most of the transcriptome. The below-crust and trampled subsurface soils both showed a dominance of sequence belonging to cyanobacteria (Figure ), although metagenomic sequencing and other metrics indicate heterotrophic bacteria and Achaea are dominant in the subsurface soils (Steven et al., 2013). The high proportion of cyanobacteria in these soils can be easily explained by an incomplete separation of the biocrusts and the below-crust soils, or more likely that the small proportion of cyanobacteria that do inhabit the sub-crust soils account for most of the transcriptional activity. These data highlight the need to look at longer term pulsed experiments to investigate if and when other members of the below-crust soils become active.

Response of Arid Soil Communities to a Wetting Pulse

In general, mRNA molecules are considered to be short-lived and of modest abundance in the cell. For example, physiological studies of Escherichia coli suggest that actively growing cells contain approximately 1,000 mRNA molecules with a half-life of about 5 min (Bernstein et al., 2002). The mRNA content of bacterial cells in the open ocean is more modest, with approximately 200 mRNA transcripts per cell, also with half-lives in the range of a few minutes (Moran et al., 2013). However, transcriptional activity in active microbial cultures or other ecosystems may not be an appropriate comparison to biocrusts. Arid lands are characterized by extended dry periods, during which biocrust respiration is virtually undetectable suggesting that biocrusts are largely dormant during this time (Garcia-Pichel and Belnap, 1996; Huxman et al., 2004; Sponseller, 2007). Periods of metabolic activity, stimulated primarily by water availability, occur infrequently and are short in duration. In this regard, biocrust response to wetting may be more similar to the germination of plant seeds. Plant seeds are dormant (or quiescent) until the initiation of germination, which is triggered by the uptake of water (Nonogaki et al., 2010). The resumption of metabolic activity for plant seeds (or arid soil microbes) is dependent on mRNA produced by the organism prior to the period of dormancy. In seeds mRNA transcripts are protected as RNA/protein complexes, also referred to as mRNPs (Bewley, 1997; Reichel et al., 2016). A similar mechanism could serve to protect transcripts in dry biocrusts, allowing them to be stable over the duration of dry periods. Several cyanobacteria encode functional homologs to eukaryotic RNA-binding proteins (Mulligan et al., 1994). The specific functions of these genes have not been conclusively demonstrated, but they appear to play a role in cold-stress or nitrogen starvation (Sugita et al., 1999; Mori et al., 2003). Here, we propose another potential role for cyanobacteria mRNPs; that is protecting mRNA molecules from desiccation damage during biocrust drying. The observation of complex pools of mRNA transcripts in dry conditions across the soil types, including transcripts more abundant in the dry samples (Figures ), points to mechanisms for maintaining RNA integrity in dry arid soils. The rapid turnover of mRNA is also supported by the differences in the metatranscriptomes after only 1 h of wetting, particularly for the biocrusts (Figure ). This matches other observations that short-term transcriptional responses may be overlooked in linking transcriptional data to ecosystem process rates (Albright et al., 2018). Future work investigating the lifespan of arid soil mRNA molecules may shed light on transcript stability as a regulator of community responses to desiccation and rewetting.

A common observation across the dry samples was an enrichment of transcripts for hydrogenases (Figures ), specifically uptake hydrogenases (see ). Uptake hydrogenases are bi-directional, capable of both producing and consuming H2. Uptake hydrogenases have been found among all heterocystous cyanobacteria studied to date, as the enzyme captures H2 produced during nitrogen-fixation, although non-nitrogen fixing cyanobacteria also contain uptake hydrogenases (Tamagnini et al., 2005). Several cyanobacteria within the order Oscillatoriales have been found to be prodigious producers of H2, although H2 evolution was absent from M. vaginatus (Kothari et al., 2012). Thus, these transcriptomic data point to cyanobacteria possibly producing and/or consuming H2 in dry conditions. However, the majority of the uptake hydrogenase transcripts enriched in the dry soils were found from various Actinobacteria, specifically the Solirubrobacterales (see ). Scavenging of atmospheric hydrogen has been proposed as a mechanism by which bacteria can maintain metabolic activity under low carbon availability (Greening et al., 2015). Thus an active H2 based metabolism appears to be operating in the dry biocrusts that may involve a give and take between the Cyanobacteria and Actinobacteria.

By the metric of the number of differentially abundant transcripts, the wetting pulse was most pronounced in the biocrust soils (Figure ). Many of the subsystems to which the contigs were annotated were within central metabolic pathways and machinery, such a translation elongation factors, ATP synthase, ribosomal proteins, and protein chaperones (Figure ).

In the biocrusts most of the differentially abundant subsystems were made up of contigs with members that were identified as both more abundant in the wet and dry conditions. For example, contigs annotated to bacterial small subunit (SSU) ribosomal proteins (Figure ). It is well established that bacteria scale their ribosomal counts with metabolic activity (Elser et al., 2003; Steven et al., 2017). In the dry crusts the enriched contigs for SSU ribosomal proteins were primarily from groups of Actinobacteria whereas SSU genes enriched in the wet crusts were primarily cyanobacterial in origin (Supplementary Table ). Thus, these data point to different populations being active in the wet and dry conditions. Specifically, this is consistent with cyanobacteria becoming more active in the wet biocrusts. Concurrently, an enrichment of photosystems I, II, and phycobillisome genes supports a rapid recovery of the cyanobacteria and the initiation of photosynthesis (Figure ).

Metatranscriptomes and Biocrust Biogeochemical Cycling

A large driver for studying the microbial ecology of soils is to include soil microorganisms in ecosystem, carbon, and climate models (Reid, 2012). In this study, we focused on carbon and nitrogen cycling, given the large role biocrusts play in these processes in arid ecosystems (Barger et al., 2016; Sancho et al., 2016). Previous studies at this field site have demonstrated that both total organic carbon and total nitrogen are reduced under chronic trampling (Kuske et al., 2012). Furthermore, simulated rainfall events have detected that losses of dissolved organic carbon, dissolved organic nitrogen, and ammonium were significantly higher from trampled plots in comparison to undisturbed plots (Barger et al., 2006). However, these measurements have been scaled to investigate landscape level differences, not the environmental fluxes that are likely to occur with a short-term wetting pulse. Our data clearly demonstrates that significant shifts in gene expression occur within an hour of a wetting pulse, particularly for the control biocrust soils (Figure ). Studies have documented biogeochemical changes can be observed in these brief time periods. For example, the consumption of CO2 and production of chlorophyll are observed in minutes after wetting (Rajeev et al., 2013; Tucker et al., 2017). Future studies that are able to connect biogeochemical processes happening at the timescales relevant to transcriptomics may be the key to linking molecular sequence data to relevant ecosystem process rates.

Linking metatranscriptomes to biogeochemical cycling requires that transcriptional abundances correlate with ecological process rates, a relationship that has not been well established (Moran et al., 2013; Prosser, 2015). The data in these transcriptomes highlight some of the potential challenges in linking transcripts to flux rates. First, multiple transcripts can map to the same functions (Table ) and transcripts that map to the same function do not always respond in concert (Figure ). Multiple transcripts for the same function may arise from functionally redundant species encoding the same pathways, as well as the potential for a single genome to encode multiple gene copies that may be differentially regulated under different environmental conditions (Misra and Tuli, 2000; Allison and Martiny, 2008). Thus, we observe for example, that transcripts for functions such as denitrification and ammonification are identified as both more abundant in wet and dry conditions (Figure ). In this regard, assaying transcriptional abundance of functional genes, and linking them to ecological processes, will require accounting for the response of diverse species and differential gene expression. Second, a significant response to the wetting was that particular bacteria increased their transcripts for ribosomal proteins, presumably leading to elevated cellular ribosomal content (Figure ). Because mRNA molecules are not consumed during translation, and in bacteria multiple ribosomes can translate the same mRNA molecule in polysomes, an increased ribosome content can lead to more proteins being produced from the same number of mRNA molecules (Scott et al., 2010). In this manner, cells could increase the stock of biogeochemical cycling proteins through increased ribosomal content without a concurrent increase in the mRNA encoding those genes. In this regard, linking transcript abundance to ecosystem process rates will also need to account for factors such as ribosome content, RNA degradation, and translation efficiency (Pedersen, 1984; Kerkhof and Ward, 1993; Taylor et al., 2013).

Conclusion

The data presented here suggest that metatranscriptome sequencing can identify differences in biocrust communities in response to perturbations that occur over dramatically different time scales. Chronic physical trampling carried out over 17 years, alters the metabolic potential of the community as well as shaping the response of the community to a wetting pulse. In contrast, a short durational wetting pulse of 1 h is also sufficient to alter the metatranscriptome. The biocrust communities respond rapidly to a pulse of water, with the cyanobacteria ramping up protein synthesis and photosynthesis genes. Taken together, these observations support that chronic trampling induces shifts in the structure and function of the soil communities that could not be assessed at a visual scale alone, which is often used to asses biocrust health and function, thus highlighting the importance of molecular microbial ecology tools in characterizing biocrusts and the response and resiliency to ecosystem disturbances.

Author Contributions

BS, JB, and CK contributed to research design and writing of the manuscript. BS was primarily responsible for bioinformatic analyses.

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.

Table 1

Biogeochemical cycling genes assessed in biocrusts.

Protein name	SEED annotation	Biogeochemical role	Number of contigs
Nitrogen cycle
Nitrogenase	Nitrogenase (molybdenum iron) alpha chain (EC:1.18.6.1)	Nitrogen fixation	2
Nitrate reductase	Assimilatory nitrate reductase large subunit (EC:1.7.99.4)	Ammonification	24
Nitrite reductase	Nitrite reductase [NAD(P)H] large subunit (EC:1.7.1.4) (assimilatory)	Ammonification	20
Cytochrome c nitrite reductase	Cytochrome c nitrite reductase small subunit NrfH	Ammonification	1
Nitrite reductase	Copper containing nitrite reductase (EC:1.7.2.1)	Denitrification	15
Nitric oxide	Nitric oxide reductase (EC:1.7.99.7) quinol dependent	Denitrification	5
Nitrous oxide	Nitrous oxide reductase (EC:1.7.99.6)	Denitrification	5
Carbon cycle
RuBisCo	Ribulose bisphosphate carboxylase large chain (EC:4.1.1.39)	Carbon fixation	96
Malate dehydrogenase	Malate dehydrogenase (EC:1.1.1.37)	TCA cycle	42
Pyruvate kinase	Pyruvate kinase (EC:2.7.1.40)	Entner–Doudoroff pathway	43
Lactate dehydrogenase	D lactate dehydrogenase (EC:1.1.1.28)	Fermentations	14
Methane monooxygenase	Methane monooxygenase component A beta chain (EC:1.14.13.25)	Methanotrophy	13