Transcriptional Landscape and Regulatory Roles of Small Noncoding RNAs in the Oxidative Stress Response of the Haloarchaeon Haloferax volcanii

ABSTRACT Haloarchaea in their natural environment are exposed to hypersalinity, intense solar radiation, and desiccation, all of which generate high levels of oxidative stress. Previous work has shown that haloarchaea are an order of magnitude more resistant to oxidative stress than most mesophilic organisms. Despite this resistance, the pathways haloarchaea use to respond to oxidative stress damage are similar to those of nonresistant organisms, suggesting that regulatory processes might be key to their robustness. Recently, small regulatory noncoding RNAs (sRNAs) were discovered in Archaea under a variety of environmental conditions. We report here the transcriptional landscape and functional roles of sRNAs in the regulation of the oxidative stress response of the model haloarchaeon Haloferax volcanii. Thousands of sRNAs, both intergenic and antisense, were discovered using strand-specific sRNA sequencing (sRNA-seq), comprising 25 to 30% of the total transcriptome under no-challenge and oxidative stress conditions, respectively. We identified hundreds of differentially expressed sRNAs in response to hydrogen peroxide-induced oxidative stress in H. volcanii. The targets of a group of antisense sRNAs decreased in expression when these sRNAs were upregulated, suggesting that sRNAs are potentially playing a negative regulatory role on mRNA targets at the transcript level. Target enrichment of these antisense sRNAs included mRNAs involved in transposon mobility, chemotaxis signaling, peptidase activity, and transcription factors. IMPORTANCE While a substantial body of experimental work has been done to uncover the functions of small regulatory noncoding RNAs (sRNAs) in gene regulation in Bacteria and Eukarya, the functional roles of sRNAs in Archaea are still poorly understood. This study is the first to establish the regulatory effects of sRNAs on mRNAs during the oxidative stress response in the haloarchaeon Haloferax volcanii. Our work demonstrates that common principles for the response to a major cellular stress exist across the 3 domains of life while uncovering pathways that might be specific to the Archaea. This work also underscores the relevance of sRNAs in adaptation to extreme environmental conditions.

Archaea. Haloarchaea have previously been shown to be highly resistant to ROS damage, withstanding many times what Escherichia coli and other radiation-sensitive organisms can survive (5)(6)(7). The haloarchaeon Halobacterium salinarum has been shown to use a wide range of strategies to combat damage from oxidative stress, including multiple copies of genomes (polyploidy) as the substrate for DNA repair, functional redundancy of DNA repair and detoxification enzymes (e.g., catalase), increased cytosolic manganese complexes to scavenge ROS, and differential regulation of genes in response to stress (5)(6)(7)(8)(9). However, pathways for DNA repair and protein turnover in haloarchaea are nearly identical to those in nonresistant bacteria and eukarya, suggesting that the regulation of these processes in response to oxidative stress might be key to their robustness. Previous work with H. salinarum oxidative stress gene regulatory networks revealed that a single transcription factor, RosR, regulates the appropriate dynamic response of nearly 300 genes to reactive oxygen species stress (5). This work demonstrated that the oxidative stress response in H. salinarum impacted a wide array of cellular processes, engaging at least 50% of all the genes (2). These results underline the importance of gene regulation in haloarchaea for responding to and counteracting the damage caused by oxidative stress.
In addition to transcription factors, small regulatory RNAs (sRNAs) similarly act as global gene regulators (10). Small RNAs (sRNAs) are ubiquitously found in Bacteria and Eukarya, playing large-scale roles in gene regulation, transposable element silencing, defense against disease state, and foreign elements (11)(12)(13)(14). Several types of sRNAs have been identified in the Eukarya (microRNAs [miRNAs], small interfering RNAs [siRNAs], and Piwi-interacting RNAs [piRNAs]) and they are typically 20 to 25 nucleotides (nt) long. Their major mode of interaction is through base pairing to the 3= untranslated region (UTR) of their target mRNAs, inhibiting translation or triggering target degradation with associated protein components (Argonautes) (10). Bacterial sRNAs have been shown to modulate core metabolic functions and stress-related responses, such as nutrient deprivation, by binding target mRNAs and causing their degradation or preventing translation (11,15). Most of the functionally characterized sRNAs in Bacteria bind the 5= UTR of their target mRNA and are longer than their eukaryal counterparts, with sizes ranging from 50 to 500 nt. These sRNAs can target multiple genes, including key transcription factors and regulators (11,15,16). As a consequence, a single sRNA can modulate the expression of large regulons and thus have a significant effect on metabolic processes. For example, the bacterial sRNA OxyS, which is dramatically induced by oxidative stress, regulates the expression of approximately 40 genes and interacts directly with eight target mRNAs (11).
sRNAs have been discovered to be abundant in Archaea, more specifically in haloarchaea, in response to a variety of environmental conditions, but the functional roles of these RNAs still remain poorly understood and a link between sRNAs expression and oxidative stress response has not been established (13,(17)(18)(19)(20)(21)(22)(23)(24). Only a few studies on sRNAs in hyperthermophiles, methanogens, and the haloarchaeon Haloferax volcanii have been reported so far (13,(17)(18)(19)(20)(21)(22)(23)(24). In H. volcanii, large numbers of intergenic and antisense-encoded sRNAs, 145 and 45, respectively, were discovered by using microarray in addition to a novel class of sRNAs recently described in eukaryotes, tRNA-derived fragments (tRFs), and a new study found thousands of sRNAs present in this organism (19,24,25). In Sulfolobus solfataricus, 125 trans-encoded sRNAs and 185 cis-antisense sRNAs were identified using high-throughput sequencing (HTS), suggesting that 6.1% of all genes in S. solfataricus are associated with sRNAs (26). A comparative genome analysis of Methanosarcina mazei, M. barkeri, and M. acetivorans revealed that 30% of the antisense and 21% of the intergenic sRNAs identified were conserved across the 3 species (27). Coimmunoprecipitation with the Lsm protein (archaeal Hfq homolog) was used to "capture" sRNAs (17), but its functional role remains to be elucidated. While Ago homologs have also been found in archaeal genomes, there is no evidence for eukaryote-like RNA interference in these organisms (28). Rather, a defensive role against foreign genetic material was recently proposed, whereby archaeal Ago proteins direct guide-dependent cleavage of foreign DNA (29)(30)(31). Altogether, these studies suggest that sRNAs are as widespread and abundant in the Archaea as in the Bacteria and Eukarya.
Target mRNA identification of sRNAs has proven to be difficult within the Archaea but is a necessary task for uncovering sRNA functionality. Transcriptome sequencing (RNA-seq) in M. mazei cultures, grown under nitrogen starvation conditions, showed the differential expression of a number of sRNAs in response to nitrogen availability and enabled the identification of the first in vivo target for archaeal intergenic sRNAs (27,32). The potential target for sRNA 162 is a bicistronic mRNA encoding a transcription factor involved in regulating the switch between carbon sources and a protein of unknown function (32). In Pyrobaculum, 3 antisense sRNAs were found opposite a ferric uptake regulator, a triose-phosphate isomerase, and transcription factor B, supporting a potential role for archaeal antisense sRNA in the regulation of iron, transcription, and core metabolism (33). sRNA deletion mutants can be used to identify potential biological functions and target genes. Deletion strains were successfully generated for H. volcanii, and phenotyping of the sRNAs deletion mutants revealed several severe growth defects under high temperatures, low salt concentrations, or specific carbon sources (9,22). While these studies revealed that sRNAs likely play essential roles in the physiological responses to environmental challenges in the Archaea, the functional roles and mechanisms of action of these important posttranscriptional regulators still remain unknown. Furthermore, no work has been done to investigate archaeal sRNAs in response to oxidative stress, a universal and frequent stressor in all domains of life that results in extensive cellular damage. To determine the impact of sRNAs during the oxidative stress response, we assessed the H. volcanii transcriptional landscape under no-challenge and oxidative stress conditions using comparative strand-specific small RNA sequencing (sRNA-seq).

RESULTS
To identify globally small noncoding RNAs differentially expressed in response to oxidative stress in H. volcanii, we exposed 5 replicate cultures of H. volcanii to 2 mM H 2 O 2 , a dose that resulted in the survival of 80% of the cells (see Fig. S1 in the supplemental material). RNA from these H 2 O 2 -treated cultures and from no-challenge cultures (controls) was sequenced using a strand-specific size-selected (50 to 500 nt) sRNA library preparation essential for sRNA discovery.
Small noncoding RNA discovery in H. volcanii and normalized expression values. We obtained at total of 137 million sequence reads (41 Gb) across all replicates and conditions. Following quality control and reference-based read mapping, we intersected the mapped reads against the H. volcanii reference genome to discover sRNA transcripts that we classified as antisense (overlapping a gene and/or its regulatory elements on the opposite strand) (Fig. 1A) and intergenic (in a noncoding region between two genes) (Fig. 1B). We were unable to identify previously described cisinternal sRNAs (24), because with our RNA-seq strategies, those transcripts would be confounded with the corresponding gene transcripts. To further validate the sRNAs we identified and reduce transcriptional noise in our data, we applied a rigorous twopronged in silico approach. First, we used a stringent threshold requiring an sRNA (i) to be present in at least four of five biological replicate libraries, and (ii) to have a minimum expression of 40 transcripts per million ([TPM] average with standard deviation from replicates) for antisense sRNAs (asRNAs) and 14 TPM for intergenic sRNAs. Second, we used the genome browser IGV to inspect visually and confirm each sRNA. These novel transcripts represented 25 to 30% of the total transcriptome and were 50 to 1,000 nt in length ( Fig. 2; see also Table S1). Putative mRNA targets for asRNAs were identified as the cis mRNAs encoded on the opposite strands with a minimum overlap of 8 nt, on the basis of IGV confirmation and reports for bacterial asRNAs (34). However, we found in our data that the median overlap between sRNAs and putative cis-mRNA targets was 221 nt, with only a small number having a minimum overlap length of 8 nt (see Fig. S2). Analyzing the upstream regions of sRNAs enabled the discovery that 30% of sRNAs contained a BRE and TATA box with centroids at nt Ϫ38 and Ϫ29 (see Fig. S3).
The use of less-conservative parameters (nt Ϫ3 and ϩ3) for BRE and TATA box centroids resulted in 70% of sRNAs having transcriptional motifs.
Normalized expression values in RNA-seq analyses are often reported as reads or fragments per kilobase of transcript per million mapped reads (RPKM or FPKM). However, RPKM/FPKM has been shown to be inconsistent for comparisons between samples due to variable transcript lengths. Another expression value, transcripts per million (TPM), was found to be preferable, because it is independent of the mean expressed transcript length and TPM normalization performs better in multiple library comparisons (35)(36)(37)(38). To minimize transcript length bias (sRNAs are generally smaller in length than protein-encoding genes) and to compare sRNA and mRNA expression levels, we chose to use TPM in our analysis.

Noncoding sRNA characterization in H. volcanii under no-challenge conditions.
H. volcanii grown under no-challenge conditions (42°C, complex medium) expressed a total of 1,533 sRNAs after quality control (as described above) ( Fig. 2; Table S1), ranging from 49 to 1,000 nt in size and with an average length of 373 nt. A majority of these sRNAs, 1,480 sRNAs (97%), were antisense to coding regions (Fig. 2). The H. volcanii H53 auxotroph genome is 4 Mbp and contains 4,130 genes. The genome constitutes a chromosome stably integrated with small chromosome pHV4 and 2 small chromosomes (pHV1 and pHV3), and it has been cured of plasmid pHV2. A majority of sRNAs (68%) were carried on the chromosome and integrated small chromosome pHV4 (18%). No sRNAs were found on plasmid pHV2, as expected, while sRNAs were carried on the remaining small chromosomes pHV1 (2%) and pHV3 (12%).
We compared mRNA expression to sRNA expression by constructing mRNA sequencing (mRNA-seq) libraries, using the same RNA pool and library preparation as the sRNA-seq libraries (omitting size selection) and calculating transcript expression as TPM (see Fig. S4A). Relative to mRNA expression levels (average, 312.1 Ϯ 1,079.1 TPM), the expression of sRNAs was on average higher (average, 1,107.9 Ϯ 137.6 TPM). We also found that a majority of sRNAs had higher expression levels than mRNAs (sRNA: range, 14.1 to 905,191.0 TPM; median, 108.5 TPM; mRNA: range, 1.0 to 210,162.0 TPM; median, 15.1 TPM). Overall, 75% of sRNAs had expression values less than or equal to 320 TPM, 15% had expression levels similar to those of highly expressed mRNAs, (500 to 10,000 TPM), and 16 sRNAs (sRNAs no. 1771 to 1786) (Table S1) had very robust expression levels with TPMs ranging from 10,000 to 60,000 TPM (Fig. S4A). Lastly, 1 intergenic sRNA (STRG.2577.4), 111 nt in size, exhibited expression levels higher than any mRNA with a TPM of 905,191 (Table S1). Transcript length did not correlate with expression levels, indicating that when we observed sRNAs with low expression levels, it was not an artifact of sequencing (i.e., longer transcripts receiving more read coverage thus skewing coverage based on length) (see Fig. S5). We found that 4 of the 5 most highly expressed sRNAs (TPM Ͼ 30,000) were antisense to coding regions (Table S1).
Regulatory effects and differential expression of sRNAs during oxidative stress. To investigate the regulatory effects sRNAs on their target mRNAs, we compared the expression levels (TPM) of all antisense sRNAs against the expression levels (TPM) of the in silico predicted putative cis-mRNA targets (Fig. 3). We found that a large population of asRNA-mRNA cis pairs exhibited lower mRNA target expression than the sRNA, under both experimental conditions ( Fig. 3; Table S1). We conducted a pairwise t test between these cis pairs on expression level differences between sRNAs and putative cis-target mRNAs and found that 755 of sRNAs had significantly (P Ͻ 0.05) higher expression, potentially indicating a negative regulatory relationship between sRNAs and putative cis-mRNA targets (Fig. 3).
To further investigate this negative regulatory relationship between sRNAs and putative mRNA targets, we probed for differentially expressed sRNAs between the no-challenge and the oxidative stress conditions. Candidate sRNAs were considered FIG 3 Transcript per million (TPM) expression levels between sRNAs and their putative cis-mRNA targets during oxidative stress. Each point represents the TPM ratio between an sRNA and its putative cis-mRNA target. TPM values are the averages from sRNA and mRNA replicates with error bars representing standard deviations among replicates. A pairwise t test was conducted between sRNA and putative cis-mRNA target replicates to infer significant difference in TPM expression between the cis pairs. Orange points indicate a P value Ͻ0.05, and gray points indicate a P value Ͼ0.05. The red line represents no change in the ratio of sRNA-cis-mRNA expression (slope ϭ 1).
significantly up-or downregulated by oxidative stress using a false-discovery rate (FDR) of less than 5%. Using this statistical framework, we identified a core set of differentially expressed sRNAs specific to oxidative stress (Fig. 4A). Both intergenic and antisense sRNAs were differentially expressed. Of the intergenic sRNAs, 48 were significantly differentially expressed (FDR Ͻ 0.05), with 23 upregulated and 25 downregulated (see Fig. S6; Table S2). Of these upregulated intergenic sRNAs, 79% had an increase in expression greater than or equal to a fold change of 2 (log 2 fold change ϭ 1) during oxidative stress, with the most upregulated intergenic sRNA (STRG.277.2) having a 16-fold increase. On the other hand, a majority of downregulated intergenic sRNAs had large fold changes in expression (Յ2-fold change), and 5 exhibited very robust downregulation (ՅϪ4-fold change). A total of 605 antisense sRNAs were significantly (FDR Ͻ 0.05) differentially expressed. These sRNAs were either upregulated (309) or downregulated (296) during oxidative stress, indicating two populations of antisense sRNAs ( Fig.  4A; Table S2). Fifty percent (302 sRNAs) of these differentially expressed sRNAs demonstrated a fold change in expression of Ϯ2 or greater; the most upregulated sRNA had a fold change of 15 and the most downregulated sRNA had a fold change of Ϫ9, Journal of Bacteriology indicating a role for these sRNAs in the cellular response to oxidative stress. More antisense sRNAs were upregulated with a fold change in expression of 4 or greater (28 sRNAs) than were downregulated (20 sRNAs). We then compared differential expression levels between all significantly upregulated antisense sRNAs and their putative cis-mRNA targets and found that a population (133 sRNAs [22%]) of these upregulated antisense sRNAs had putative mRNA targets that were downregulated during oxidative stress (Fig. 4B). For example, during oxidative stress, 36 upregulated antisense sRNAs targeted transposase mRNAs, and 24 of these antisense sRNAs had putative cis-target transposase mRNAs downregulated (Fig. 4D). While a smaller subset of significantly (FDR Ͻ 0.05) downregulated antisense sRNAs had their putative cis-mRNA target upregulated during oxidative stress (72 [11%]), it also indicated a correlative potential negative regulatory effect (Fig. 4C).
Oxidative stress-responsive asRNAs were predicted to overlap both the 5= and 3= UTRs of mRNAs (Fig. 5). We found that 7% of asRNAs overlapped at the 5= UTRs and 26% overlapped at the 3= UTRs. However, the majority of the asRNAs (67%) overlapped the coding sequences (CDSs) of mRNAs rather than the UTRs, which has not been previously reported (Fig. 5). We calculated that, on average, the overlap between sRNAs and their putative cis-target mRNAs was 265 nt, the range was 8 to 992 nt, and the peak overlap was between 150 and 200 nt (Fig. S2). Using Northern blot analysis, we recapitulated the in vivo differential expression patterns of selected candidate sRNAs, further confirming the transcript sizes and differential expression levels for oxidative stress sRNAs (STRG.3823.1, -4700.1, -8.6, -277.2, -3733.1, -2983.4, and -4213.4) ( Fig. 6A and B). We also showed that the strandedness (the strand on which the sRNA was carried) predicted by our sRNA-seq analysis was confirmed by our in vivo data using oligonucleotide probe Northern blotting of 5= UTRs, 3= UTRs, CDS antisense sRNAs, and intergenic sRNAs (Fig. 6B).
Target enrichment of sRNAs. We identified in silico targets for the differentially expressed oxidative-stress responsive antisense sRNAs. Genes for the putative target mRNAs were categorized according to cellular function by using archaeal Cluster of Orthologous Genes (arCOGs) and according to pathways by using gene ontologies (GO) from Database for Annotation, Visualization, and Integrated Discovery (DAVID). For sRNAs upregulated during H 2 O 2 stress, we found a functional enrichment of target genes encoding transposases, involved in chemotaxis methyl-receptor signaling and in transcriptional regulation (transcription factors) (P Ͻ 0.05). Genes from many other pathways that were not enriched were also the targets of antisense sRNAs, including peptidase activity genes and serine and threonine biosynthesis genes (Fig. 7A). Twentythree of these sRNAs targeted transposase genes. Each of these transposase mRNAs was downregulated, while their cognate sRNA was upregulated, and the sRNA was always located at the 5= UTR of its target. Most transposases belonged to the IS family of transposases, except for one DDE transposase. Three transcription factor families (IclR, ArcR, and AsnC) were also targeted by asRNAs. A functional enrichment gene ontology analysis found that downregulated sRNAs targeted genes involved in membrane transport (ABC transporters) and the biosynthesis of secondary metabolites, as well as hydrolases (Fig. 7B). A significant proportion of enriched targets for both upand downregulated sRNAs were genes encoding hypothetical proteins.
mRNA transcriptional response to oxidative stress in H. volcanii. To determine the transcriptional landscape of mRNAs during oxidative stress, especially for mRNAs that were predicted targets of sRNAs, we sequenced rRNA-depleted mRNA-seq libraries in parallel with the previously described sRNA-seq libraries (derived from the same pool of total RNA). During H 2 O 2 -induced oxidative stress, one-fourth of all genes (1,176) were significantly differentially expressed with a false-discovery rate of less than 5% (see Table S4). Both catalase and superoxide dismutase, known ROS detoxification enzymes, were upregulated at the mRNA level, thus validating our experimental approach and

DISCUSSION
Previous studies of sRNAs in Archaea revealed the abundance of sRNAs within the third domain of life and have been pivotal in establishing a working hypothesis on archaeal sRNA functionality. These studies have been limited to (i) microarray studies that do not enable de novo discovery of sRNAs, (ii) differential RNA-seq (dRNA-seq) approaches, which select only for primary transcripts and do not provide length (nt) information nor expression information (only coverage), and (iii) individual sRNA studies, which do not give a holistic view of the pathways being regulated within the cell. Using a custom strand-specific sRNA-seq library preparation and analysis pipeline, we have developed a method to perform high-throughput analysis of the sRNA transcriptional landscape, expression, and regulatory effects, and to identify regulated gene pathways in response to environmental stressors within the Archaea. In this study, we propose that sRNA-mediated transcriptional regulation is key in regulating stress responses to environmental challenges, such as oxidative stress, in the haloarchaea. sRNAs have the potential to carry out large-scale regulation of genes involved in the oxidative stress response, resulting in increased resistance to extreme environmental stressors.
The discovery that sRNAs comprised nearly one-third of the total transcriptome of H. volcanii and included basal transcriptional promoters, under both no-challenge and oxidative stress conditions, suggests that sRNAs have an important functional role under a variety of environmental conditions. We discovered thousands of sRNAs expressed in H. volcanii, with the majority being antisense to genes, indicating that antisense transcription was ubiquitous within the cell. This is in stark contrast to most of the literature reporting that a majority of sRNAs discovered in Archaea were intergenic (9,13,22,39). This discrepancy is likely due to previous studies being limited to microarray approaches. Indeed, a recent study using directional RNA-seq (dRNA-seq) to map all transcription start sites (TSS) in H. volcanii found thousands of novel TSS, with 1,244 of these TSS being antisense to mRNAs (24). Most of the TSS (75%) of the sRNAs we discovered in H. volcanii had the same TSS (Ϯ10 nt) as those found in the dRNA-seq study by Babski et al. (24), further confirming our results. This underlines the importance of HTS studies, especially strand-specific RNA-seq such as in our study, to discover the full extent of antisense sRNA expression in Archaea.
Our finding suggests that a priori cis-acting sRNAs may play a larger role than trans-acting sRNAs within the cell. It should not be overlooked that the difficulty in finding in silico targets for intergenic sRNAs, because these sRNAs do not form 100% complementarity with their targets, might suggest that they have multiple mRNA targets. Antisense and intergenic sRNAs are broad classifications used in the archaeal small noncoding RNA field, but our data revealed that further classifications can be made on the basis of sRNA-mRNA binding characteristics (5= UTR, 3= UTR, and CDS), differential expression, and regulatory effects. We found that only a small fraction of asRNAs targeted the 5= UTRs of mRNAs, which is in concurrence with work demonstrating that most mRNAs in H. volcanii are leaderless (lacking a 5= UTR). A majority of the 5= UTR-binding sRNAs targeted transposons, providing further evidence that they may constitute their own class of sRNAs. Within this context, 3= UTR-binding sRNAs should also be considered another class of sRNAs, resembling eukaryotic sRNAs, which likely lead to the degradation of their target transcripts. The majority of the asRNAs we identified in H. volcanii had 100% complementarity within the CDS of their target mRNAs. This is the first report of such a finding in any domain of life and might constitute an attribute unique to archaeal sRNAs. We could not identify any "seed" binding region for these CDS-binding sRNAs, indicating that they likely have full occupancy on the mRNA. It is worth noting that there were only a few instances (Ͻ20) where sRNAs overlapped more than the full length of the target mRNA (i.e., the asRNA is longer than the mRNA) or overlapped multiple cis-mRNA targets (i.e., the asRNA overlaps the 3= UTR of one mRNA and the 5= UTR of an adjacent mRNA), which is novel in the Archaea.
Most H. volcanii sRNAs had a normalized expression value of 200 TPM or less, indicating that sRNA transcripts are relatively abundant in the cell. The use of a stringent thresholding approach resulted in a smaller number of more highly expressed sRNAs but avoided potentially reporting false positives with low TPM values or transcriptional noise. In comparison, most mRNAs within H. volcanii had TPM values of 20 or less. We found a population of antisense sRNAs with significantly higher expression than their putative cis-mRNA targets, suggesting a potential negative regulatory role in sRNA-mRNA interactions ( Fig. 3 and 4B) (40)(41)(42). This trend extended to many asRNAs (under both no-challenge and oxidative stress conditions), down to sRNAs with expression levels of 40 TPM. Stronger evidence for a negative regulatory effect lies with a population of upregulated sRNAs. A group of significantly upregulated asRNAs had target mRNAs that were downregulated during oxidative stress, indicating these asRNAs could potentially negatively regulate mRNA targets at the transcript level. Whether this negative regulation is occurring during transcription initiation/elongation or if these asRNAs are causing mRNA degradation is currently unknown. All the asRNAs targeting transposons at the 5= UTR were upregulated, and the transposon mRNAs were downregulated (Fig. 4D), suggesting that these asRNAs might have a similar mechanistic function (6,43,44). If indeed these asRNAs are negatively regulating their target mRNAs in H. volcanii, we expected to find that downregulated asRNAs have upregulated target mRNAs. Some downregulated asRNA targets exhibited this pattern, further supporting a potential negative regulation (Fig. 4C). Many mRNA targets of both up-and downregulated asRNAs exhibited the same regulatory pattern, suggesting a potential nonnegative regulatory role for a group of asRNAs (Fig. 4). Alternative hypotheses, reflecting the complexity of transcriptional regulation in the Archaea, can be formed: (i) these asRNAs may have a positive regulatory effect, such as stabilizing target mRNAs and masking them from degradation as seen with an intergenic sRNA (itsRNA) in M. mazei (23), (ii) trans-acting intergenic sRNAs might be targeting these mRNAs, negatively regulating them, and (iii) some may have an unknown function (23).
The most enriched negatively regulated sRNA targets were transposases, chemotaxis proteins, and transcription factors. It has been shown that transposons are opportunistic, and under stress conditions, they can wreak havoc by hopping around in the genome causing double-strand breaks and hence need to be silenced (6,43,44). A functional enrichment of IS4 transposon genes that are downregulated during oxidative stress supports our observation that the upregulated sRNAs can potentially negatively regulate transposons and suggests that transposon activity is tightly regulated during oxidative stress in H. volcanii. sRNA-mediated regulation of chemotaxis transducer proteins during oxidative stress suggests interesting implications in sensing ROS and motility. H. volcanii expresses genes coding for a flagellum analog named "archaellum," which are organized into an operon and are regulated by a network of regulators called the archaellum regulatory network (arn) (identified in crenarchaea) (45,46). The regulation of these motility genes is still under investigation and so far is restricted to a few examples, such as under H 2 or nitrogen limitation conditions in Methanococcus jannaschii and Methanococcus maripaludis (46)(47)(48)(49). No direct transcriptional regulators of the archaellum have been identified in Euryarchaeota, but the deletion of archaellin genes, the presence of the H-domain set of type IV pillins, and Agl proteins have been shown to affect the assembly of the archaellum in H. volcanii (46,(50)(51)(52). Integral to how microorganisms maintain homeostasis in stressful and fluctuating environments are the gene regulatory networks composed of interacting regulatory transcription factors and their target gene promoters (53). Our discovery that sRNAs are targeting transcription factors provides evidence that sRNAs are likely deeply interlaced within complex gene regulatory networks of H. volcanii, and these sRNAs are key to maintaining homeostasis during environmental stresses such as oxidative stress. Many mRNA targets of differentially regulated sRNAs were hypothetical proteins, indicating that much remains to be elucidated in this organism.
Our whole transcriptional analysis demonstrated that close to one-third of the genes (ϳ1,100 [30%]) in H. volcanii were differentially regulated during constant H 2 O 2induced oxidative stress at ϳ80% survival, which is in agreement with the transcriptional response of H. salinarum during constant H 2 O 2 -induced (929 genes [38%]) and paraquat-induced (1,099 genes [45%]) oxidative stresses at ϳ80% survival (2). This indicates that transcriptional regulation is crucial in order to mount this oxidative stress response via gene activation and repression. Two single-stranded DNA binding proteins (RpaB and RpaC) were found necessary for the increased survival of H. volcanii under ionizing radiation (a proxy for desiccation) and UV radiation, stressors that both cause oxidative stress (54, 55) (J. DiRuggiero, unpublished data). In H. salinarum, Rpa operons were upregulated under ionizing radiation as well and contributed to resistance (56,57). In conjunction with previous findings, we observed that two of the most upregulated genes during H 2 O 2 oxidative stress were RpaB1 (HVO_RS10725, 23-fold change) and RpaB2 (HVO_RS06105, 8-fold change), confirming their roles in oxidative stress resistance in H. volcanii and likely other haloarchaea. One gene, a reactive intermediate/ imine deaminase RidA-homolog (HVO_RS12485), was upregulated orders of magnitude more than any other gene. The encoded protein is known to be involved in the synthesis of branched-chain amino acids by speeding up the IlvA-catalyzed deamination of threonine into 2-ketobutyrate (58,59). Previous work has shown that in the presence of reactive chlorine species (RCS), such as HOCl, imine deaminase seemed to inhibit IlvA activity suggesting that imine deaminase may have a different function in the presence of RCS (58,60). Further studies found that imine deaminase can sense RCS and, in doing so, becomes a chaperone that prevents protein aggregation (58). Reactive oxygen species in hypersaline environments produce RCS (61). In addition, ROS cause extensive irreversible protein damage such as carbonylation, which in turn causes protein aggregation (62,63). This reactive intermediate/imine deaminase is the most upregulated protein-encoding gene, suggesting that it may be playing a similar chaperone role to prevent protein aggregation, by sensing either ROS or RCS produced by H 2 O 2 (58,60).
This study is the first to report on the transcriptional response of H. volcanii to oxidative stress, and while we found similar responses to H 2 O 2 exposure as previously reported for H. salinarum (2), further validating our work and providing evidence that haloarchaea have evolved similar strategies to survive their environmental stresses, we also found responses that were unique to H. volcanii. Similarities to H. salinarum include the upregulation of genes encoding ROS scavenging proteins (catalase and superoxide dismutase) and iron sulfur assembly proteins (SufB and SufD), as well as proteasome genes, indicating high protein turnover, and many DNA repair genes (2). Most of the downregulated genes were involved with metabolism, such as sugar/phosphate/ peptide ABC transporters, electron carriers (halocyanin), and tricarboxylic acid (TCA) cycle enzymes, possibly to halt growth until damage is repaired (2,64). The most downregulated gene was for an iron ABC transporter, most likely to limit further production of ROS via Fenton reactions (2). Of the unique responses to oxidative stress in H. volcanii, we found that genes encoding all of the RNA polymerase subunits and transcription elongation factors and seven basal transcription initiation factors (HVO_RS12755, 18.4-fold change; HVO_RS01380, 8.6-fold; HVO_RS09745, 6.5-fold; HVO_RS01840, 2.1-fold; HVO_RS05475, 1.74-fold; HVO_RS11835, 1.5-fold; HVO_RS05475, 1.3-fold change) were significantly upregulated in response to oxidative stress. The increase in sRNAs during oxidative stress could be attributed to this increase in transcription machinery. The majority of the 30S and 50S ribosomal subunits were downregulated, in contrast to that in H. salinarum. The upregulations of genes coding for histidine biosynthesis and catabolism into glutamate and 2-oxocarboxylic acid metabolism were unknown to be involved in the oxidative stress response, which further demonstrates there are still more mechanisms to uncover for oxidative stress resistance. RosR was identified as a global transcriptional regulator in H. salinarum, and it was strongly upregulated during oxidative stress (5). RosR demonstrated no differential expression during oxidative stress in H. volcanii, indicating that it may play another role in this organism. Cell cycle genes (parA and cdc6) involved in chromosome segregation (65) were downregulated, further suggesting that division is being arrested (halting growth) in order to repair damage.
In this study, we showed for the first time that small noncoding RNAs are specifically associated with the oxidative stress response in archaea. During oxidative stress, asRNAs were prevalently transcribed, comprising nearly 30% of the transcriptome of H. volcanii, and many upregulated asRNAs imparted a correlative negative regulatory effect on putative target mRNAs. These results support the hypothesis that some asRNAs in Archaea behave similarly to cis-acting bacterial sRNAs and eukaryotic siRNAs, which negatively regulate mRNAs by sharing extensive complementarity and facilitating RNA degradation (66,67). The precise mechanism(s) of sRNA-mRNA mediated regulation, in particular, whether proteins are required to complex with sRNAs in order to mediate gene regulation, such as in Bacteria (Hfq) and Eukarya (Ago), remains to be elucidated. We also identified several classes of asRNAs, on the basis of their mRNAbinding patterns (3= UTR, 5= UTR, and CDS) and showed that CDS targeting of mRNAs was the predominant mode of action for sRNA hybridization. The mechanistic differences between these classes of sRNAs, as well as the regulatory roles of sRNAs in archaea and their functional importance in adaption to extreme environments, still need to be investigated. Oxidative stress exposure. We exposed 5 biological replicates of H. volcanii strain H53 liquid cultures to the oxidative stress agent H 2 O 2 . Initially, cultures were grown in 80 ml of Hv-YPC under optimal conditions to an optical density (OD) of 0.4 (mid-exponential phase). To ensure homogeneity, each replicate was subsequently split into two 40-ml cultures, one used for the no-challenge (control) condition and the other for the oxidative stress condition. For the latter condition, 2 mM H 2 O 2 (80% survival rate, previously determined) was directly added to the cultures followed by a 1-h incubation at 42°C with shaking at 220 rpm. Cultures were then rapidly cooled down and centrifuged at 5,000 ϫ g for 5 min, and the pellets were resuspended in 18% seawater. The cell suspensions were then transferred to 1-ml tubes and centrifuged at 6,000 ϫ g for 3 min, and the pellets were flash frozen and stored at Ϫ80°C until ready for RNA extraction. Control no-challenge culture replicates were processed in the same manner without the addition of H 2 O 2 .

MATERIALS AND METHODS
Oxidative stress survival curves. The assessment of survival in H. volcanii under oxidative stress conditions was conducted using microdilution plating as described previously (7). The counts from replicates were averaged and the standard deviations calculated. Survival was calculated as the number of viable cells following H 2 O 2 treatment divided by the number of viable untreated cells and graphed with standard error bars.
RNA extraction. Total RNA was extracted using the Zymo Quick-RNA Miniprep kit with the following modifications: after the addition of RNA lysis buffer to the frozen cell pellets, the cells were processed with 23-gauge needles and syringes to ensure complete cell lysis. H. volcanii liquid cultures are slimy and viscous; thus, to increase cellular lysis, 23-gauge needles and syringes were used to break down the cell pellets. Total RNA was then extracted according to the standard kit protocol.
Small RNA-sequencing library preparation. Total RNA for each biological replicate and condition was size selected by denaturing polyacrylamide gel electrophoresis. Twenty micrograms of total RNA was loaded onto a 7% denaturing urea polyacrylamide gel (SequaGel; National Diagnostics) in 0.5ϫ Trisborate-EDTA (TBE) buffer and run at a constant power of 30 W until the bromophenol blue bands reached the bottom of the gel. The gel was stained with SYBR Gold and visualized on a blue light box, and bands in the 50 to 500 nucleotide range, as indicated by the RNA Century Marker plus ladder (Thermo Fisher), were excised. Small RNAs (sRNA) were eluted by rotating overnight in 1.2 ml 0.3 M NaCl and were ethanol precipitated and DNase I (NEB) treated (37°C for 2 h) as previously described (69). Strand-specific libraries were prepared using the SMART-seq Ultralow RNA input kit (TaKaRa), and insert sizes were checked with the Bioanalyzer RNA pico kit (Agilent). Paired-end sequencing (2 ϫ 150 bp) was carried out on the Illumina HiSeq 2500 platform at the Johns Hopkins University Genetic Resources Core Facility (GRCF).
mRNA-sequencing library preparation. Individual mRNA-seq libraries were made from the same RNA pools as the sRNA-seq libraries above. Total RNA was rRNA depleted using the Ribo-zero bacteria kit (Illumina) and sequenced using the sRNA-seq library preparation method described above but omitting the size selection by denaturing urea polyacrylamide gel electrophoresis.
sRNA-and RNA-seq read quality control and reference-based read mapping. The assessment of the quality of each sequencing library read was conducted using fastqc. The program trim galore was used with base settings to trim adapter sequences from reads and to filter out low phred score reads (Ͻ20). Short-length reads were preserved. Reads from each replicate were aggregated together per condition to get a set of consensus sRNAs and were mapped against the H. volcanii NCBI RefSeq genome (taxonomy identification [taxid] 2246; 1 chromosome, 4 plasmids) using the hisat2 aligner with strandspecific options turned on and splice aware options turned off, in paired-end mode (70).
sRNA-and RNA-seq transcriptome assembly. The reference-based alignments were assembled into transcriptomes using the program stringtie to build full-length transcripts and calculate coverage and expression values (TPM). The assembly was guided by a gene annotation file from the H. volcanii DS2 (NCBI RefSeq taxid 2246) genome to build transcripts and annotate them as either a gene or novel transcript (71). The minimum distance between reads for transcript assembly was specified at 30 nucleotides. gffcompare under default options was used to compare the assembled transcriptomes against the gene annotation file from H. volcanii DS2 (NCBI RefSeq taxid 2246) to annotate transcripts as genes or noncoding RNA (antisense or intergenic) (72,73). In-house python scripts were used to bin transcripts that were annotated as genes, transcripts annotated as antisense (classified as a noncoding region opposite from a coding region), and transcripts annotated as intergenic (classified as a noncoding region between two coding regions), and antisense sRNAs were subsequently binned as 3= UTR, 5= UTR, or CDS overlapping.
sRNA-and RNA-seq differential expression analysis. We used a read count-based differential expression analysis to identify differentially expressed sRNAs during oxidative stress. The program featureCounts was used to rapidly count reads that map to the assembled sRNA transcripts (described above) (74). featureCounts was run with strand-specific options on, paired-end mode on, and multimapping off (74). The read counts were then used in the R differential expression software package DESeq2 (75). Briefly, read counts were converted into a data matrix and normalized by sequencing depth and geometric mean. Differential expression was calculated by finding the difference in read counts between the no-challenge state (control) normalized read counts and the oxidative stress normalized read counts (75). The differentially expressed sRNAs were filtered on the basis of the statistical parameter of false-discovery rate (FDR), and those that were equal to or under an FDR of 5% were classified as true differentially expressed sRNAs.
In silico validation of sRNAs. Differentially expressed sRNAs were validated by two in silico methods: (i) visualization of transcripts and (ii) open reading frame protein homology search. In the first method, transcriptomes for each replicate and condition were visualized on the Integrated Genome Viewer (IGV) against the H. volcanii (NCBI RefSeq taxid 2246) genome and annotation (76). The sRNA transcript coordinates were used to locate putative sRNAs, and if it was found within an operon, it was eliminated from further analysis. In the second method, blastx (default parameters) was used to search for protein and domain homology for each sRNA, and those that had significant homology with known proteins or domains were eliminated from further analysis (77).
Regulatory element motif identification of sRNAs. One hundred nucleotides upstream and downstream from the sRNA transcript start and stop coordinates were extracted using in-house python scripts. These regions were searched for transcription motifs (BRE and TATA box) using multiple sequence alignments, visualization with WebLogo (default parameters), and motif searching with MEME and CentriMo (default parameters) (78,79).
In vivo validation of sRNAs by Northern blot analysis. Twenty micrograms of total RNA and [ 32 P]ATP end-labeled Century-Plus RNA markers was loaded onto 5% denaturing urea polyacrylamide gels (SequaGel; National Diagnostics) and run at 30 W for 1.5 h to ensure well-spaced gel migration from 50 to 1,000 nucleotides (nt). Gels were transferred onto Ultrahyb nylon membranes and hybridized with 2 types of probes. For lowly expressed sRNAs, we probed with [␥-32 P]dATP randomly primed amplicons generated with custom primers based on sRNA transcript genomic coordinates as determined by the sRNA-seq in silico analysis. Probe primers were at a minimum of 10 nt inward from the predicted genomic coordinates (start and stop) to ensure accurate transcript detection. Hybridizations were performed at 65°C. To determine strandedness of sRNAs, we used [␣-32 P]dATP end-labeled oligonucleotide probes (20 to 23 nt) that were antisense to sRNAs. Hybridizations were at 42°C. The rpl30 protein (HVO_RS16975) transcript was used as a loading control for differential expression calculation, because it was not differentially expressed under oxidative stress in this RNA-seq data set. Differential expression was calculated using ImageJ.
Gene ontology enrichment analysis of mRNA-targets. NCBI gene names for all mRNA targets of antisense sRNAs were uploaded into Database for Annotation, Visualization, and Integrated Discovery (DAVID) to determine the pathways and gene ontologies targeted by sRNAs.
Accession number(s). All raw reads and processed data from these experiments are available at NCBI under BioProject PRJNA407425. Illumina raw sequence data (.fastq) for each replicate and condition are deposited in NCBI Sequence Read Archive with accession number SRP117726.

SUPPLEMENTAL MATERIAL
Supplemental material for this article may be found at https://doi.org/10.1128/JB .00779-17.

ACKNOWLEDGMENTS
This work was supported by grant FA9950-14-1-0118 from the AFOSR. We thank Madeline Cassani, Zhao Zhang, and German Uritskiy for advice and guidance on sRNA-seq library preparation, Evan Hass and Vidya Balagopal for advice on Northern blotting, Jacques Ravel, Mike Humphrys, and David Mohr for sequencing efforts and technical advice, and John Kim and Sarah Woodson for helpful discussions.