**Figure 6.** Provirus elements PVH3A1 and PVHS1 compared to ChaoS9. PVH3A1 from haloarchaeon strain 3A1\_DGR (accession NZ\_KK033114; nt 1475908-1520418), and PVHS1 from *Haa. sulfurireducens* strain M27-SA2 (accession NZ\_CP011564; nt 1971038- 1927125), are compared to ChaoS9. Grey bands connecting genome maps represent similarity (tBLASTx, >30% identity) between the inferred proteins of ChaoS9 and each provirus. Positions of several annotated proteins and their genes are indicated by color and name; TerL, large subunit terminase (red); Por, portal protein (light green); Head, SPP1\_gp7 family head assembly protein, (light brown); capsid proteins such as Mcp (major capsid protein) and Tail sheath protein (brown); Tpm, tape measure protein (yellow); Dam, adenine methyltransferase (light grey); Bpj, baseplate J family protein (blue); Orc1, replication protein Orc1 (green); RepH, plasmid replication protein (turquoise). Scale bar, in kbp, is shown at top.

The right arm (29–55 kbp) of ChaoS9 shows little similarity to either provirus except for Orc1. Across the left arm (0–28 kbp), the region covering virus structural and assembly proteins, ChaoS9 maintains good synteny and protein similarity to PVHS1, but with PVH3A1 there is a distinct break in similarity within the virus structural/assembly gene module. PVH3A1 neck and tail proteins are similar to ChaoS9 but the head/assembly proteins are unrelated, except for the head assembly

protein (labeled Head; Figure 6, panel a). PhiH1, and to a lesser extent phiCh1, show a comparable break in similarity to ChaoS9 head and tail genes (Figure 5). All of these examples suggest that ancestral recombination events between ChaoS9-like viruses have occurred within genes located between the head morphogenesis genes and the tail morphogenesis genes, and that these often result in viable progeny.

The distinct right halves of the proviruses suggest that they have also undergone extensive recombination relative to ChaoS9, perhaps reflecting differences in virulence and/or provirus state.

#### *3.7. A Diverse Family of Haloviruses*

In order to examine evolutionary relationships between ChaoS9 and other viruses and proviruses, genes were sought that were both sufficiently conserved and present in all examples. As shown in Figures 5 and 6, relatively few genes matched these criteria. Even the large subunit terminase (TerL) and major capsid protein (Mcp) were poorly conserved, and these have been widely used in previous studies of bacteriophages for delineating virus taxa [4]. The Bpj (baseplate J family protein) and the tail sheath proteins were selected to infer phylogenetic relationships as they were relatively long, conserved in sequence and present in all examples. BLASTp searches (accessed 20 November 2018) with ChaoS9 Bpj retrieved over 50 high scoring matches (E value < 10−30), with the closest matches all being from haloarchaea or haloviruses, while less similar matches included proteins from Bacteria and bacteriophages. An inferred phylogenetic tree based on alignment of Bpj (Figure 7a) shows that ChaoS9 is part of a robust clade (100% bootstrap confidence) which includes haloviruses phiCh1, phiH1, and provirus-like elements of four haloarchaea belonging to at least three different genera. A separate clade, also with high bootstrap confidence, contains four haloviruses (HF1, HF2, HRTV-8, and HSTV-2) and a provirus element of *Haloferax larsenii*.

BLASTp searches of the NCBI database with the ChaoS9 tail sheath protein (ChaoS9\_080) retrieved only nine high scoring matches (E value < 10−35), and these were all from haloarchaea or haloviruses present in the Bpj tree, and most closely related to ChaoS9. The inferred phylogenetic tree based on this protein (Figure 7b) reveals a topology similar to that of Bpj proteins.

The major capsid protein (Mcp) of ChaoS9 was used to search the NCBI database (BLASTp, nr database, accessed 22 November 2018) and retrieved only five matches (Figure 7c), which varied in similarity from 36 to 76% (aa identity). Three were from organisms previously identified as specifying ChaoS9-related proteins (*Haa. sulfurireducens*, *Saliphagus* sp. LR7 and *Hpt. malekzadehii*) and are present in Figures 6 and 7, while the other two sequences were from *Salinigranum rubrum* and the tailed halovirus HHTV-1 [5,54].

The DNA sequences of HHTV-1 and ChaoS9 share no significant similarity (Figure 4), but a BLASTp comparison of all ChaoS9 and HHTV-1 proteins found that Mcp was the only protein with significant similarity (36%) shared between these haloviruses. Such a pattern of similarity between tailed viruses that is strictly limited to the Mcp appears to be uncommon in the published literature. A less clear-cut but comparable example occurs between actinophages Jace and Tweety (accessions EF536069 and MH153804), which share similar Mcp (32%) and integrase (39%) protein sequences but weak or negligible similarity between all other proteins. Cases where head genes and tail genes derive from different virus lineages are slightly more common [55,56] (see Discussion).

**Figure 7.** Phylogenetic tree reconstructions using (**a**) baseplate J family proteins, (**b**) major tail sheath proteins and (**c**), major capsid proteins. Halovirus proteins are marked by stars and bacteriophage proteins are marked by triangles. The consensus trees were produced using the Neighbor–Joining algorithm and 100 bootstrap replications. Scale bars represent the estimated number of amino acid replacements per position.

#### *3.8. CRISPR Spacer Matches to ChaoS9*

The ChaoS9 genome was used to search for matching CRISPR spacer sequences available from public databases (see methods). A total of 277 spacers were identified that closely matched ChaoS9, with the majority of spacers originating from Antarctic lake metagenomes (Deep Lake, Rauer Lake and Club Lake). After removing duplicates, the number of distinct spacers was reduced to 39, and in Table S3 the matching spacers have been ordered by their position along the ChaoS9 genome. The distribution of spacers is highly skewed, with most (34/39) targeting sequences within the right arm, and particularly the *repH* gene, for which there were 14 distinct spacer sequences. The direct repeats (DR) of the matching spacers were most similar or identical to those found in CRISPR arrays of sequenced haloarchaea, particularly antarctic isolates such as *Halorubrum lacusprofundi*, *Halobacterium* sp. DL1 and haloarchaeon DL3. For comparison, the phiH1 and phiCh1 genomes were also scanned for matching spacers in the same antarctic lake metagenomes (IMG/VR webserver;

(https://img.jgi.doe.gov/cgi-bin/vr/main.cgi; accessed 20 December 2018). This returned only 2 (phiCh1) and 3 (phiH1) significant matches (supplementary Table S4), of which all of the phiH1 matching spacers were to the tail fiber gene (PhiH1\_165), as was one of the spacers matching phiCh1 (PhiCh1\_155). The remaining spacer to phiCh1 matched a sequence within the prohead protease gene (gpB, PhiCh1\_045). All six spacers were from Deep Lake and Rauer Lake metagenomes.

#### **4. Discussion**

This study focused on identifying the cause of a lysis event in a large, laboratory culture of *Hbt. salinarum* S9. A novel myovirus was recovered, ChaoS9, with morphological and molecular characteristics specifically resembling myoviruses phiH and phiCh1, but which differed significantly in sequence from both of these. An earlier lysis event in the same laboratory that affected a culture of *Hbt. salinarum* R1 was shown to have been caused by phiH, however, the 2007 event was not a recurrent infection. The sources of both infections are unknown, possibly raw salt, but these events highlight the need for preventative measures, even though the high salt conditions used for cultivation of haloarchaea are generally regarded as providing a strong barrier to contamination by non-halophilic microorganisms. However, tailed viruses (caudoviruses) such as ChaoS9 are not only the most common type of prokaryotic virus, but together with other bacterial and archaeal viruses, they represent the most abundant biological entity on Earth, estimated to be 1031 virions [57]. It is not surprising, then, that virus contamination and lysis events are a constant threat in large-scale commercial fermentations [22], and cause such significant losses that systematic preventative programs have been formulated, such as PACCP (phage analysis and critical control point) [58].

ChaoS9 was most similar in morphology, genome type, genome length, and sequence to the tailed haloviruses phiH1, and phiCh1 [26,43]. The dsDNA genome was terminally redundant and circularly permuted, consistent with headful packaging, as has been shown for the related viruses. The likely *pac* site was identified by sequence homology. Gene synteny of the virus morphogenesis genes of the left arm of the genome was similar to those of phiH1 and phiCh1, while the right arm comprised genes for DNA replication, plasmid partition, and a tRNA, as well as many genes specifying proteins of unknown function. The right arm corresponds to the lysogeny, replication, and accessory gene region of phiH1 and phiCh1, and probably serves the same general function in ChaoS9. A characteristic feature of phiH1 and phiCh1 is an invertible tail fiber gene module, which was also present in ChaoS9, except that it also included an ISH12 transposon.

Comparison of the ChaoS9 genome with those of phiH1 and phiCh1 revealed a distinctive pattern of similarity and difference, suggesting an evolutionary history involving large recombination events. While the tail gene region of ChaoS9 was similar to the other two viruses, the head/assembly genes of the left arm of the genome, as well as most of right arm (the replication/lysogeny region), were not. This means that the major capsid protein, terminase (large subunit) and portal proteins of ChaoS9 are all unrelated to those of phiH1 and phiCh1. The most parsimonious explanation is that a recombination event has replaced the head morphogenesis module while leaving the tail morphogenesis module intact. It appears to be a chimeric recombinant, with the point of similarity disjuncture occurring between the putative head closure (*hco*, ChaoS9\_055) and neck protein (*nep*, ChaoS9\_070) genes. Similar examples where head morphogenesis and tail morphogenesis modules appear to derive from different virus types are not uncommon in the literature, such as the *Gordonia* spp. phage Kita (Figure 4 of [55]), and *Xylella fastidiosa* phage Xfas53 [56], phi80 and Gifsy-2 [59]. The likelihood of such recombination events producing viable progeny is increased because of the way caudoviruses are assembled, with heads and tails produced as separate structures that are then joined together [60].

Provirus elements related to ChaoS9 were found among haloarchaeal genome sequences, integrated into their host chromosomes via tRNA genes, which supports the view that ChaoS9, like phiH1 and phiCh1, has a temperate lifestyle [61]. One provirus (PVHS1) showed good sequence similarity to ChaoS9 proteins across most of the virus morphogenesis genes, while a second provirus (PVH3A1) showed a general pattern of sequence similarity to ChaoS9 like that of phiH1 and phiCh1; with mainly unrelated head morphogenesis genes (including dissimilar major capsid proteins), similar tail morphogenesis genes, and largely unrelated right arms (lysogeny/replication). The right arms of both proviruses differed significantly from the right arm of the ChaoS9 genome, but all carried a putative *orc1* gene, presumably for use in replication as a (provirus) plasmid if the circularized genome is unable to integrate into the host chromosome via *attP*. These two examples again point to frequent recombination events that exchange the head morphogenesis module while retaining the tail morphogenesis genes.

Curiously, the ChaoS9 major capsid protein is most similar to that of the haloarchaeal siphovirus HHTV-1, but this is the only protein between them that is similar, and they are otherwise unrelated. In studies of enterobacterial viruses, major capsid protein (MCP) gene exchange was estimated to be extremely rare between phage clusters or types [61]. How this occurred in the evolutionary history of these viruses will be intriguing to resolve, but presumably indicates widespread recombination between tailed viruses of haloarchaea.

A surprisingly large number of identical or near identical CRISPR spacers to ChaoS9 were detected, and even more surprising was their identification in the metagenomes from Antarctic salt lakes. The majority of target sites were located in the right arm of the genome, in the lysogeny/replication region. A few were to sites within the left arm of the genome, such as the tail fiber gene. In contrast to ChaoS9, the genomes of phiH1 and phiCh1 matched only a few spacers from the same Antarctic metagenomes, and all but one targeted the tail fiber gene. This suggests that these lakes harbor tailed viruses with a lysogeny/replication module, similar to ChaoS9. Their tail fiber genes also show significant sequence similarity to ChaoS9, phiH1, and phiCh1, but the other virus morphogenesis genes are distinct.

The precise relationship between ChaoS9, phiCh1, and phiH1 was explored in several ways. They share average nucleotide identity (ANIb) values of ≥74%, and phylogenetic tree reconstructions using baseplate J virus protein sequences clustered ChaoS9, phiH1, and phiCh1 into a well-supported clade, distinct from other tailed haloviruses. The close relationship is consistent with the other comparative data (morphology, genome size and packaging, gene synteny, and inferred protein sequences). *Halobacterium virus phiH* is the type species of the genus *Myohalovirus* and a previous study has shown that phiCh1 should be included in the same genus, as their genomes share 63% nucleotide identity, are mostly colinear, and their proteins show, on average, 70% amino acid identity [26]. PhiCh1 proteins affected by genome revision became more similar to those of phiH1. ChaoS9 shows many similarities and correspondences to phiH1 and phiCh1 but also considerable differences, including a distinct major capsid protein, terminase (large subunit) and portal protein, and is certainly a distinct species to phiH1 and phiCh1, but whether it should be classified outside of the *Myohalovirus* genus requires further consideration.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/10/3/194/s1, Figure S1: Experimental and in silico restriction digests of ChaoS9 DNA. Table S1; DNA primers used for PCR or sequencing viral DNA. Table S2; Revisions of phiCh1 genes after re-sequencing the genome. Table S3; CRISPR spacers matching ChaoS9; Table S4, CRISPR spacers matching phiH1 and phiCh1.

**Author Contributions:** Conceptualization, visualization, investigation, M.D.-S.; investigation (electron microscopy) G.W.; investigation (genome sequencing) P.P.; data curation, F.P., M.D.-S.; funding acquisition, D.O.; project administration, F.P.; resources, A.W.; writing—original draft, F.P., M.D.-S.; writing—review & editing, A.W., D.O., F.P., M.D.-S.; validation, A.W.

**Funding:** This research received no external funding.

**Acknowledgments:** We wish to thank Cyril Boulegue, Beatrix Scheffer, and Snezan Marinkovic (core facility, MPI) for their contribution to genome sequencing. We thank Graham Hatfull for advice on bacteriophage diversity.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Complete Genome Sequence of the Model Halovirus PhiH1 (ΦH1)**

### **Mike Dyall-Smith 1,2, Felicitas Pfeifer 3, Angela Witte 4, Dieter Oesterhelt 1 and Friedhelm Pfeiffer 1,\***

Received: 11 September 2018; Accepted: 8 October 2018; Published: 12 October 2018

**Abstract:** The halophilic myohalovirus *Halobacterium virus phiH* (ΦH) was first described in 1982 and was isolated from a spontaneously lysed culture of *Halobacterium salinarum* strain R1. Until 1994, it was used extensively as a model to study the molecular genetics of haloarchaea, but only parts of the viral genome were sequenced during this period. Using Sanger sequencing combined with high-coverage Illumina sequencing, the full genome sequence of the major variant (phiH1) of this halovirus has been determined. The dsDNA genome is 58,072 bp in length and carries 97 protein-coding genes. We have integrated this information with the previously described transcription mapping data. PhiH could be classified into Myoviridae Type1, Cluster 4 based on capsid assembly and structural proteins (VIRFAM). The closest relative was *Natrialba* virus phiCh1 (ϕCh1), which shared 63% nucleotide identity and displayed a high level of gene synteny. This close relationship was supported by phylogenetic tree reconstructions. The complete sequence of this historically important virus will allow its inclusion in studies of comparative genomics and virus diversity.

**Keywords:** halovirus; virus; halophage; *Halobacterium salinarum*; Archaea; haloarchaea; halobacteria; genome inversion

#### **1. Introduction**

The temperate myovirus *Halobacterium virus phiH* (ΦH) infects the extremely halophilic archaeon *Halobacterium salinarum* strain R1 (DSM 671) and was isolated after the spontaneous lysis of a culture of its host [1]. Purified virions require 3.5 M NaCl for stability, have an isometric head of 64 nm diameter and a long, contractile tail (170 × 18 nm) with short tail fibres [1,2]. Virus preparations contain 3 major and 10 minor proteins [3]. The virus genome is linear dsDNA with a G+C content of 64%, contains a *pac* site, is about 3% terminally redundant and partially circularly permuted, and estimated to be 59 kb in length [4,5]. In the provirus state, the genome is extrachromosomal, covalently closed and circular, and 57 kb in length [4]. While always classified within the *Myoviridae*, the genus name has changed over the years from phiH-like viruses to *Phihlikevirus*, and most recently to *Myohalovirus* [6,7]. The species name itself has changed from *Halobacterium phage phiH* to *Halobacterium virus phiH* [6,7] but for convenience we will refer to it from here onwards simply as phiH, and the analysed variant as phiH1 or halovirus phiH1.

The original lysate of phiH was found to consist of a mixture of several distinct variants that appeared to have arisen from the activity of insertion sequences. The predominant variant, phiH1, was plaque-purified and a restriction map determined [5]. This was used for further study [3]. PhiH1 became a key model in the study of gene expression and regulation in haloarchaea and was instrumental in the development of genetic tools and methods in these extremophiles. Examples include the polyethylene glycol (PEG)-mediated transfection method [8], the pUBP1 cloning/expression vector [9], the identification of archaeal promoters, mapping transcription start and stop sites [10] and the analysis of gene regulation via repression [11,12]. The presence and function of antisense RNA in haloarchaea was first described in this virus [13]. An 11 kb invertible segment of the virus genome, called the L-region, was found to be flanked on one or both sides by the insertion sequence ISH1.8, and could also circularize to form a 12 kb plasmid (including one copy of ISH1.8), with subsequent loss of the remaining phage DNA [14]. A strain carrying this plasmid was immune to infection [14].

Unfortunately, work on phiH stopped in 1994 [15,16] but a related virus, *Natrialba* virus phiCh1 (ϕCh1), was described a few years later [17] and continues to be studied in the Witte laboratory [18,19]. PhiCh1 infects a haloalkaliphilic archaeon, *Natrialba magadii*, and the genomes of both host and virus are fully sequenced [18,20]. The provirus state of phiCh1 corresponds to plasmid pNMAG03 carried by *Nab. magadii*. A full comparison between phiCh1 and phiH1 was prevented as only parts of the phiH genome were ever determined. This deficit also prevented the inclusion of phiH in broad-ranging studies of virus diversity, taxonomy, and evolution. The aim of this study was to complete the phiH1 genome sequence and provide a thorough annotation. This will not only provide a better understanding of the results from previous studies on this virus but also allow complete genomic comparisons with a wealth of other datasets, including other sequenced viruses, haloarchaeal proviruses, metaviromic/metagenomic and environmental RNA sequences.

#### **2. Materials and Methods**

#### *2.1. Virus DNA and Sequencing Methods*

Purified phiH1 DNA [1] was originally provided to F. Pfeifer by Hans-Peter Klenk while both were working in the department of W. Zillig [21]. The DNA was stored frozen at −80 ◦C until use. Sequencing was performed in two stages. For the first stage, all available sequences of phiH1 were downloaded from National Center for Biotechnology Information (NCBI) [22] and imported into the Phred–Phrap–Consed package [23]. Overlapping sequences were assembled and primers designed to gather additional sequences using Sanger technology. This consisted either of primer-walking directly on virus DNA, or on polymerase chain reaction (PCR) amplimers, or PCR-sequencing across gaps. The resulting sequence reads were progressively assembled into contigs, base calls inspected manually and corrected where needed, and new primers designed for further rounds of sequencing until all gaps were closed. Except for overlaps, this approach left most of the previously published sequences unchecked.

In the second stage, short-read Illumina HiSeq sequencing of phiH1 DNA was performed (Max-Planck Genome Centre, Cologne, Germany). This returned 243 Mb of high quality sequence data (coverage = 4200-fold). De-novo assembly did not produce a single contig, due to short read-lengths and the presence of repeat sequences within the viral genome, but reads could be confidently mapped to the genome sequence obtained in the first stage (*Map to Reference* option; Geneious mapper method) in order to improve the sequence reliability.

#### *2.2. CRISPR Spacer Searches*

The crass v0.3.12 software [24] was used to extract CRISPR spacer sequences from genomic/metagenomic data available at the NCBI SRA database (accessed 27 July 2018) [25], as described previously [26]. These included all available genomes of members of the class Halobacteria, and metagenomes of hypersaline environments. CRISPR direct repeats (DR) identified by crass were used to search the CRISPRfinder database (accessed 25 July 2018) [27] for haloarchaea with matching or closely matching DR.

#### *2.3. Bioinformatic Methods*

Gene annotation used a combination of gene prediction with GeneMarkS-2 [28] and manual refinement using database searches (BLASTp/BLASTn; nr databases) at the NCBI webserver [29]. Repeats were identified by BLASTn, dot-plot comparison in Yass [30], and with tools within the Geneious software suite [31]. Circos plots were performed via the circoletto webserver [32]. Plots are coloured by the 'score/max' ratio of tBLASTx bitscores (real score/maximal score). Colours are: blue ≤ 0.25, green ≤ 0.50, orange ≤ 0.75, red > 0.75. Sequence mapping, alignments, editing and phylogenetic tree reconstructions were performed with Geneious software version 10.2 [31]. For phylogenetic tree reconstructions, protein sequences were first aligned using CLUSTALW, and trees inferred using the Neighbor-Joining algorithm (within Geneious). Consensus trees were determined after 100 bootstrap repetitions. Protein structural modelling used the I-Tasser webserver [33]. Identification of the *pac* site utilised the program PhageTerm [34] as implemented on the CPT Phage Galaxy [35]. The VIRFAM webserver [36] uses proteins of the phage head-neck-tail module to cluster phages into related groups, and was used to classify phiH1.

#### *2.4. Data Availability*

The phiH1 genome sequence has been deposited at Genbank under the accession MK002701. Raw reads were submitted to the SRA archive under accession SRP159490.

#### **3. Results and Discussion**

#### *3.1. Sequence and Annotation of PhiH1*

The previously sequenced regions of the phiH1 genome represented about 50% of the complete sequence (Figure 1, red lines). Using virus DNA as template, the gaps between these sequences were PCR amplified and Sanger sequenced. However, the quality of the previously sequenced regions was of uncertain reliability. High-coverage Illumina sequencing (ca. 4200-fold) was then used to enhance sequence confidence. Sequence revisions were only found to be required in previously deposited sequences but not to Sanger sequencing results of the first stage of the project. While the virus DNA found in capsid particles is linear, the head-full packaging process produces a population of molecules that are terminally redundant and partially circularly permuted [1]. The complete genome sequence determined in the current study is represented as the provirus form; a circular sequence of 58,072 bp. This value is close to the published size of 57 kb, estimated from restriction fragment sizes [4,37]. The G+C content of the genome was 63.7%, almost identical to the published value of 64% [3] but slightly lower than that of the host chromosome (68.0%) [38].

**Figure 1.** Diagram of the phiH1 genome with lines below showing regions previously sequenced (red) along with their database accessions. The blue lines (NEW) indicate regions sequenced in the present study by Sanger sequencing. Tick marks (dark green) below the blue lines show the positions of oligonucleotide primers used for PCR and primer-walking. Dots at the right and left contig ends indicate sequence continuity between them. Scale bar at top shows position in bp.

The original restriction map of phiH1 DNA, as determined by [5], corresponded closely with the *in silico* map inferred from the phiH1 genome sequence (Figure S1). The *pac* site located at the left end of the restriction map matched closely to the corresponding *pac* sequence of phiCh1. While the *pac* site of phiCh1 had been localized by restriction mapping [18], it had not been precisely mapped. For consistency, the start point of phiH1 was set to the corresponding start of phiCh1 even though this splits the *terS* gene. Using this numbering, the program PhageTerm [34] was used to analyse the mapping of Illumina reads to the phiH1 genome, and this located the *pac* site terminal base at nt 46, with high probability (*p* = 2.5 × 10−238). This is within the *terS* coding sequence (CDS) close to the stop codon and within a GC-rich region that is strongly conserved between phiH1 and phiCh1.

Annotation of the phiH1 genome resulted in 97 CDS (Table 1), most of which were encoded on the plus strand (86/97, Figure 2 panel b), and were frequently closely spaced, with 45 overlapping at start/stop codons and 23 separated by ≤8 nt. Many genes were in functional groupings typical of bacteriophages (Figure 2 panel b). The left end of the genome encodes DNA packaging proteins (e.g., terminase, portal protein), then virus assembly and structural proteins (e.g., major capsid protein, tape-measure protein, tail proteins). The three main proteins of purified virus were originally labelled by their estimated sizes on sodium dodecylsufate (SDS)-polyacrylamide gels (22, 53 and 80 kDa) [1], which were later revised to 27, 46 and 80 kDa [3] but in 1994, Stolt et al. [39] determined their N-terminal amino acid sequences and used this information to map the proteins (HP20, HP32 and HP67) to their genes (*hp20*, *hp32*, *hp67*) and sequence them. The inferred molecular weights (MWs) of proteins HP20 and HP67 were noted by these authors to be much smaller than previous estimates. In the present study, the locations of these genes on the full genome sequence have been resolved, an error in the *hp20* (accession X80161) coding sequence was corrected, and the MWs of the inferred proteins calculated (11.6, 35.4 and 45.5 kDa). For consistency we have retained the original gene names (Table 1).

The next genomic region is a replication/regulatory module (the L-region) that encodes RepR (repressor), a ParA-family protein (partition) and RepH (replication). There is also a VapC-like protein that together with the small overlapping upstream CDS may form a toxin-antitoxin pair that could be involved in plasmid maintenance [40]. The right end of the genome carries many genes with unknown function but includes genes specifying DNA methylases and cell lysis proteins. The taxonomic position of phiH1 was assessed using the VIRFAM webserver [36], which classifies bacteriophages and archaeal viruses based on the order and similarity of capsid assembly/structural proteins. Consistent with previous studies [6], phiH1 was classified by this system as a member of the Myoviridae (Type1, Cluster 4).

**Table 1.** *Cont.* 1PhiCh1/pNMAG03 homologs of phiH1 proteins show BLASTp E-values < 10−20. For phiCh1 proteins, both the PhiCh1p and originally assigned ORF codes (ORF for open reading frame) are shown (e.g., PhiCh1p02, ORF1). Codes starting with ORF represent the original annotation of the phiCh1 genome [17] (GB accession AF440695.1); and codes starting with PhiCh1p represent the RefSeq version of the annotation of the same genome sequence (GB accession NC\_004084). The number shift is due to the *terS* gene, the N-terminal part being encoded at the end of the genome, and the C-terminal part at its beginning. This ORF is complete in the provirus state due to circularization and in the linear virus state due to terminal redundancy. This gene is *ORF98* in the original annotation and PhiCh1p01 in the RefSeq annotation. Codes starting with Nmag\_ represent the annotation of the *Natrialba magadii* plasmid pNMAG03 [20] (accession CP001935.1). The point of ring opening in pNMAG03 was set between Nmag\_4303 and Nmag\_4211. Codes in square brackets represent NCBI accessions referring to homologous proteins (BLASTp E-values ≤ 10−11), which are from other sources. 2 Gene PhiH1\_185 is encoded on an invertible segment. In the current sequence version, it is inactivated because it is uncoupled from a start codon. By genome inversion, it becomes activated while its partner gene PhiH1\_165 becomes inactivated. Overall, this results in tail fibre protein switching. 3 This protein (PhiH1\_465) has three predicted transmembrane domains and has been suspected to function as a holin [18]. 4 In these cases, the phiCh1 gene is split into two CDS but is continuous in phiH1. 5 These proteins are more distantly related (show less than 39% sequence identity or fall above BLASTp E-values of 10−20). In these cases, a similar genetic context supports their stated relationship.

A GC-profile plot [41] of the phiH1 genome shows a major low point inflection within the L-region (Figure 2, panel a), indicating a potential replication origin. The L-region is ~12 kb in length, can replicate as a plasmid in *Halobacterium* [14], and carries genes encoding a replication protein (RepH), and a DNA-binding repressor (RepR). It can also provide cells with immunity to infection by phiH1 virus. The transcription program of phiH1 during lytic growth (panels c, d and e) has been summarized from previous studies, and shows temporal changes (early, middle and late transcripts). The broad directions of transcription reflect the closely spaced and similarly directed gene clusters as well as the correspondence with functional gene groupings (panel b). The lowest two panels (d, e) summarize the results of hybridizing labelled transcripts from infected cells to Southern blots of restricted phiH1 DNA [42], so mapping transcripts to fragments of the virus genome. Panel c shows a summary of the virus-specific transcripts that were sized by agarose gel electrophoresis and had 5 start sites mapped. While transcription across the L-region has been examined in more detail compared to the rest of the genome, there remains much that is incomplete or uncertain. For example, the 3 end of the late transcript labelled TLL, which is depicted ending in a dotted line and question mark (at ~21 kb), has not been determined. This transcript could potentially extend for another 5.5 kb. Counter-transcripts are commonly produced by prokaryotes and their viruses and play important roles in gene regulation. Their presence and activity in phiH1 gene expression has been studied and was one of the first reports of antisense RNA in Archaea [13]. However, this interesting topic remains to be fully explored.

**Figure 2.** PhiH1 GC-profile, genetic map, and corresponding transcription program (adapted from [16,42]). (**a**) GC-profile of the phiH1 genome. (**b**) Genetic map of the phiH1 genome, showing coding sequences as red, blue or grey arrows. Dotted lines above indicate gene clusters involved in particular functions. Some CDS are labelled above the map, e.g., TerL, terminase large subunit; Portal, portal protein; Tape, tape-measure protein; RepH, replicase (label within CDS arrow); Mt, DNA methylases; TerS, terminase small subunit. Some genes are shown below the map, such as *hp32*, encoding the major capsid protein HP32. Panels c, d and e summarise transcription data from previously published studies, and above them is a colour key that indicates the time of appearance of early (0–1 h, blue), middle (1–2 h, green) and late (>2 h, pink) transcripts. (**c**) Precise mapping of viral transcripts, including start and termination sites [16,39]. (**d**) Summary transcription program of lytic

infection based on hybridisation of labelled infected-cell transcripts to restriction fragments of virus DNA [42]. Thin coloured lines indicate whether continuing transcription persists over time. (**e**) The transcription map data of [42] are shown, projected onto the in silico restriction map of phiH1, as determined from the complete genome sequence (this study). Enzymes are indicated at the left. Numbers on the restriction map refer to those of the original publication of [42] (see also Figure S1). Coloured shading follows that of panels c and d. Dotted pattern shown beyond the right-hand *pac* site indicates terminal redundancy of virus DNA. Scale bars (in kb) are shown below panels a and e.

Corrections to the previously sequenced regions resulted in significant changes to several coding sequences. For example, the *tnpB* gene of transposon ISH1.8 (nt 41,906–43,789) was thought to be inactive as it was split into three CDS by multiple mutations [43]. The high-quality Illumina sequence data show, however, that the gene is intact and that the previously reported transposon ISH1.8 (X00805) is actually an exact copy of transposon ISH12 from the host *Hbt. salinarum* strain R1 [44]. The element plays a key role in the mobilisation of the L-region of the genome to form the 12 kb plasmid, pΦHL [14,43]. Another case is the Dcm5 cytosine methylase, which was also reported as being split [39]. The revised sequence shows that the gene codes for a single, probably functional protein (PhiH1\_405, nt 47,732–49,618) and not for the two parts (dcm5a, dcm5b) as previously reported. Although phiH1 carries three potentially active DNA methylase genes (*dcm5*, *yhdJ* and *ycdA*), the presence of modified bases in phiH1 DNA was not detected in the chromatographic (high-pressure liquid chromatography) profiles of deoxyribonucleosides released by enzymatic hydrolysis [45]. In that study, the genomes of phiH and another, unrelated virus (phiN) were analysed, and while unmodified dC was detected in phiH, the phiN genome contained only methylated dC (Figure 3 in [45]). The related halovirus phiCh1 carries homologs of two of the phiH1 methylases, and one of them, N6-adenine methylase (ORF94/M·ϕCh1-I, corresponding to YcdA of phiH1) has been shown to methylate DNA at GATC motifs [46] but the proportion of sites found to be modified in virus DNA by M·ϕCh1-I varies from 5% to 50%, depending upon the infection conditions. Modifying only some of the available sites is presumably advantageous to avoid host restriction, as distinct enzymes may target either unmethylated or methylated sites.

#### *3.2. Matches to CRISPR Spacers*

The phiH1 genome was used to search for matching CRISPR spacers among metagenomic datasets of hypersaline environments downloaded from the NCBI sequence read archive (SRA; see methods). Only four spacers showing close to moderate similarity to phiH1 were detected (Table 2). These spacers match to virus genes encoding structural and non-structural proteins, and the DRs of these spacers show that they are carried by haloarchaea. The datasets include metagenomes from the USA and Iran, as well as an isolate from the Andaman Islands, India. The results suggest that phiH1-like viruses are geographically widespread.

**Table2.**CRISPRspacersmatchingphiH1.

phiH1, and asterisks indicated identical bases. DR: direct repeat (with haloarchaea containing most closely matching DR shown in brackets). 2 Symbols under alignment (\*:.)identical,similarandweaklysimilarresidues,respectively(basedonGonnetPAM250matrix).

 indicate

#### *Genes* **2018** , *9*, 493

#### *3.3. Relatives and Phylogeny of PhiH1*

The only close matches to the phiH1 genome in the GenBank database were phiCh1 and the corresponding *Nab. magadii* plasmid pNMAG03 (BLASTn, accessed 20 July 2018). A dot-plot comparison of phiH with phiCh1 (Figure S2) revealed a largely colinear relationship (green line) and an overall nucleotide similarity of 63%. The plot also highlights several indels (line gaps) and two regions showing inversions (red lines). Inversion 1 (nt 24,227–27,767) corresponds closely in sequence and arrangement to the invertible region described in phiCh1 (ORF34-36) that has a central XerD type integrase/recombinase gene flanked by inverted repeats, and facilitates switching between two related tail fibre genes, each containing numerous short repeats [18]. The phiH1 orthologous integrase is PhiH1\_175. In the current sequence version, PhiH1\_165 is active while PhiH1\_185 is uncoupled from a start codon and thus is inactivated. Upon inversion of the genome segment, PhiH1\_185 is activated while PhiH1\_165 is inactivated. Overall, this results in tail fibre protein switching, which may affect receptor binding specificity and host range of phiH1. The similarity in the tail fibre protein repeats is high enough to be detectable at the DNA level, which results in the X-shaped pattern for this region in the dot-plot. Inversion 2 (nt 31,932–34,126) occurs within the phiH1 L-segment, encompasses four CDS including a ParA-domain protein, and is nearby a different integrase/recombinase gene (Int2, PhiH1\_240). Protein searches (BLASTp) of the phiH1 genome returned matches to phiCh1, a limited number of haloarchaeal genomes (often 5–10, which may flag proviral regions) and the haloarchaeal caudoviruses BJ1 [47] and CGphi46 (NC\_021537), both of which infect *Halorubrum* spp. Pairwise alignments (BLASTp) between all phiH1 and matching phiCh1 proteins gave an average protein sequence identity of 70% (range 39–95%; with a few exceptions, see footnote 5, Table 1). Figure 3 is a graphical comparison of phiH1 proteins (tBLASTx) with those of phiCh1, BJ1, CGphi46 and, as an outlier, HSTV-1. The *Haloarcula* caudovirus HSTV-1 [48] shows very low similarity to phiH1. The figure summarizes the close similarity of phiH1 and phiCh1 proteins. BJ1 and CGphi46 show far fewer matching regions, mainly to proteins encoded near the left end of the phiH1 genome, a region specifying portal and capsid proteins. The three significant matches to HSTV-1 were to a methyltransferase (HSTV1\_52), a hypothetical protein (HSTV1\_53), and a DNA polymerase sliding clamp protein (HSTV1\_40).

While several caudovirus proteins have been used to infer virus phylogenies, the major capsid protein (MCP) is often used because of its functional constraints maintaining a conserved structure [49]. Figure 4 shows a tree reconstruction using an alignment of phiH1 MCP (HP32, 35.4 kDa) and related sequences. Haloarchaeal proteins are seen to branch together (pink shading) and within this cluster the phiH1 and phiCh1 MCPs form a distinct and closely branching clade. These two proteins share 82% amino acid identity. The MCPs of CGphi46 and BJ1 branch at distant locations from each other and from phiH1 MCP. Structures of close homologs of phiH1 HP32 have not yet been determined. However, the major capsid proteins of bacterial caudoviruses and eukaryotic herpesviruses share a common folding structure, the archetype of which is the phage HK97 MCP (gp5) [50]. Consistent with this, modelling of the phiH1 MCP (I-Tasser) returned bacteriophage HK97 gp5 (PDB 2fs3A) as the closest matching structure (Template Modeling (TM)-score = 0.848, Root-Mean-Square Deviation (RMSD) = 1.17). Based on structure prediction and homology modelling, the HK97-fold may also be present in the MCP of phiCh1 [49]. The structure of the MCP of the haloarchaeal podovirus, HSTV-1, has recently been shown to be of the HK97 type [48].

**Figure 3.** Circos plot of amino acid similarity (tBLASTx) between phiH1 and the haloviruses phiCh1, BJ1, CGphi46 and HSTV-1. The threshold for connecting lines was E-value ≤ 10−40, with line colours reflecting the ratio of actual tBLASTx bitscore to the maximal score (using 'score/max' ratio colouring with blue ≤ 0.25, green ≤ 0.50, orange ≤ 0.75, red > 0.75). The outer histogram counts how many times each colour has hit the specific part of the sequence and uses an equivalent colouring scheme. The distance between successive tick marks shown along each virus genome represents 0.1 of the full genome length. Protein names shown along the phiH1 genome indicate the positions of the corresponding genes.

PhiH1 and phiCh1 display a close sequence similarity across most of their genomes yet infect physiologically and biochemically different haloarchaeal hosts. *Hbt. salinarum* is a widely distributed neutrophilic heterotroph with glycolipid-containing membranes, and has often been isolated from spoilage of salted products while *Nab. magadii* is a haloalkaliphile (optimum pH 9.5) that lacks glycolipids [51] and is restricted in its distribution to highly alkaline salt lakes [52]. Looking more widely, the presence of phiH1 MCP homologs in diverse genera of haloarchaea and two haloviruses (Figure 4) indicates that the *Myohalovirus* genus and related viruses are a highly successful group, the reasons for which are worthy of more detailed study, particularly when large-scale cultivation of *Halobacterium* becomes more common [53]. PhiH1 has been well studied in the past, and the completion of its genome sequence now allows it to be included in much of the sequence-based studies used today, including comparative virology, detection of proviruses in archaeal genomes, virus evolution and the microbial ecology of hypersaline environments.

**Figure 4.** Phylogenetic tree reconstruction (NJ method) of major capsid proteins (MCP) of phiH1, other haloviruses and related proteins of haloarchaea. Species names of haloarchaeal species are shown, with accession numbers given at the right side. Bootstrap confidence values (100 repetitions) are shown at branch points. The pink shading highlights taxa belonging to the class Halobacteria. Scale bar (expected changes per site) is shown at top. The outgroup (not shown) consisted of distantly related MCP sequences of *Bacillus* spp. (WP\_001060157.1, WP\_098773561.1, WP\_001064748.1 and WP\_000178926.1).

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4425/9/10/493/s1, Figure S1: Original phiH1 restriction map compared to in silico map from genome sequence. Figure S2: Dot plot sequence comparison of phiH1 and phiCh1 genomes.

**Author Contributions:** Conceptualization, visualization, investigation, M.D.-S.; data curation, Fr.P., M.D.-S.; funding acquisition, D.O., Fr.P.; project administration, Fr.P.; resources, A.W., Fe.P.; writing—original draft, Fr.P., M.D.-S.; writing—review & editing, A.W., D.O., Fe.P., Fr.P., M.D.-S.; validation, A.W.

**Funding:** This research was funded by the Max Planck Society, Germany, to Dieter Oesterhelt, emeritus Director of the Department of Membrane Biochemistry, Max Planck Institute, Martinsried, Germany.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Back to the Salt Mines: Genome and Transcriptome Comparisons of the Halophilic Fungus** *Aspergillus salisburgensis* **and Its Halotolerant Relative** *Aspergillus sclerotialis*

#### **Hakim Tafer, Caroline Poyntner \*, Ksenija Lopandic, Katja Sterflinger and Guadalupe Piñar**

VIBT EQ Extremophile Center, Department of Biotechnology, University of Natural Resources and Life Sciences, Muthgasse 18, 1190 Vienna, Austria<EMAIL_ADDRESS>(H.T<EMAIL_ADDRESS>(K.L<EMAIL_ADDRESS>(K.S<EMAIL_ADDRESS>(G.P.)

**\*** Correspondence<EMAIL_ADDRESS>
Received: 22 April 2019; Accepted: 15 May 2019; Published: 20 May 2019

**Abstract:** Salt mines are among the most extreme environments as they combine darkness, low nutrient availability, and hypersaline conditions. Based on comparative genomics and transcriptomics, we describe in this work the adaptive strategies of the true halophilic fungus *Aspergillus salisburgensis*, found in a salt mine in Austria, and compare this strain to the ex-type halotolerant fungal strain *Aspergillus sclerotialis*. On a genomic level, *A. salisburgensis* exhibits a reduced genome size compared to *A. sclerotialis*, as well as a contraction of genes involved in transport processes. The proteome of *A. sclerotialis* exhibits an increased proportion of alanine, glycine, and proline compared to the proteome of non-halophilic species. Transcriptome analyses of both strains growing at 5% and 20% NaCl show that *A. salisburgensis* regulates three-times fewer genes than *A. sclerotialis* in order to adapt to the higher salt concentration. In *A. sclerotialis*, the increased osmotic stress impacted processes related to translation, transcription, transport, and energy. In contrast, membrane-related and lignolytic proteins were significantly affected in *A. salisburgensis*.

**Keywords:** halophiles; halotolerant; fungi; transcriptomics; genomics

#### **1. Introduction**

Saline and hypersaline environments like salt marshes, saline soil, and salt water, as well as the Dead Sea harbor a diverse community of fungi [1]. Most of these fungi do not require salt for growth and have optimum growth in the absence of salt. Nevertheless, they are halotolerant and can withstand a salt concentration up to 30% [2].

Real halophilic microorganisms are adapted to conditions of high salinity and require a certain concentration of NaCl for their optimum growth. Halophiles are phylogenetically quite diverse and can be found in the domains of Archaea, Bacteria, and Eukarya. Until almost a decade ago, it was a general belief in mycology and in food microbiology that fungi growing on substrates with low water activity (aw) have a general xerophilic phenotype [3] that is determined by the water potential of the medium, rather than by the chemical nature of the solute [4,5]. Therefore, fungi were considered xerophilic if they grew well at an aw of 0.85, corresponding to 17% NaCl or 50% glucose added to their growth medium. Xerophilic Aspergilli, together with Penicillia, dominate the actively-growing mycoflora on dried food and also on materials stored in museums and archives. *Aspergillus vitricola* was first isolated from glass surfaces [6] and occurs frequently in house dust. Furthermore, *Aspergillus glabripes* occurs on books and archive material, as well as in house dust [7]. Xerophilic halotolerant Aspergilli were also shown to be common invaders of pipe organs in churches [8].

The term "halophile" for fungi was introduced in 1975 for those few xerophilic food-borne species that exhibit quite superior growth on media with NaCl as the controlling solute. Fungi have subsequently been described in moderately-saline environments, such as salt marshes, saline soil, and sea water, but were considered to be unable to grow in highly saline waters [2]. In 2000, fungi were isolated for the first time from the brine of solar salterns [2,9]. However, halophilic fungi are uncommon [2]. Currently, only *Wallemia ichthyophaga*, *Wallemia muriae*, *Aspergillus baarnensis*, *Aspergillus salinarium*, and *Aspergillus ruber* [10] are classified as halophiles. Recently, a fungus isolated from a historical wooden staircase in a salt mine in Austria was described as a new species, *Aspergillus saliburgiensis* [11,12], and was added to the list of halophilic fungi.

Cells living in natural saline systems, where high salt amounts cause high osmotic pressure, must maintain lower water potential than their surroundings in order to survive and proliferate. The ability to survive osmotic stress requires several adaptations. *Wallemia ichthyophaga* withstands high salt concentration by increasing the intracellular concentration of polyols, using high-affinity K+ transporters and increasing the thickness of its cell-wall by a factor of three [13]. Under high salt concentration, *Aspergillus ruber* increases the number of ion transporters, exhibits a higher proportion of acidic amino acids, restructures its cell wall, and uses glycerol as a compatible solute [10]. Similar adaptation strategies are found in halotolerant species, where the production of hydrophilic compounds, such as amino acids, sugar alcohols, and soluble sugars, have been reported [14–18].

The mechanisms behind the adaptive capacities of halotolerant and halophilic fungi have only begun to be understood in recent years, with the first study of the transcriptome of *Wallemia ichthyophaga* dating from 2013 [13]. Since then, an additional halophile, *Aspergillus ruber*, has been published [10]. In order to extend the knowledge in the field, we present the results of the comparative sequencing of the genomes of two species, *A. salisburgensis*, a halophilic fungus from a salt mine, and *A. sclerotialis*, as a representative of halotolerant fungi. Additionally, the two species were exposed to two concentrations of salt (5% and 20%), and the cellular response was studied on the transcriptome level. The genetic content, as well as the transcriptome, were compared to other halophilic, halotolerant, and to other well-studied, less salt-tolerant fungi (named as control) in order to better understand which mechanisms were involved in the resistance against osmotic stress.

#### **2. Materials and Methods**

#### *2.1. Strains, Media, and Cultivation*

The strain *Aspergillus salisburgensis* (EXF-10247/MA6005) was isolated in a previous study from a wooden staircase, built in the year 1108 B.C. and discovered in 2003 in a salt mine in the Austrian region "Salzkammergut", Upper Austria [12]. The strain *A. sclerotialis* (Strain No. MA 5985, CBS 366.77) was supplied by the Fungal Biodiversity Centre (CBS), Utrecht, Netherlands. The two examined strains were grown on salt glucose media for the whole genome sequencing (5% glucose, 0.1% peptone, 0.1% malt extract, agar, and either 5% or 20% NaCl) and salt yeast media for the transcriptome sequencing (0.4% glucose, 0.4% yeast extract, 1.0% malt extract, 1.2% agar, and either 5% or 20 % NaCl). The two strains were incubated at their optimal temperature (*A. salisburgensis* at room temperature, *A. sclerotialis* at 37 ◦C) for 35 days.

#### *2.2. Whole Genome Sequencing*

Genomic DNA extraction of pure fungal strains was performed as previously described [19]. Approximately one cm2 of fungal colonies, consisting of both hyphal and conidiogenous structures, were collected and placed in a 1.5-mL bead beater tube with 0.2-g glass beads (0.75–1 mm, Carl Roth GmbH Co., KG, Karlsruhe, Germany) and 500 μL of CTAB-buffer (1.2 g Tris-HCl, 8.2 g NaCl, 0.81 g EDTA × 2 H2O, 2.0 g CTAB, and 0.2 g β-mercaptoethanol, pH 8.0) and processed twice in the Fast Prep FP120 Ribolyzer (Thermo Savant, Holbrook, AZ, USA) for 40 s at a speed of 4 m/s. Between these ribolyzing steps, samples were incubated at 65 ◦C for 10 min at 300 rpm. Further DNA extraction

was performed with phenol/chloroform/isoamyl alcohol. Genome sequencing was carried out using the ION Proton Technology (Ion AmpliSeq Library Preparation kit, Template OT2 200 kit, Ion PI Sequencing 200 kit, Ion PI chip kit V2, Life Technologies, Carlsbad, CA, USA) following the instructions of the manufacturers.

#### *2.3. Total RNA Extraction*

After 35 days of cultivation, colonies were collected as described in Section 2.2, and RNA was extracted from the biomass using the FastRNA Pro RED Kit (MP Biomedicals, Santa Ana, CA, USA). For each condition and strain, three replicates of total RNA were extracted following the manufacture's protocol. The quality and quantity was measured using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and a Qubit 2.0 (Life Technologies, Carlsbad, CA, USA).

#### *2.4. RNA Library Preparation*

The total RNA library was prepared as described before [20]. mRNA was isolated from the total RNA using the Dynabeads mRNA DIRECT Micro Kit (Ambion, Life Technologies, Carlsbad, CA, USA). Subsequently, the RNA library for sequencing was constructed using the Ion Total RNA-Seq kit v2 (Life Technologies, Carlsbad, CA, USA). Quality was checked using the Agilent 2100 Bioanalyzer instrument (Agilent Technologies, Santa Clara, CA, USA). The final library was size selected to 290 bp using Pippin Prep (Sage Science, Beverly, MA, USA). The sequencing was done using the Ion Torrent Chef instrument, Ion Torrent Proton instrument, and the HiQ sequencing kit (Life Technologies, Carlsbad, CA, USA).

#### *2.5. Genome Annotation and Assembly*

*A. salisburgensis* and *A. sclerotialis* were assembled with Newbler. Genome completeness was assessed with BUSCO [21]. The final protein annotation was obtained by passing the ab initio protein annotation generated by BUSCO with the help of Augustus [22] to BRAKER2 [23,24] together with the RNA sequencing data. The functional annotation was done with InterProScan 5.25.64 [25] and EggNog [26] with default settings. Homologies with proteins in the Transporter Classification Database (TCDB) [27], the peptidase database (MEROPS) [28], and the Carbohydrate-Active Enzymes Database (CAZY) [29] were assessed with Blastp [30] (E-value <0.001).

#### *2.6. Comparative Genomics*

Protein orthology was assessed with ProteinOrtho 5.16b [31], where BLAST was replaced with DIAMOND [32]. The species phylogenomic tree was generated with iqtree [33–35] based on the MAFFT-alignment of one-to-one homologues. The TimeTree was inferred using the Reltime [36,37] method and ordinary least squares estimates of branch length. The TimeTree was computed using 8 calibration constraints retrieved from TimeTree.org. All position containing gaps and missing data were eliminated. There was a total of 37,607 positions in the final dataset. Evolutionary analyses were conducted in MEGA X [38].

Protein family expansion and contraction were computed with CAFE [39] and the previously-obtained TimeTree. The significance of protein families over- and under-representation was computed with the Fisher exact test.

Amino acids in the one-to-one conserved proteins and proteins containing a signal peptide, as annotated by SignalP [40], were tested for distribution bias with Wilcoxon tests in R [41].

#### *2.7. Differential Expression*

Genes' differential expression was obtained with kallisto [42], tximport [43] and DeSeq2 [44]. All figures were made in R with ggplot2 [45], ggtree [46], and deeptime.

#### **3. Results and Discussion**

#### *3.1. Highly Complete Genome and Gene Set*

The genome assemblies of *A. salisburgensis* and *A. sclerotialis* used for this study were sequenced on the Ion Proton sequencing platform (Sections 2.2–2.4). The assembled genome of *A. salisburgensis* had a size of 21.9 Mb and contained a total of 8895 protein-coding loci for a total of 1351 contigs (>500 bp). *A. sclerotialis* had a 27% larger genome (27.9 Mb) with 27% more coding-genes (11307), compared to *A. salisburgensis*, and 1704 contigs (Supplementary Table S1). Among the 4406 highly-conserved genes in Eurohiomycetes [21], 96% were found in *A. sclerotialis* and 94% were present in *A. salisburgensis*. It was reported that compact genomes with a reduced number of protein-coding genes might be a characteristic of obligate halophiles, such as *Wallemia ichthyophaga* [15]. This trend could not be confirmed in Aspergilli where the genome of the halophile *Aspergillus ruber* [10] has a similar size and gene content compared to *A. sclerotialis*.

#### *3.2. Gene Conservation*

The gene conservation pattern between *A. salisburgensis*, *A. sclerotialis*, two halophiles (Table 1, marked as HH), six halotolerant (marked as H), and three control strains exhibiting no salt resistance (marked as C) were studied with ProteinOrtho [31].

**Table 1.** List of genomes compared in this study. For each species, the name, the range of salt tolerance, the type of tolerance, and reference are listed. H: halotolerant, HH: halophile, C: control.

ProteinOrtho identified 1169 groups of orthologs present in all 13 species, 1091 genes specific to *A. salisburgensis*, and 2362 genes specific to *A. sclerotialis*. The genes specific to *A. salisburgensis* were enriched in functional terms related to protein refolding (GO:0042046, GO:0061077, *HSP20, HSP70, HSP90, CLP, DNAJ/HSP40*), DNA modification (GO:0015074, GO:0005727, GO:0046821, IPR01584, IPR025668, PF13683, GO:0005727, GO:0046821, TCDB 3.A.7), tripartite ATP-independent periplasmic (TRAP), and ATP-binding cassette transporter (IPR010656, IPR000515, TCDB 3.A.1, 2.A.56), as well as chitin degradation (CAZY:GH23) (Figure 1). TRAP transporters were shown to mediate the uptake of compatible solutes in halophilic Proteobacteria such as *Halomonas elongata* [50], while the proteins related to refolding are known to be involved in stress-response [51].

**Figure 1.** Graphical representation of the −*log*10(*f dr*) of the enriched categories in CAZY (Carbohydrate-Active Enzyme Database), Gene Ontology (GO), InterProScan (IPR), Protein Family (PFAM) and Transporter Classification Database (TCDB) in the set of *A. salisburgensis*-specific genes. Functional categories related to chitin degradation, protein folding, transport, and DNA integration are overrepresented.

In *A. sclerotialis*, the set of species-specific genes are functionally enriched in terms related to protein folding, heat shock protein 70, and helix-turn-helix (HTH)-binding domain protein (Supplementary Figure S1 and Supplementary Table S1). The helix-turn-helix binding-motif is found in many proteins involved in gene expression regulation [52].

#### *3.3. Gene Family Evolution*

A TimeTree of life for the 13 fungi (Table 1) was constructed from the multiple sequence alignments of the proteins shared by all strains and from nine time anchors extracted from the TimeTree database [53] (Figure 2). The results were consistent with the currently-accepted taxonomy [53]. *A. salisburgensis* and *A. sclerotialis* separated approximately 41 MYA (+/−40 MYA), 20 MY before the split between *Debaryomyces fabryi* and *Debaryomyces hansenii* and 40 MY after the split between *A. salisburgensis*, *A. sclerotialis*, and *Aspergillus ruber*. Interestingly, *A. salisburgensis* and *A. sclerotialis* exhibited the largest evolutionary rates, followed by the two *Debaryomyces* species.

This tree was used to estimate the numbers of rapidly-evolving annotation elements with CAFE [39]. The functional annotations used were derived for the carbohydrate enzymes (CAZY) [29], general enzymes from the EggNog annotation [26], InterProScan protein domains (IPR) [54], KEGG pathways [55], proteolytic enzymes (MEROPS) [28], PFAM protein domains [56], and transporters from TCDB [27].

**Figure 2.** Phylogenomics analyses. (**A**) TimeTree generated from the merged alignment of the one-to-one homologues of the 13 species. The green bars represent the 95% confidence interval. Geologic periods are color-coded and indicated at the bottom of the tree. Evolutionary rate is color encoded. (**B**) CAFE analysis of the expansion and contraction of the annotation elements of KEGG, CAZY, TCDB, Enzyme, and MEROPSfor the 13 species. The suffixes Mand Pafter the annotation categories stand for contraction and expansion, respectively.

The internal node with the largest numbers of rapidly-evolving annotation families was the one leading to the common ancestor of *Hortaea* and *Aspergillus* (*leotiomyceta*) (Supplementary Table S1). CAFE reported 106 IPR and 105 PFAM rapidly-evolving families. The most significantly expanding families were related to fungal-specific transcription factors (*p*-value 4.75 × 10−19 and cytochrome P450 (*p*-value 1.127 × 10−16). Cytochrome P450 was also shown to improve salt tolerance in plants [57] and is upregulated under salt stress in the fungus *Piriformosa indica*, a plant endophyte that confers salinity stress tolerance in rice plants [58]. An additional significantly-expanding family was DUF3468 (DUF, domain of unknown function, *p*-value 1.59 × 10−15), a protein domain involved in transcription activation of genes related to asexual conidiation and sexual differentiation in *Aspergillus nidulans* and *Aspergillus flavus* [59].

The branch leading to *A. salisburgensis* exhibited the largest contraction of gene families from all species analyzed: 91 IPR, 74 PFAM, 7 CAZY, 8 TCDB, 1 MEROPS, and 8 Enzyme families contracted significantly, which is in line with the reduction in genome size of 27% compared to *A. sclerotialis*. Fungal transcription factors (PF11951) exhibited the largest loss in *A. salisburgensis* (−17). Protein families transporting amino acid and oligopeptide (2A3, −11; IPR002293, −9; 2A18, −5; PF03169, −3), fatty acid (4C1, −7), and iron (2A108, −5) were contracting in *A. salisburgensis*. Interestingly, the amino acid polyamine organocation superfamily (APC) contracted also in *Aspergillus ruber* (2A3, −8). APC is one of the largest families of secondary active transporters, and it is found in all living organisms [60]. Another family depleted in both *A. salisburgensis* and *Aspergillus ruber* is aminoglycoside phosphotransferase, a bacterial antibiotic resistance protein (IPR002575) with a reduction of 15 and 13 members, respectively.

A few annotation families expanded in *A. salisburgensis*. These families were related to the type IV secretory pathway, which exports proteins or DNA-protein complexes out of the cell [61] (3A7 +9), the TRAP-C4 transporter (IPR010656 +6, PF06808 +6), the fungal cellulose binding domain (PF00734, +5), the integrase catalytic core (IPR001584), cutinase (PF01083), sodium:bile-acid symporter/arsenical resistance protein ACR3 (PF01758 +2), and secretory lipase (PF03583, +2). The increase in cellulose

degrading ability due to the fungal cellulose binding domain [62] and cutinase [63] is probably related to the environment, a wooden staircase, from which *A. salisburgensis* was isolated. In yeasts, proteins belonging to the PF01758 family confer arsenic resistance by extrusion of sodium arsenate and sodium arsenite [64].

#### Gene Family Enrichment

The significance of size variations of functional annotation families was studied by comparing species or a group of fungal species (Table 1) with the help of chi-squared tests. The most significant bias was seen for the major facilitator superfamily (MFS) transporters, which showed a peculiar pattern of enrichment. *A. salisburgensis* had significantly more MFS (2A1, IPR020846) transporters than *Wallemia ichthyophaga*, but significantly less than *A. sclerotialis* and *Aspergillus ruber* (Figure 3 and Supplementary Table S1). Further, there were significantly more MFS transporters in halotolerant fungi than in the halophiles (fdr = 0.018). Still, compared to the group of control fungi, both the halophilic (fdr = 3.76e−15) and the halotolerant (fdr = 7.19e−37) fungi were enriched in MFS.

Genes involved in protein translation, such as ribosomes (GO:0003735, GO:0005840), RNA binding (GO:0003723), and reverse transcriptase (GO:013103, GO:0015074), were depleted in halophiles and halotolerant compared to the control group. In contrast, GO terms involved in transcription, such as RNA polymerase II (GO:0000981), transcription (GO:0006531, GO:0006535), and transcription factor (IPR007219, IPR001138), were enriched in both H and HH compared to C. The terms related to cytochrome P450 E-class group I (GO:0016491, GO:0055114, IPR001128, IPR023753, IPR002401, IPR017972) were enriched in both groups of salt-resistant fungi compared to the control group.

**Figure 3.** Heat map of the −*log*(*f dr*) of the significantly-enriched or -depleted functional annotations in *A. salisburgensis* compared to *A. sclerotialis*, *Aspergillus ruber*, and *Wallemia ichthyophaga*. Depletion/enrichment found in HH vs. H, HH vs. C, and H vs. C is also shown. Only entries with an enrichment/depletion with an fdr < 1e−5 are shown. (**A**) Depletion. (**B**) Enrichment.

#### *3.4. Amino Acid Composition*

The protein amino-acid (AA) composition of the C, H, and HH groups for the set of orthologous genes and for the set of exported proteins, i.e., containing a signal-P annotation, was compared (Figure 4 and Supplementary Table S1). With respect to the control group, the conserved proteins in halophilic fungi were significantly enriched in glycine, proline, and arginine and depleted in isoleucine and lysine (Figure 4A). The same bias pattern for those 5 AA was previously reported in the extremely halophilic bacteria *Salinibacter ruber* and *Halomonas elongata* and in the Archaea *Halobacterium salinarum* and *Haloarcula marismortui*, compared to *Escherichia coli* [9]. A low occurrence of lysine and increase of arginine in halophilic proteins is a signature of halophilic microorganisms [65,66]. Psychrophilic organisms, which are similar to halophiles, face reduced water availability and exhibit an increased glycine protein content [67].

Among the set of exported proteins, the acidic AA aspartate and glutamate were overrepresented in halophiles compared to the control group. This is a general characteristic found in many osmotolerant bacteria [65,68]. Similarly, in the three studied HH species, the conserved proteins were enriched in glycine and proline. Finally, the serine depletion was previously reported in halophilic Archaea [66].

**Figure 4.** Bias in the amino-acid (AA) distribution between the sets of control, halotolerant, and halophilic fungi. (**A**) Boxplot representing the distribution of AA for the 5 AA with the most significant bias between the control and halophilic group for the set of conserved (top) and secreted proteins (bottom). (**B**) Scatter plot of the −*log*(*p*-value) computed with the Wilcoxon test for the AA distribution bias for the set of conserved proteins and the set of secreted proteins.

#### *3.5. Differential Expression between 5*% *and 20*% *NaCl*

The differential expression of genes between the low and high salt concentrations was studied. Upregulated genes are defined as the set of transcripts that exhibit an increase in expression at 20% salt concentration, while downregulated genes exhibit a reduced expression under high salinity.

#### 3.5.1. *Aspergillus Sclerotialis*

*A. sclerotialis* differentially regulated 2097 genes (1121 up- and 976 downregulated fdr < 0.05, |*log*2*FC*| > 1) between both salt concentrations. The most strongly upregulated gene (x3536) belonged to the major facilitator superfamily (MFS) transporter, which are membrane transport proteins that facilitate movement of small solutes across cell membranes in response to osmotic gradients [69]. Among the 96 genes upregulated by a factor larger than 50 in the halotolerant fungus, 14 genes were related to transmembrane transport (Supplementary Table S1). Cerato-ulmin hydrophobin, a parasitic fitness factor of the agents of Dutch elm disease [70], was upregulated by a factor of 2641. The regulation of this hydrophobin might indicate a role in salt resistance, as previously reported in *Wallemia ichthyophaga* [15]. Given the fact that the halotolerant strain is a dog pathogen, this hydrophobin might be involved in the pathogenicity of the strain [71].

The C-terminal dimerization domain found in transposases of elements belonging to the activator superfamily was increased by a factor 280 under high salt conditions, indicating that DNA modification might take place. A similar feature was seen in sunflower exposed to salt and drought stress [72]. Transcriptional regulators are also found among the most upregulated genes [73]. A gene involved in protein neddylation [74] had its expression increased by a factor of 1571. Inositol monophosphatase, which is upregulated by a factor of 571, is involved in the phosphatidylinositol signaling pathway and was shown to increase Na+ resistance in yeast [75]. An F-box domain-containing protein was upregulated by a factor of 286 upon salt-induced stress, which is in line with previous reports in *Schizosaccharomyces pombe* [76]. Interestingly, two genes involved in the synthesis of antibiotics, Acyl-CoA 6-aminopenicillanic-acid-acyltransferase and pristinamycin IIA synthase subunit A, were upregulated by a factor of 369 and 338, respectively.

Among the genes downregulated by a factor larger than 50, eight were related to transporters, and one was related to SAM-dependent methyltransferase. Chorismatase-degrading enzymes were downregulated, in agreement with the salinity stress experiment in *Carthamus tinctorius* [77], where chorismatase synthase was upregulated upon salinity stress, indicating that chorismate plays a role in salt resistance. Two cupredoxin were downregulated by a factor of 247. A reverse transcriptase and a integrase were downregulated 95- and 129-times, respectively, while one Zn2-Cys6 transcription factor was downregulated by a factor of 124.

The functional enrichment analysis for the downregulated genes indicates that transcriptional factors, MFS transporters, and urease activity (Enzyme <IP_ADDRESS> fdr = 0.0083, GO:0009039 fdr = 3.15e−4) were negatively impacted by the high salt concentration (Supplementary Table S1 and Figure 5A). It was previously shown that a depletion of urease improves salt resistance in *Arabidopsis thaliana* [78].

For the set of upregulated genes (Figure 5B), an enrichment in translation, glycine metabolism, oxidation-reduction process, mitochondrial ATP synthesis, and microtubule was seen. Osmotic stress was previously shown to increase the free glycine concentration in wheat transiently [79]. The five upregulated genes related to microtubules were two tubulins (jg1358, jg9339), one kinesin (jg11911), one dynein (jg2827), and one cytoskeleton-associated protein (jg228). These genes are involved in the transportation of vacuoles and organelles along the microtubules [80].

#### 3.5.2. *A. salisburgensis*

In *A. salisburgensis*, 305 and 306 gene were up- and down-regulated, respectively. The most upregulated gene (x9013) was a peptidase S8 subtilisin protein with putative keratinolytic activity. Although the biological functions of most fungal subtilases are not yet described, some attempts led to the assumption that subtilisins (S8) could be related to a saprotrophic lifestyle [81]. More recently, a study showed striking correlations for subtilisins' (S8) expansion in pathogenic and soil-/dung-inhabiting fungi [82]. Subtilases have further been reported to be involved in drought and salt resistance mechanisms in plants [83]. An example is the *Arabidopsis thaliana* subtilase ATSBT6.1, which is associated with the unfolded protein response on salt stress through the cleavage of an ER-resident type II membrane protein (BZIP28). The cleavage of this protein is essential for the activation of genes associated with the salt stress response [84].

Among the proteins with at least 50-fold upregulation, eight proteins were glycoside hydrolases, three were involved in oxido-reduction processes, four exhibited a cupin-like domain, and two were involved in post-translational modification (Supplementary Table S1). Among the set of genes

downregulated at least 50-fold, four ribitol dehydrogenases, seven transporters, two cyclic amino-acid related enzymes, two transcription factors, and two alcohol dehydrogenase were found.

**Figure 5.** Bar plots showing overrepresented functional annotations in the set of regulated genes in *A. sclerotialis*. (**A**) Enriched functional terms in the set of downregulated genes. (**B**) Enriched functional terms in the set of upregulated genes.

Functional enrichment analyses of the downregulated genes indicated that the expression of membrane-located proteins, especially serine-/threonine-rich proteins, were downregulated (Figure 6A and Supplementary Table S1). Similarly, the gene coding for glutathione S-transferase, which was previously shown to play a negative role in salt stress tolerance in *A. thaliana* [85], was downregulated under high salinity condition. Chitin synthase genes were downregulated, which is in line with proteomic results in yeast [86].

The set of upregulated genes covered a broader spectrum of functions than that of the downregulated genes. Genes located at the cell periphery and extracellular region were overrepresented in the set of upregulated genes (Figure 6B and Supplementary Table S1). The over-representation of functional terms related to cell wall production like the D-alanine–D-alanine ligase and mannose-6-phosphate isomerase might indicate that cell wall damages are happening under high salinity [87]. Enrichment of superoxide-dismutase (K04564), which is an enzyme involved in the catalysis of superoxide, is a known marker of oxidative and salt stress [88]. Finally, the over-representation of MFS transporters is a paramount osmotic stress response, as these transporters are transporting small molecules in response to chemiosmotic ion gradients [69].

**Figure 6.** Bar plots showing overrepresented functional annotations in the set of regulated genes in *A. salisburgensis*. (**A**) Enriched functional terms in the set of downregulated genes. (**B**) Enriched functional terms in the set of upregulated genes.

#### 3.5.3. Comparative Transcriptomics

Difference and similarities in expression of genes found in both studied fungal strains were assessed. A total of 84 genes was one-to-one homologous in both strains and was significantly regulated. Genes upregulated in the halophilic and halotolerant species were related to polyketide synthase, chloroperoxidase, haloacid dehalogenase, fructose-biphosphate aldolase, copper transport, and pectinolytic enzyme (Supplementary Table S1). Chloroperoxidases were shown to be potential chlorinators of lignin in plant materials and to contribute to lignin degradation [89]. This is similar to the use of chlorine for wood pulp delignification [90]. In the fungi *Fusarium fujikuroi* and *Neurospora crassa*, genes belonging to the haloacid dehalogenase family were shown to be involved in the osmotic stress and glycerol metabolism [91]. Conserved genes downregulated in both species are related to amino acid transport, trehalose-phosphatase, thiolase, and the p53 transcription factor.

Genes solely upregulated in *A. sclerotialis* were related to oxido-reductive processes, zinc transport, and polyketide synthase. In contrast, genes upregulated in the halophilic fungus and downregulated in the halotolerant fungus were related to the cation/H+ exchanger, threonine/serine exporter, voltage-dependent channel, intradiol cleavage, the putative sensor protein SUR7/RIM9, which is a pH-response regulator protein, and the isochorismatase, a gene involved in salicylic acid metabolism. The upregulation of the last gene only in the halophilic strain may indicate an adaptive advantage to survive permanently in a hypersaline environment. Salicylic acid (SA) is a natural phenolic

compound known to control many physiological and biochemical functions in plants such as growth, development, and responses to abiotic stresses [92]. Some studies have shown that SA plays a role in the response to salinity in plants and that its exogenous application improves tolerance to salt stress in several species [93].

#### *3.6. Comparison to Fungal Osmoadaptation and Osmoregulation*

The literature on osmoadaptation and osmoregulation in fungi was compared to the genome content and transcription patterns of *A. sclerotialis* and *A. salisburgensis*.

#### 3.6.1. Osmosensing

All members of both branches of the high osmolarity signaling pathway were found in *A. salisburgensis* and *A. sclerotialis* [94,95]. Histidine kinase 7A/B [96] from *Hortaea werneckii*, the mitogen activated protein kinase (KSS1) [97] from *Saccharomyces cerevisiae*, both involved in signal transduction, and the cell wall integrity sensor (MID2) [98] were found in the halotolerant, but not in the halophile genome. Interestingly, the irritation of the cell wall and membrane by shrinking triggers induction of the HOG signaling pathway (e.g., in *S. cerevisiae* [99,100]). Inversely, *A. salisburgensis* has a homologue of the tyrosine-protein phosphatase 3 (PTP3), which is involved in the repression of the *HOG1* gene [101] (Supplementary Table S1).

Protein kinase A (PKA1, TPK1), an enzyme generally involved in the regulation of the glycogen, sugar, and lipid metabolism that negatively regulates the transcription factors (MSN2/MSN4) [102], was upregulated six-fold in *A. salisburgensis*, but downregulated by a factor of three in *A. sclerotialis*. *TOR1*, another gene that negatively regulates the MSN2/MSN4 transcription factors [102], was also downregulated by a factor five in *A. sclerotialis*, but not in *A. salisburgensis*. The PBS2 MAP kinase kinase, which is involved in the HOG pathway, and GSP2, a histidine kinase involved in the nuclear import of HOG1 [103], were upregulated 5.5- and 4.8-fold at a high salt concentration in *A. sclerotialis*, but not in *A. salisburgensis*. NIK1, a histidine kinase involved in the SLN1 branch of the HOG pathway [104], was downregulated nine-fold in *A. sclerotialis*.

#### 3.6.2. Ion Homeostasis

The genome of *A. salisburgensis* contains most of the ion transporters reported to be involved in ion trafficking during osmotic stress in *Saccharomyces cerevisiae* and *Hortaea werneckii*, with the exception of the high affinity potassium transporter (HAK1) [104], the potassium antiporter (KHA1) [15], and the outward-rectifier potassium channel TOK1 [15]. In *A. sclerotialis*, the endosomal/prevacuolar sodium/hydrogen exchanger (NHX1) and HAK1 are missing. TOK1 and the mitochondrial exchanger system, while present in the halotolerant strain, were downregulated four-fold at 20% NaCl in this strain. PAM1, a proton ATPase involved in salt resistance in *Debaryomyces hansenii* [105], was upregulated by a factor of 13 in *A. sclerotialis*. Similarly, NHA1, the Na+/H+ antiporter, was upregulated 6.3-fold under high salt concentration. In the obligate halophile, ENA1, a P-type ATPase sodium pump involved in Na+ and Li+ efflux [106] was upregulated by a factor of 39.

#### 3.6.3. Cellular Respiration

In *A. sclerotialis*, seven genes annotated as involved in respiration (GO:0045333) were upregulated under high salt concentration. The electron carrier protein CYC1 was upregulated by a factor of 4.62; the cytochrome BC1 complex subunit 7 (QCR7) expression increased by a factor of 3.94; and the cytochrome c oxidase subunit VIa (COX13) was induced by a factor of 3.63. Further, an ADP/ATP carrier protein (AAC) was upregulated by a factor of 6.14. Mitochondrial isocitrate dehydrogenase (IDH1), malate dehydrogenase (MDH1), and mitochondrial trans-2-enoyl-CoA reductase were upregulated 4.82-, 12.82-, and 4.62-fold, respectively. This strongly indicates that the ATP production is increased during salt stress in the halotolerant strain. In contrast, no genes related to cellular respiration were significantly regulated in *A. salisburgensis*.

#### 3.6.4. Stress Response

*A. sclerotialis* had a stress response more pronounced than *A. salisburgensis*, both in terms of the number of regulated genes, as well as in the regulation intensity. Markers of oxidative stress response were strongly upregulated in *A. sclerotialis* at a high salt concentration, which is in line with the increased respiratory process seen in the halotolerant strain. The mitochondrial superoxide dismutase (SOD2), an important antioxidant defense, catalyzes the dismutation of superoxide O− 2 , a mitochondrial byproduct of respiration, into oxygen and hydrogen peroxide. This gene was upregulated by a factor of 8.3 in *A. sclerotialis* at high salinity. H2O2 is then reduced to H2O by peroxiredoxin (PRX1) using electrons from the reduced form of thioredoxin (TRX) [107]. In *A. salisburgensis*, TRX was upregulated by a factor of 3.6 at a 20% salt concentration. Two mitochondrial PRX1 were upregulated by 9.5- and 6.5-fold, respectively. Further, the nuclear thioredoxin peroxidase DOT5 and the AHP1 peroxiredoxin were upregulated by a factor 4.8 and 5.9, respectively. H2O2 can also be scavenged through the glutathione system with the help of glutathione redoxins (GRX) [107]. A homologue to GRX1 and GRX2 and a homologue to GRX3 and GRX4 were upregulated 6.31 and 4.3 in *A. sclerotialis*, respectively. Two catalase paralogues (CTT1), an enzyme that is involved in peroxide scavenging, were upregulated by a factor of 10.1 and 6.9, respectively. The homologue to MCA1, a regulator of apoptosis upon H2O2 and clearance of insoluble protein aggregates, was upregulated by a factor of 3.71 [108,109]. Additional studies have shown that oxidative stress originating from intensive mitochondrial respiration [110] can pose a further threat to the survival of aerobic microorganisms living in high-salinity environments. Therefore, this stress can be one of the limiting factors for growth in such environments. In the halotolerant strain *Hortaea werneckii*, the levels of ROS degradation and resistance determine the upper limit of salt tolerance [111]. The peptidyl-prolyl cis/trans-isomerase, a gene involved in protein folding [112], underwent an upregulation under salt stress by a factor of 3.43 in the halotolerant strain.

In *A. salisburgensis*, the oxidative stress response was moderate. Similar to the halotolerant strain, CTT1 was upregulated by a factor of 8.8. Methionine sulfoxide reductase (MXR2), a gene protecting against methionine oxidation by catalyzing thio-dependent reduction of oxidized methionine residues [113], saw a 5.29-fold increase in transcription. Interestingly, the obligate halophile reduced the production of superoxide radicals by downregulating NADPH oxidase and the NADPH oxidase regulator by a factor five and 10, respectively. Beyond genes involved in oxidative stress response, *A. salisburgensis* upregulated HSP12, a chaperone regulated in osmotic stress conditions in *Saccharomyces cerevisiae*, by a factor of five [114].

#### 3.6.5. Cell Interface

In the extracellular region, *A. salisburgensis* strongly upregulated the subtilisin peptidase (S8 x8964), two cellulose degradation enzymes (GH51, x62), a pectinolytic enzyme (EC:<IP_ADDRESS> x23), and a cellulase (GH5 x22.62). Together with the upregulation of a chloroperoxidase, which might be involved in the degradation of lignin, and the environment where this fungus has been isolated, i.e., a wooden staircase in a salt mine, it seems reasonable to assume that salt triggers a cellulolytic response in the obligate halophile. *A. sclerotialis* strongly upregulated a hydrophobin (PF06766 x2646).

In the cell wall of the halotolerant fungus, a beta glucosidase involved in cell wall remodeling (UTH1) [115] was upregulated by a factor of 17, while a membrane-bound HSP70 involved in selective cation trafficking was upregulated by a factor of 6.33 [116]. In the halophilic fungus, the expression of a lysophospholipase increased 18-fold. This enzyme is involved in phospholipid degradation and was previously reported to be upregulated in the micro-algae *Dunaliella salina* under high salt concentration [117].

At the membrane level, besides the previously discussed transporters, *A. sclerotialis* increased the expression of a putative ammonia transporter (PKC1 x22.31), a dipeptidase (DPE1, x9.91), the SUR7 membrane protein involved in membrane organization and cell wall stress [118] (x7.78), a sulfate permease (SULP x6.77), and a cellulolytic enzyme (EggNog: 0PG84 x11.15). The halophile strain increased the transcription of two guanine nucleotide proteins (G protein) [119] (ion-translocating rhodopsin x10.63, RHEB x4.00) and one G protein-coupled receptor-like GPR1 (x12.46), a gene involved in endocytosis, the eisosome component PIL1/LSP1 (x32.11) [120], and two SUR7 membrane proteins (x4.37, x25.73).

#### 3.6.6. Compatible Solute Management

Known genes involved in compatible solutes' management from halotolerant and obligate halophilic fungi were reviewed. Compatible solutes are known to be accumulated or synthesized during changed osmotic conditions. One of these is D-mannitol, which can be synthesized by a reduction of fructose, catalyzed by NADP-dependent mannitol dehydrogenases [13]. In *A. salisburgensis*, two homologues of NADP-dependent mannitol dehydrogenases were found in the genome, but were not regulated in high salt concentration, while one homologue was found in *A. sclerotialis* and was upregulated by a factor of 29. Still, *A. salisburgensis* upregulated GRE3, an aldole reductase involved in polyol metabolism [114], by a factor of 28.

Three homologous genes of STL1, a glycerol/proton symporter of the plasma membrane, were found in *A. salisburgensis* and one in *A. sclerotialis*. In comparison, four homologues were found in the genome of the obligate halophilic fungus *Wallemia ichthyophaga* [15]. In *A. salisburgensis*, one copy of the genes was upregulated by a factor of eight, and in *A. sclerotialis*, STL1 was upregulated by a factor of 91. In *Saccharomyces cerevisiae*, it was reported to be upregulated during hyperosmotic shock [121]. The aquaglyceroporin channel FPS1, which is opened in hypoosmotic conditions to release glycerol [122], did not receive any regulation in *A. salisburgensis* and was downregulated in *A. sclerotialis* by a factor of 30.

Glycine betaine is an osmoprotectant found in plant, animals, bacteria, and fungi and is involved in reactive oxygen scavenging [123]. In *Aspergillus fumigatus*, glycine betaine is produced by converting first choline to betaine aldehyde with a choline oxidase. Betaine aldehyde or its rapidly-forming hydrated equivalent gem-diol-choline are then converted to glycine betaine with the help of betaine aldehyde dehydrogenase or choline oxidase, respectively [123]. In *A. salisburgensis*, two choline oxidases and one betaine aldehyde dehydrogenase were found, while in *A. sclerotialis*, two copies and one copy were found, respectively. While *A. salisburgensis* did not significantly regulate these enzymes, *A. sclerotialis* increased the transcription of betaine aldehyde dehydrogenase by a factor of 4.5 in the high salinity condition.

Genes involved in polyols metabolism were also strongly upregulated. Polyols (also called sugar alcohols) compensate for differences between the extracellular and intracellular water potential without affecting the integrity and function of proteins [124]. In *A. sclerotialis*, gene jg11062 belongs to the same group as the L-arabitol 4-dehydrogenase from *Aspergillus fumigatus* and was upregulated by a factor of 1063 at high salt concentration (Supplementary Table S1). A sorbitol dehydrogenase and a ribitol dehydrogenase were upregulated by a factor of 62 and 67, respectively. The important role of polyols for the resistance of extremotolerant fungi under osmotic stress was already described by Sterflinger [125]. However, this general protection strategy is widespread also in bacteria and Archaea [126].

In the halotolerant strain, four SAM-dependent methyltransferases showed an increase in transcription by a factor of 455, 218, 205, and 52, respectively, similar to previous reports in algae [127], where sequence homologues of the SAM-dependent methyltransferase are involved in the production of the compatible solute homoserine betaine [128]. Glutamine synthetase was upregulated by a factor of 181, indicating that *A. sclerotialis* might use glutamine as an osmolyte, something previously reported for *Halobacillus halophilus* [129].

#### **4. Conclusions**

In this work, the genomic adaptation and gene regulation of the obligate halophile *A. salisburgensis*, a fungus isolated from a wooden staircase in a salt mine, were studied by employing comparative genomics and transcriptomics approaches. The halotolerant relative *A. sclerotialis*, a pathogen, was used to gain insight into the environment-specific osmoadaptation and regulation of the obligate halophile.

On the genomic level, *A. salisburgensis* exhibited a 27% decrease in genome size and gene content compared to *A. sclerotialis*. Considering the unique, extremely stable niche where the obligate halophile has been found, it can be hypothesized that the fungus optimized its genome content by dumping genes unnecessary for its survival, a known phenomenon in bacteria [130]. Niche adaptation is further seen in the enrichment of genes involved in cellulose degradation. Halophile-specific depletion of MFS transporters previously reported in *Wallemia ichthyophaga* and *Aspergillus ruber* [10,15] was also found in *A. salisburgensis* and might therefore be a strategy to survive in the high saline environment. This study further confirmed the specific amino-acid enrichment and depletion patterns found in other halophilic species compared to the control species.

The fact that the obligate halophile regulated three-times fewer genes than the halotolerant strain between both salinities further underlines the adaptation of the former to high salt concentration. Besides the regulation of a few transporters like ENA1, STL1, a hydrophobin, and an aldol-reductase, a gene reported to be involved in the production of compatible solutes, there was almost no sign revealing that the obligate halophile was under stress at a 20% salt concentration. Instead, the cellulolytic activity was an indication that the high salt concentration was beneficial to *A. salisburgensis*, as it triggered the expression of a battery of wood-degrading enzymes. Among them, the chloroperoxidase was making use of the chloride ion to chlorinate lignin and improved its degradation.

In contrast, the halotolerant strain exhibited both an osmotic- and oxidative-stress response. The link between both stresses might probably be found in the respiration processes. In fact, it is probable that under high salinity, homeostatic regulation requires an increased supply of ATP, which is produced mainly from respiration. At the transcriptome level, the increased respiration is seen in the upregulation of genes involved in the electron transport chain and the citrate cycle. ROS originating from the respiration induced an oxidative stress response, as can be seen from the increased transcription of genes involved in oxidative-stress, such as superoxide dismutase, thioredoxin peroxidases, glutathione oxidoreductases, and catalases. *A. sclerotialis* further upregulated a hydrophobin, a family of proteins that plays a role in salt resistance in *Wallemia ichthyophaga* and in pathogenicity.

**Supplementary Materials:** The supplementary table and the supplementary figure are available online at http://www.mdpi.com/2073-4425/10/5/381/s1.

**Author Contributions:** The conceptualization was done by G.P. and K.S.; methods were planned by G.P., H.T., K.L. and C.P.; Labwork was done by K.L., C.P. and G.P.; Bioinformatics and Data analysis was done by H.T.; original draft preparation and data interpretation was done by H.T; Editing, Data interpretation and reviewing was done by C.P., K.S., K.L. and G.P.; Funding aquisition was done by K.S.

**Funding:** This research was funded by BOKU-VIBT Equipment GmbH.

**Acknowledgments:** The computational results presented have been achieved using the Vienna Scientific Cluster (VSC).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **The Primary Antisense Transcriptome of** *Halobacterium salinarum* **NRC-1**

**João Paulo Pereira de Almeida 1,**†**, Ricardo Z. N. Vêncio 2,**†**, Alan P. R. Lorenzetti 1, Felipe ten-Caten 1, José Vicente Gomes-Filho 1 and Tie Koide 1,\***

Received: 6 February 2019; Accepted: 1 April 2019; Published: 5 April 2019

**Abstract:** Antisense RNAs (asRNAs) are present in diverse organisms and play important roles in gene regulation. In this work, we mapped the primary antisense transcriptome in the halophilic archaeon *Halobacterium salinarum* NRC-1. By reanalyzing publicly available data, we mapped antisense transcription start sites (aTSSs) and inferred the probable 3 ends of these transcripts. We analyzed the resulting asRNAs according to the size, location, function of genes on the opposite strand, expression levels and conservation. We show that at least 21% of the genes contain asRNAs in *H. salinarum.* Most of these asRNAs are expressed at low levels. They are located antisense to genes related to distinctive characteristics of *H. salinarum*, such as bacteriorhodopsin, gas vesicles, transposases and other important biological processes such as translation. We provide evidence to support asRNAs in type II toxin–antitoxin systems in archaea. We also analyzed public Ribosome profiling (Ribo-seq) data and found that ~10% of the asRNAs are ribosome-associated non-coding RNAs (rancRNAs), with asRNAs from transposases overrepresented. Using a comparative transcriptomics approach, we found that ~19% of the asRNAs annotated in *H. salinarum* belong to genes with an ortholog in *Haloferax volcanii*, in which an aTSS could be identified with positional equivalence. This shows that most asRNAs are not conserved between these halophilic archaea.

**Keywords:** antisense RNA; *Halobacterium salinarum*; transcriptome; dRNA-seq; archaea; transcription start site; post-transcriptional regulation; gene expression; type II toxin-antitoxin systems; Ribo-seq

#### **1. Introduction**

Antisense RNAs (asRNAs) are non-coding RNAs (ncRNAs) transcribed from the opposite strand of a given gene. Intuitively, asRNAs can be assumed to be cis-acting as they are complementary to the messenger RNA (mRNA) of the gene from which they derive. This does not restrict asRNA action to the gene on the opposite strand, they can also act in trans [1]. AsRNAs can act as regulators at different stages of gene expression [2]. They can modulate the stability and lifespan of RNAs by either occluding degradation sites or forming double-stranded complexes (dsRNAs) that are targets for RNases [3]. Furthermore, they may directly influence translation, inhibiting this process by occluding ribosome binding sites or promoting conformational alterations that might increase or decrease the frequency of translation of its target [2,4,5].

Simultaneous transcription of opposite DNA strands was first reported in λ phage [6], later it was identified in bacteria in processes related to plasmid replication, phage repression and transposases [7], and finally in eukaryotes [8]. Technical advances in transcriptome analysis allowed greater confidence in strand-specific datasets, allowing the identification of asRNAs in different organisms [9–12]. Overlapping antisense transcription is a ubiquitous phenomenon [13–15], while the functionality of these molecules is still under debate, since it can arise from spurious transcription [16]. In bacteria the percentage of protein coding genes with asRNAs varies from 2% in *Salmonella* [17] to up to 80% in *Pseudomonas aeruginosa* [18]. In archaea, asRNAs have been reported for a while [19,20]. Some studies assured that they were primary transcripts: 1244 in *Haloferax volcanii* [21], 1110 in *Methanolobus psychrophilus* [22], 1018 in *Thermococcus kodakarensis* [23], 48 in *Methanosarcina mazei* Gö1 [24], 29 in *Thermococcus onnurineus* NA1 [25], and 12 in *Methanocaldococcus jannaschii* DSM [26]. In *H. volcanii*, the only halophile represented, asRNAs are present in genes related to chemotaxis, transcription regulation, and insertion sequences [27].

*Halobacterium salinarum* is an extreme halophile archaeon, thriving in 4.3M NaCl concentrations [28]. Historically known for its bacteriorhodopsin light-dependent proton pump [29], *H. salinarum* has a wide metabolic versatility for energy production such as amino acid oxidation, anaerobic respiration using different electron acceptors, and arginine fermentation [30,31]; it also shows high tolerance to diverse environmental stresses [32–34]. The plethora of transcriptional data in different environmental conditions and genetic backgrounds allowed the construction of a global gene regulatory network based primarily on mRNA levels [35,36]. Many overlapping features in *H. salinarum* genome, such as alternative promoter usage in operons [37] and internal RNAs that produce protein isoforms [38], have been reported, highlighting the complexity of transcriptional regulation. In addition, ncRNAs in intergenic regions [37], RNAs associated to transcription start sites (TSSaRNAs) [39], and associated to insertion sequences (sotRNAs) [40] have also been identified. However, there is little information on antisense RNAs. The first asRNA reported in this organism was a 151 nucleotide (nt) molecule complementary to the 5 end of a phage related gene [41,42]; asRNAs regulating gas vesicles genes were also identified [19], but no further global analysis of asRNAs have been performed since then.

In this work, we report a genome-wide primary asRNA mapping in *H. salinarum* NRC-1. To do that, we took advantage of published differential RNA-sequencing (dRNA-seq) data for the identification of transcription start sites antisense to annotated genes (aTSS) [38]. To annotate the minimum length of asRNAs, we reanalyzed published strand-specific RNA-seq data acquired in paired-end mode. Data was manually inspected to guide asRNA annotation of probable 3 ends. Publicly available Ribosome profiling (Ribo-seq) data (Raman et al. submitted) were reanalyzed to identify asRNAs bound to ribosomes. We investigated asRNA global properties: size distribution, location relative to the cognate genes, expression levels, and cognate gene functions. Finally, we did a comparison with dRNA-seq results of *H. volcanii* to evaluate the cross-species conservation of asRNAs.

#### **2. Materials and Methods**

#### *2.1. Antisense Transcription Start Sites Annotation*

We reanalyzed publicly available raw dRNA-seq data from *H. salinarum* NRC-1 grown in complex media (250 g/L NaCl, 20 g/L MgSO4,2g/L KCl, 3 g/L sodium citrate, and 10 g/L bacteriological peptone (Oxoid)) over a growth curve sampled at 17 h, 37 h, and 86 h, and grown under standard conditions (37 ◦C, 225 rpm, constant light) sampled at mid-log phase [38]. We used our in-house workflow ("Caloi-seq", https://github.com/alanlorenzetti/frtc/), described in [38], for *H. salinarum* dRNA-seq data. Briefly, libraries were downloaded from NCBI Sequence Read Archive (SRA) [43] and trimmed using Trimmomatic [44]. Reads surviving as a pair were aligned to reference genome (NCBI Assembly ASM680v1) in paired-end mode using HISAT2 [45], suppressing alignments resulting in fragments longer than 1000 nt. Orphan R1 and R2 sequences were aligned using the single-end mode. Multimappers aligning up to 1000 times were allowed to be reported. SAM files were converted to BAM using SAMtools [46] and input in MMR to find the most likely position for each multimapper [47]. The resultant BAM files were filtered to keep only R1 reads, which serve as input for the TSS inference

step. BEDTools [48] was used to compute the fragment 5 accumulation profile, employing all the aligned R1 reads. Data visualization was performed using IGV [49] and Gaggle Genome Browser [50].

TSSs were identified from dRNA-Seq experiments using TSSAR [51] with the following parameters: *p*-value *p*<0.005, a minimum of four reads, and a distance of TSS grouping of at least five nt. An antisense TSS (aTSS) was defined based on genome annotation of *H. salinarum* NRC-1 available at RefSeq updated in 2017 with additional manual curation using *H. salinarum* R1 annotation as reference [52] (Table S1). TSSs located inside genes on the opposite strand were classified as aTSS. TSSs on the opposite strand and up to 200 nt downstream of a gene 3 end were considered downstream aTSS (daTSS).

Potentially structured regions were filtered out by calculating folding minimum free energy (MFE) along the whole genome sequence as previously described [38]. A sliding window (51 nt) with an offset of 10 nt was used to tile the genome and all subsequences were subjected to secondary structure prediction using RNAfold [53] with default parameters. The distribution of MFE obtained for the tiled genome was compared with the distribution obtained for only subsequences immediately downstream of aTSS. The 33.3% quantile in the whole-genome MFE distribution was arbitrarily chosen as cutoff for potentially forming structures and thus putative false positive aTSS.

#### *2.2. Antisense RNA Loci Inference*

To annotate asRNA loci we used aTSS positions together with all TEX- strand-specific paired-end RNA-seq libraries publicly available for *H. salinarum* NRC-1. Library accessions and number of sequenced, trimmed, and aligned reads used for asRNA loci annotation are shown in Table S2. These libraries are the control libraries for a dRNA-seq experiment, representing four replicates of *H. salinarum* NRC-1 grown under standard conditions and sampled at mid-log phase and biological duplicates of three different time points over a growth curve (17 h, 37 h, and 86 h) [38]. For these libraries, we took the MMR adjusted BAM files, cited in the last section, and computed the genome-wide coverage (transcriptional profile) using deepTools [54], considering the extension of full fragments for paired-end alignments and the proper strand orientation for alignments of R1 and R2 orphan reads. BEDTools [48] was used to compute the fragment 3 accumulation profile, employing all the aligned R2 reads. Visual inspection was performed using IGV [49] and Gaggle Genome Browser [50] by looking for a mapped aTSS followed by a sharp decrease in read coverage to infer an asRNA locus in at least four different libraries. To aid the inference of locus ending, we used the fragment 3 end accumulation profile, requiring at least 10 observations to demarcate the minimum end position of an asRNA. The steps for processing dRNA-seq data and computing the transcriptional profiles are depicted as a workflow in Figure S1, and asRNA loci inference is schematically shown in Figure S2.

#### *2.3. Promoter and 3 End Sequence Analysis*

The computation of nucleotide frequency for the promoter region of aTSSs and 3 end of asRNAs was performed using Weblogo [55]. For the promoter region we analyzed sequences 40 nt upstream and 10 nt downstream of the aTSS; for the 3 end, we used sequences 10 nt upstream and downstream of the last nucleotide.

#### *2.4. Gene Functions*

Gene functional classification and Gene Ontology enrichment analysis were performed using the PANTHER system [56]. In the absence of an associated PANTHER annotation, MicrobesOnLine [57] and BlastCD [58] tools were used to identify conserved protein domains.

#### *2.5. Type II Toxin-Antitoxin Systems Annotation and Antisense Transcription Start Sites Identification in Other Archaea*

TA finder 2.0 [59] was used to annotate type II Toxin-Antitoxin (TA) systems in *E. coli* and archaeal organisms with available dRNA-seq data. aTSS positions were mapped as described in Section 2.1, reanalyzing dRNA-seq data from NCBI BioProject database for *Haloferax volcanii* DS2 (PRJNA324298), *Methanocaldococcus jannaschii* DSM 2661 (PRJNA342613), *Thermococcus kodakarensis* KOD1 (PRJNA242777), and *Thermococcus onnurineus* NA1 (PRJNA339284). aTSS positions for *E. coli* were obtained from Thomason et al. [60], and positions of asRNAs in type II TA systems found in an immunoprecipitation ofdsRNAs experiment were obtained from Lybecker et al. [61].

#### *2.6. Di*ff*erential Expression Analysis*

A GFF file including *H. salinarum* NRC-1 annotated genes, in addition to the 846 asRNAs identified in this study, was built, generating a total of 3680 features. Matrices for read counts per feature were generated using HT-Seq [62] from BAM files. Differential expression analysis was performed using DEseq2 [63]. Genes with more than two-fold change (FC) up or down regulation (log2 FC > 1 or < −1) and adjusted *p*-values *padj* < 0.01 were considered differentially expressed.

#### *2.7. Antisense Transcription Start Sites Comparison between Halobacterium salinarum and Haloferax volcanii*

Ortholog genes in *H. salinarum* and *H. volcanii* were identified using OrtholugeDB 2.1 [64]. Pairs of genes with aTSSs were sub-selected and gene sizes were normalized by their length in a 0 to 100 normalized scale (*d* = 0 at start codon, *d* = 100 at stop codon). The length-normalized genes were divided in 10 partitions (tenths) depending on *d*: [0;10), [10;20), ... , [90;100), and aTSSs were considered conserved if present in the same partition in the orthologous genes. Differences between relative positions, *D* = |*dHsal* − *dHvol*|, were used to estimate how "equivalent" aTSS are between both organisms.

#### *2.8. Ribo-Seq Data Analysis*

*H. salinarum* NRC-1 Ribo-seq data is publicly available ahead of publication (Raman et al. submitted) at NCBI BioProject database under accession PRJNA413990 and was reanalyzed using the same pipeline described in Sections 2.1 and 2.2. AsRNAs presenting coverage greater than 50% and at least 20 reads were considered putative ribosome-associated ncRNAs (rancRNAs).

#### **3. Results**

#### *3.1. Mapping Primary Antisense RNAs*

We identified the primary transcription start sites (TSS) for asRNAs (aTSS) in *H. salinarum* NRC-1 by reanalyzing dRNA-seq data sampled at three time points over growth [38]. We were able to identify 2146 aTSS, located in 1231 genes. Probable false positives due to secondary structure were filtered out (see Methods 2.1), resulting in 1626 aTSSs located in 963 genes (Table S3). Figure 1a shows the position of aTSS in a length-normalized scale (0 to 100 scale) relative to the gene on the opposite strand. The overrepresentation of aTSS at the 3 end is stronger than at the 5 end.

**Figure 1.** Antisense RNA (asRNA) properties. (**a**) Antisense transcription start site (aTSS) positions relative to cognate genes. (**b**) Size distribution of mapped asRNAs.

Given that (i) asRNAs can be located in regions downstream of the cognate gene, (ii) 3 untranslated regions (UTR) can be targeted by ncRNAs, and (iii) asRNAs from 3 UTRs are important regulatory RNAs in archaea with short 5 UTRs [65], we searched for downstream antisense TSS (daTSS). They were defined as transcript start sites located on the opposite strand of annotated genes, located up to 200 nt downstream of the 3 end of a gene. We found 80 daTSS (Table S4).

Having strand-specific paired-end data available, we could search for genome positions at which antisense transcripts preferentially terminate. Such presumed RNA 3 ends, if downstream of an aTSS (RNA 5 end), were inferred to represent the end of an asRNA locus (see Methods 2.2). We were able to define 846 asRNAs loci distributed in 613 genes, starting from aTSS or daTSS (Table S4). Most of the genes contain only one asRNA locus and 26 of the asRNAs are antisense to two genes. The size distribution of asRNAs is shown in Figure 1b, most of them being smaller than 100 nt.

We analyzed the frequency of nucleotides around aTSSs and daTSSs (Figure S3a) and observed the predominance of purines at the TSS, pyrimidines at position −1 and a BRE-TATA signature, showing a characteristic archaeal promoter region for the mapped asRNAs. The dinucleotide composition around the TSS is postulated as important for the transcription initiation among organisms from all domains [66–68], as well as the GC enrichment around the −36 position for TFIIB binding and TATA-box at −26 [69,70]. As previously shown, dRNA-seq data reliably identifies TSS genome-wide in *H. salinarum* NRC-1 [38], recapitulating known features of the promoter region. Nevertheless, this does not exclude the possibility that many of the asRNAs identified are products of spurious transcription.

For asRNAs with at least 10 fragments that corroborated the 3 end annotation, we analyzed the frequency of nucleotides 10 nt upstream and downstream of the 3 end. We observed an enrichment of pyrimidines, specially uracil at positions 0 and −1 (Figure S3b). Comparing this signature with poly-U signatures found in other archaea [71], our data shows a much shorter region, which is unlikely to form a secondary structure involved in transcription termination of asRNAs.

Most of the characterized sRNAs in prokaryotes bind to the 5 region of the mRNA, pairing with the start codon or ribosome binding site (RBS) in the 5 UTR [17], and thus inhibiting the access of the ribosomal machinery. asRNAs that overlap these regions can be candidates to act through a similar mechanism. Halophiles such as *H. salinarum* NRC-1 and *H. volcanii* have mostly leaderless transcripts [21,37,72]. However, it has been experimentally shown in *Salmonella* that ncRNAs overlapping nucleotides up to the 5th codon in the mRNA are capable of blocking the translation machinery, even without AUG or Shine–Dalgarno pairing [73]. We thus looked for asRNAs that overlap their cognate gene s start codon or at least 12 nt (four codons) downstream of it. We found 145 asRNAs in *H. salinarum* overlapping the 5 end of genes (Table S5), which could, in principle, impair

mRNA translation by occluding ribosome binding. Figure 2 and Figure S4 show examples of putative RBS occlusion by an asRNA: *gcvP1*, encoding a glycine dehydrogenase subunit (VNG\_RS06215); *cdcH*, encoding an AAA-type ATPase (VNG\_RS06465); and *rpl1*, which encodes the 50S ribosomal protein L1 (VNG\_RS04315). The gene *gcvP1* is the first from an operon and has a strong Shine–Dalgarno-like signature at −19 upstream of the start codon, colocalized with asRNA VNG\_as06215\_888 (Figure 2). Similarly, RBS at −11 and −13 upstream of *rpl1* and *cdcH* are colocalized with asRNA VNG\_as04315\_654 and VNG\_as06465\_925, respectively (Figure S4).

**Figure 2.** Antisense RNA in 5 untranslated region (UTR). *gcvP1*, encoding a glycine dehydrogenase subunit (VNG\_RS06215). VNG\_RS06215 locus (orange arrow) is in reverse strand (5 →3 right to left), neighbor gene VNG\_RS06210 (orange arrow) is also in reverse strand. Differential RNA-sequencing (dRNA-seq) read coverage signal is shown in dark and light green for TEX+ and TEX− libraries, respectively. Coverage signals below and above the central axis are for reverse and forward strands, respectively. VNG\_as06215\_888 asRNA (pink arrow) encompasses Shine–Dalgarno-like signature (\* light blue highlight in genome coordinates and zoomed in sequence).

All data are available for browsing in Gaggle Genome Browser format (interactive versions of Figure 2 and similar outputs) at http://labpib.fmrp.usp.br/~{}rvencio/asrna/.

#### *3.2. Antisense RNAs Expression Levels*

Some of the asRNAs with characterized functions are expressed at levels equivalent or higher than the mRNA on the opposite strand, which would be expected if a dsRNA is necessary for post-transcriptional regulation [74]. However, most of the asRNAs identified in eukaryotes and prokaryotes are expressed at low levels, which has challenged their identification before high-resolution sequencing methods were available [2,16,60]. Low expression levels might indicate that these asRNAs are products of spurious transcription from low complexity promoter regions [16]. However, this does not exclude the possibility of a functional RNA, given that even at low levels, asRNAs can present a buffering effect and fine-tune gene regulation [2,3].

We observe a negative correlation between asRNA and mRNA transcripts when we compare the average read counts for these transcripts (Figure S5a), as observed in *H. volcanii* [21]. Using RNA-seq data, we analyzed the relationship between the fold change of the read counts of asRNAs relative to the mRNA on the opposite strand (Figure 3). We observed that most of the asRNAs annotated in our study present low expression levels relative to the gene on the opposite strand. Only 112 asRNAs (~13%, Table S6) present expression levels equal or greater than the gene on the opposite strand in at least one of the conditions analyzed. These molecules could be candidates for *cis*-regulators of their respective cognate genes in the considered experimental conditions.

**Figure 3.** Relative expression levels of asRNAs and messenger RNA (mRNAs) on the opposite strand for arbitrarily selected representative libraries: (**a**) stationary phase (17 h) and (**b**) gas vesicle release phase (86 h). Vertical dotted lines mark 1:1 expression levels. Expanded version in Figure S5.

We also evaluated the expression profiles of the asRNAs and cognate genes over the growth curve of *H. salinarum* NRC-1 using RNA-seq libraries sampled at three time points [38]. We compared the expression levels from (i) the stationary phase relative to exponential phase (17 h vs. 37 h) and (ii) late exponential gas vesicle release phase relative to stationary phase (37 h vs. 86 h).

In the first transition, we observed 93 asRNA differentially expressed (Table S7). For 26 asRNAs, the gene on the opposite strand is also differentially expressed (Figure S6a), and six asRNAs overlap the 5 UTR of the genes. In the second transition, we observed 63 asRNAs differentially expressed (Table S8). For 30 cases, the gene on the opposite strand is also differentially expressed. Twenty-seven pairs were differentially expressed in the same direction and only three in opposite directions (Figure S6b). Of them, seven asRNAs overlap the 5 UTR region. Overall, 56 pairs of gene-asRNAs are differentially expressed in at least one of the transitions analyzed. These asRNAs can be potential *cis* regulators of the cognate genes.

#### *3.3. Function of the Genes on the Opposite Strand of Antisense RNAs*

From the total of 613 genes that present at least one asRNA locus, 198 are hypothetical proteins and 32 are nonredundant transposases (i.e., some of them are multicopy, but counted only once). Gene enrichment analysis returned no overrepresented Gene Ontology (GO) term relative to the whole genome background distribution. Genes with asRNAs are associated with different functions according to GO categorization (Figure S7). Next, we present the results for some of the gene categories.

*H. salinarum* NRC-1 is known to produce intracellular gas vesicles: structures composed of proteins that are filled with gas and allow floating on the media surface. Their synthesis and regulation have been studied in details and includes multiple layers of regulation [75]. We were able to identify 11 asRNAs in the *gvp* gene cluster (Table 1), including one previously identified on the opposite strand of *gvpD* [19], which encodes a repressor of gas vesicle production. Most asRNAs related to *gvp* genes showed similar expression levels relative to their cognate genes.

**Table 1.** asRNAs in *gvp* gene cluster, located in pNRC100. (nt: Nucleotide).

AsRNAs to the *gvpD* gene have been detected by Krüger and Pfeifer [19], complementary to the 5 and 3 end of the gene and detected when GvpD protein was present at low levels. We were able to identify three RNAs antisense to *gvpD*, one of them (VNG\_as12315\_1568) overlaps the 5 end that recapitulates previously published information. This asRNA was annotated as a 164 nt molecule in our work, presenting a similar size to the 190 nt band observed in [19]. The aTSS mapped is located 3 nt downstream to the 5 end of the probe used in the previous work, showing the reliability of our asRNA mapping (Figure S8). We also identified primary asRNAs colocalized with genes *gvpC*, *gvpH*, *gvpI*, *gvpJ*, *gvpL*, and a strong signal antisense to the 5 end of *gvpA*, although no aTSS was mapped (Figure S8). These data indicate asRNAs as important players in gas vesicle regulation in *H. salinarum*, as reported in cyanobacteria *Calothrix* sp. PCC 7601 [76].

Moreover, we were able to identify RNAs possibly involved in rhodopsin regulation in *H. salinarum*: asRNAs in *bop* (VNG\_RS05715—bacteriorhodopsin) and its regulators *brz* (VNG\_RS05710—bacteriorhodopsin regulating zinc finger protein) and *brb* (OE3105F bacteriorhodopsin regulating basic protein), in addition to an asRNA in halorhodopsin (VNG\_RS00745) (Table 2). Brb fine-tunes the activation of *bop*; this activity was experimentally shown using reporter genes and mutagenesis. The Brb protein was proven to exist, but was not detected using mass spectrometry [77]. The asRNA overlapping bacteriorhodopsin regulators (VNG\_da3105F\_36) starts downstream of *brb*, overlaps it completely, and ends inside *brz* gene (Figure S9). In the RNA-seq libraries analyzed in this study, the number of reads for this asRNA was approximately 4-fold higher than for *brb* gene. If asRNAs should be present at high levels to post-transcriptionally regulate an mRNA, this data could account for the difficulty in Brb protein detection since asRNA–mRNA pairing could block translation, indicating a possible role for this asRNA. The presence of asRNAs in bacteriorhodopsin and halorhodopsin genes could indicate additional regulators to be studied for understanding the photobiology of this organism.

In bacteria, type I TA systems are known for the presence of an asRNA acting as an antitoxin, while in type II TA systems, both toxin and antitoxin are known to be proteins [78]. In 2014, Lybecker et al. [61] found asRNAs to type II TA systems in *E. coli*, indicating a possible role for asRNAs in these systems. In *H. salinarum* NRC-1, we identified asRNAs to genes VNG\_RS11240 and VNG\_RS00140, which are annotated type II antitoxins.

The presence of asRNAs in type II TA systems has not been systematically explored and could indicate another layer of regulation for these systems. To verify the compatibility of TSS identification and asRNAs identified by Lybecker et al. using dsRNA immunoprecipitation [61], we reanalyzed dRNA-seq data from *E. coli* [60], and we were able to detect the aTSS corresponding to the asRNAs reported by Lybecker et al. [61] in type II TA systems. Given that aTSS identification was reliable in *E. coli*, we looked for asRNAs in type II TA systems in archaea by precisely annotating these genes using TA finder 2.0 [59] and reanalyzing available dRNA-seq data for *H. salinarum* NRC-1 (PRJNA448992) [38], *H. volcanii* DS2 (PRJNA324298) [21], *M. jannaschii* DSM 2661 (PRJNA342613) [26], *T. kodakarensis* KOD1 (PRJNA242777) [23], and *T. onnurineus* NA1 (PRJNA339284) [25]. We annotated new type II TA systems in archaea (Table S9), including nine complete pairs in *H. salinarum*. Both genes composing one of these new pairs, VNG\_RS11890 (toxin) and VNG\_RS11895 (antitoxin), have asRNAs displaying expression levels higher than the cognate genes (Figure S10). Most of the annotated type II TA systems have at least one aTSS in one of the genes. For *T. kodakarensis* and *M. jannaschii*, aTSSs are predominantly in the toxin genes.

We found asRNAs for 37 genes related to translation process, including ribosomal proteins, translation initiation factors, transfer RNA (tRNA) ligases and tRNAs, Asn, Lys, and Ser (Table S10). There are reports of RNAs antisense to tRNAs in *S. solfataricus*[79] and *T. kodakarensis*[23], which indicate a conserved regulatory role for these molecules.

Interestingly, 58 out of 846 asRNAs are overlapping transposase genes in *H. salinarum* NRC-1 (Table S11). Since transposases are usually encoded within repetitive elements called insertion sequences (IS), we eliminated redundancy in numbers by choosing only one representative element for multicopy entities. AsRNAs in transposases have been reported in other archaea such as *H. volcanii*, *T. kodakarensis*, *S. solfataricus*, and *M. mazei* [20]. Retrieving legacy data from tiling microarray experiments performed along *H. salinarum* growth curve [37], we verified that several transcripts antisense to IS are differentially expressed (Figure S11).

In bacteria, there are examples of asRNA inhibiting the translation of transposase mRNAs by occluding the ribosomal machinery assembly at the 5 end of an mRNA [80]. We found 10 asRNAs that overlap the 5 end of transposase coding gene (highlighted in Table S11), which could be potential candidates for a similar regulatory mechanism.

#### *3.4. Ribosome-Associated Antisense RNAs*

We reanalyzed Ribo-seq data, (BioProject PRJNA413990), to identify asRNAs that are putative ribosome-associated ncRNAs (rancRNAs) [81]. By looking at asRNA loci covered by at least 20 reads along at least half its extension, we identified 91 asRNAs (~11%) with relevant signal (Table S12). Recently published rancRNA data in *H. volcanii* found 68 candidates antisense to genes, ~6% of their total [82].

Interestingly, 11 asRNAs with relevant Ribo-seq signal are located in transposases (representative instances highlighted in Table S12). This might indicate that these asRNAs can be either targets for translation or regulate/interfere with the ribosomal machinery. Given that years of *H. salinarum* proteomics studies refuted a widespread colocalization of open reading frames in both strands as spurious "overprediction" [52], the regulatory ribosome binding scenario seems more plausible.

#### *3.5. Conservation of Antisense Transcription Start Sites*

Transcriptome analysis of bacteria has shown that the conservation of asRNAs even among phylogenetically close organisms is low. The comparison between *E. coli* and *Salmonella enterica* serovar Typhimurium showed that only 14% of the asRNAs are conserved [83]. The number of conserved aTSSs varies in different organisms: comparison between *Campylobacter* strains showed 45% conservation [84]; within eight different species of *Shewanella* genus, 22% [85]; and among *Synechocystis* strains, only 4% [86].

*Genes* **2019**, *10*, 280

To evaluate the conservation of aTSS, we compared the identified positions in *H. salinarum* with dRNA-seq data reanalysis of *H. volcanii.* First we sub-selected pairs of orthologous genes in both halophiles (1554 pairs) [64]. Then, from these groups of orthologous genes, we sub-selected pairs with an annotated asRNA in *H. salinarum* s genes and at least one aTSS in its correspondent ortholog in *H. volcanii*, obtaining 244 pairs. We normalized genes sizes by coding sequences (CDS) length defining *d* = 0 at start codons and *d* = 100 at stop codons. We arbitrarily partitioned CDS in 10 regions depending on *d*: (0;10), (10;20), and so on up to (90;100). Then, we searched for annotated asRNAs in *H. salinarum*'s genes and aTSSs in their correspondent ortholog in *H. volcanii* located in the same "equivalent" region. We were able to identify 160 asRNAs, distributed in 110 *H. salinarum*'s genes that contain at least one aTSS in the same tenth region of its ortholog gene in *H. volcanii* (Table S13), representing ~19% of the annotated *H. salinarum* asRNAs. We deemed an aTSS as a conserved feature if its relative position falls into the same tenth partition in both organisms. The distribution of differences between relative positions (*D* = |*dHsal* − *dHvol*|) shows that the majority of such 160 conserved aTSS are in fact positionally equivalent (*D* < 3) (Figure 4).

**Figure 4.** Distribution of differences (*D*) between relative positions among aTSS conserved in *H. salinarum* NRC-1 and *H. volcanii* DS2. Positions are coding sequences (CDS) length-normalized thus *D* = 0 mean same position in both organisms.

We found only seven genes in which aTSS are found exactly at the same relative position (*D* = 0): *aglM2*, *rpl10e*, *eEF1A*, *rpoB1*, *ndhG4*, *nirH*, and *atpI.* Two of them are also equal in absolute positions since they have the same length in both species (highlighted in Table S13).

Currently, *H. volcanii* is the only other halophilic archaeon for which primary transcriptome is available (dRNA-seq data). Regular strand-specific comparative transcriptome analysis (RNA-seq) is available for other halophiles but only one, in *Natrinema* sp. J7-2 (formerly *H. salinarum* J7), focused on salinity adaptation questions [87]. Since regular RNA-seq does not have the resolution power to pinpoint aTSS as dRNA-seq, we inspected differential read distributions between *Natrinema* sp. J7-2 cultures grown on relatively low vs. high salt concentrations (15% vs. 30% NaCl) guided by the *H. salinarum* and *H. volcanii* conservation results. From the seven aforementioned genes bearing conserved aTSS, only *nirH* (which encodes a sirohaem decarboxylase) is clearly upregulated in low salt relative to high salt, while its putative conserved asRNA is upregulated in high salt relative to low salt concentrations (Figure S12).

Recently published oxidative stress data in *H. volcanii* [27] was investigated to find an intersection between differentially expressed asRNAs in *H. salinarum* NRC-1. Among the 160 conserved asRNAs, four and 15 were down- and upregulated in the *H. volcanii* oxidative stress response (Table S14). From these, two were also found differentially expressed in 37 h/17 h transition: antisense to *rmeS*,