Determining virus-host interactions and glycerol metabolism profiles in geographically diverse solar salterns with metagenomics

Obtaining metagenomic data sets and haloviral genomes

The Santa Pola, Isla Cristina, and Cahuil microbial metagenomes were downloaded from the NCBI Sequence Read Archive (SRA). Technical details about these microbial metagenomes are provided in Table 1. Further details about sampling are described in associated publications (Ghai et al., 2011; Fernandez et al., 2013; Fernández et al., 2014b; Fernández et al., 2014c; Plominsky et al., 2014). Chula Vista microbial metagenomes were obtained from the iMicrobe Collaborative (Hurwitz, 2014) and have been previously examined in a larger study of several aquatic biomes (Rodriguez-Brito et al., 2010). Chula Vista metagenomes were grouped by corresponding salt concentration (low, medium, or high) and concatenated to form single low (CV6-8), medium (CV12-14), and high (CV27-30) salt concentration metagenomes, as shown in Table 2. The ways combined microbial metagenomes were constructed and utilized in later analyses are outlined in Table 3. Santa Pola and Isla Cristina metagenomes (SS13, SS19, IC21, SS33, and SS37) were analyzed individually except for assembly and de novo CRISPR detection steps. The Cahuil (C34) metagenome was analyzed individually throughout. All analyses performed on the microbial metagenomes are outlined in Fig. S1.

Haloviral genomes as well the Lake Tyrell viral contigs were downloaded from NCBI GenBank (Table S1). Additionally, Chula Vista saltern viral metagenomes were obtained through iMicrobe (Hurwitz, 2014) and then combined into one viral metagenome (Table  3). This combined viral metagenome (Table 3) was then assembled with Newbler version 2.9 into Chula Vista viral contigs that were also included in the haloviral library (Table S1).

Profiling taxonomic and functional compositions of selected metagenomes

The selected individual metagenomes (Tables 1 and 2) were profiled for relative taxonomic abundances with MetaPhyler (Liu et al., 2011a). MetaPhyler aligns reads against a library of clade-specific marker genes. Reads that align against a marker gene above a particular threshold alignment score are then binned to the corresponding taxon (e.g., Haloquadratum, for reads aligning to Haloquadratum marker genes). Taxonomic profiles were further compared in MATLAB using principal component analysis (PCA). Clusters identified with PCA were evaluated for statistical significance using the ADONIS R function from the R package vegan (Oksanen et al., 2016) loaded in R version 3.2.2 (R Core Team, 2015). The significance of salinity in structuring the similarities between sample taxonomic profiles was determined using the envfit and ordisurf functions from the R package vegan (Oksanen et al., 2016), which fit vectors and smooth surfaces onto an ordination, respectively. Envfit was used with default parameters, while the number of knots (degrees of freedom plus one) for ordisurf was set to nine due to the nine metagenomes examined. The percent of reads mapping to markers for each metagenome were also calculated and compared to detect differences in the proportions of reads used for taxonomic profiling depending on metagenome.

The metagenomes were further profiled for functional abundances with ShotMAP (Nayfach et al., 2015) and metAnnotate (Petrenko et al., 2015). ShotMAP determination of gene family abundance was conducted with a Hidden Markov Model (HMM)-based search algorithm against an indexed library of Pfam family HMMs. Abundances were calculated from the numbers of metagenomic reads mapping to particular Pfam families normalized by the length of each Pfam target gene in the gene family. We selected abundances of nine gene families known to be related to glycerol metabolism for further analysis—dihydroxyacetone kinase (Dak1—PF02733; Dak1_2—PF13684; Dak2—PF02734), NAD⁺-dependent glycerol dehydrogenase (Fe-ADH—PF00465; Fe-ADH_2—PF13685), NAD⁺-dependent glycerol-3-phosphate dehydrogenase (both N- and C-terminal domains: NAD_Gly3p_dh_N–PF01210; NAD_Gly3p_dh_C–PF07479), and glycerol kinase (FGGY_N–PF00370; FGGY_C–PF02782), listed both with Pfam family names and ID numbers. Relative abundances were then analyzed for correlations with respect to salinity using Microsoft Excel.

metAnnotate (Petrenko et al., 2015) was used to determine the abundance of particular Pfam gene families for particular taxa. Like ShotMAP, metAnnotate uses an HMM-based search algorithm to find reads that match certain Pfam gene family HMMs. Unlike ShotMAP, however, it also determines the taxonomic affiliations of these gene family reads, and thus determines which taxa have which gene families. We used metAnnotate to taxonomically profile nine aforementioned gene families known to be related to glycerol metabolism. We then calculated the proportion of reads mapping to each taxon among all those mapping to a gene family to determine the taxonomic profiles of these gene families. These gene family profiles were then visualized as heatmaps and hierarchically clustered with the R pheatmap package (Kolde, 2015). Clusters were evaluated for statistical significance with the ADONIS R function (Oksanen et al., 2016).

De novo and reference-guided CRISPR detection

CRISPR arrays were detected using two different methods (reference-guided and de novo). All CRISPR detection analyses were treated as qualitative comparisons amongst sites because of limited sample size and inability to fully correct for differences in sequencing depth amongst metagenomes. Because CRISPR arrays’ repetitive nature makes them highly difficult to assemble, CRISPR arrays were detected purely from unassembled metagenomic reads. De novo detection was performed with Crass (Skennerton, Imelfort & Tyson, 2013) with default parameters, which identifies repeat-spacer-repeat patterns in individual reads, uses these as “seeds” to recruit other CRISPR-containing reads, and then assembles these reads into CRISPR arrays. We further clustered Crass spacers using CD-HIT (Li & Godzik, 2006) with a similarity threshold of 0.9 and aligned them against known direct repeats from CRISPRdb (Grissa, Vergnaud & Pourcel, 2007) to determine the associated host  organism.

Because Crass does not directly indicate the putative host encoding the detected CRISPRs, we developed a pipeline (MetaCRAST) to constrain spacer detection by expected host species (Moller & Liang, 2016). MetaCRAST searched each read for direct repeat (DR) sequences matching query DRs specified by the user. These query DRs were selected from known CRISPR DRs. In general, the query DR sequences might be known to be present in CRISPR arrays found in genomes expected from the taxonomic profile; for example, known genomic Haloquadratum walsbyi DRs might be used to search a metagenome for Haloquadratum walsbyi CRISPRs. Direct repeats within a certain Levenshtein edit distance (i.e., number of insertions, deletions, or substitutions necessary to convert one sequence to another) were identified in reads using the Wu-Manber multi-pattern search algorithm (Wu, Manber & Myers, 1995).

In the first step of the pipeline, reads that contain a single DR sequence matching the user-specified DR were detected with the Wu-Manber algorithm and collected for further analysis. In the second step of the pipeline, individual reads gathered in the first step were searched for two or more copies of the query direct repeat sequence within appropriate distances from each other. In the case of CRISPRs, direct repeats have been found to be between 23 and 47 bp long, while spacers, and thus the distances between DRs, have been found to be between 26 and 50 bp long (Skennerton, Imelfort & Tyson, 2013). In the third step, the sequence fragments between the direct repeats were extracted as the spacers. These spacers were then collected as a non-redundant set in a hash, and then further clustered into non-redundant sets using CD-HIT (Li & Godzik, 2006).

MetaCRAST was implemented in Perl and takes a set of query direct repeats and a metagenome (each in FASTA format) as file inputs. It is publically available for download at GitHub (https://github.com/molleraj/MetaCRAST), where all command line parameters are also documented. It has also been validated against simulated metagenomes and tested against two real metagenomes corresponding to acid mine drainage (AMD) and enhanced biological phosphorus removal (EBPR) ecosystems (Moller & Liang, 2016).

We performed reference-guided CRISPR detection with MetaCRAST using the following parameters. We used a group of 29 DR sequences previously detected in Haloferacales and Halobacteriales genomes (Table S4) in the genomic CRISPR database CRISPRdb (Grissa, Vergnaud & Pourcel, 2007) as query DRs. We thus searched for reads containing sequences matching any of these 29 direct repeat sequences. The maximum allowed distance between detected direct repeats (i.e., maximum spacer length) was set to 60 bp. Acceptable edit distance—the allowed differences (insertions, deletions, or substitutions) between query and detected read DR sequences—were adjusted from 0 to 3 to see how variation in edit distance affected reference-guided CRISPR detection. All spacers detected for each condition (for each metagenome and acceptable error number) were then clustered using CD-HIT with a clustering similarity threshold of 0.9, which effectively clusters spacers together into a non-redundant set based on a pairwise similarity of 90% or greater.

Mapping virus-host interactions with CRISPR spacers

CRISPR spacers extracted from the metagenomes with Crass and reference-guided methods and clustered with CD-HIT were then aligned against the library of haloviral genomes using NCBI BLAST (Altschul et al., 1990) with an e-value cutoff of 1e-03. All spacers detected with Crass for the combined Santa Pola dataset (see Table 3) were aligned against this library to create its corresponding interaction network. For several individual metagenomes (e.g., SS13, SS19), no Crass spacers mapped to the haloviral genomes, making the more sensitive reference-guided search (MetaCRAST) necessary. Thus, for individual metagenomes, spacers detected with MetaCRAST using a maximum Levenshtein edit distance of 3 were aligned for network creation. BLAST results from each method were then processed to form lists of virus-host pairs. Networks of virus-host interactions were drawn with Cytoscape (Shannon et al., 2003).

Comparative analysis of virus and host sequence characteristics

Viruses and hosts matched through CRISPR mapping were further compared in nucleotide k-mer frequencies to determine whether virus and host frequencies were similar. Dinucleotide and trinucleotide usages were determined using the function nmercount from the MATLAB Bioinformatics Toolbox (The MathWorks, Inc., 2015). Dinucleotide and trinucleotide frequencies were compared for a non-redundant list of viruses associated to each of the six detected hosts listed below. All viruses were associated to possible hosts using our in-house reference-based CRISPR detection method. Haloquadratum walsbyi C23, Haloferax volcanii DS2, Halogeometricum borinquense DSM 11551, Halorhabdus utahensis DSM 12940, Haloarcula marismortui ATCC 43049 (chromosome 1), and Haloarcula hispanica N601 (chromosome 1) host genomes were analyzed along with the matched haloviral genomes. To measure the sequence characteristic differences among virus-host pairs, each virus dinucleotide and trinucleotide frequency was subtracted from each corresponding host genome frequency and squared. All these squared distances were summed for possible dinucleotides and for possible trinucleotides to generate an overall Euclidean measure of distance between virus and host.

Assembling metagenomes and detecting halobacterial cas genes in metagenomic contigs

Combined Santa Pola, combined Chula Vista, and Cahuil metagenomes (Table 3 outlines combined metagenomes) were assembled into contigs to detect cas operons. Combined Santa Pola, combined Chula Vista, and Cahuil metagenomes were assembled both with Newbler and Velvet, creating a total of six sets of contigs (Zerbino & Birney, 2008). Both Newbler and Velvet were used to compare the effects of Overlap-Layout-Consensus (Newbler) and de Bruijn graph (Velvet) assembly methods. Newbler assembly was performed with default parameters, while Velvet assembly was performed with a k-mer size of 31.

Newbler and Velvet contigs assembled from each combined metagenome were first aligned against cas genes annotated in Haloferax volcanii DS2, Halorubrum lacusprofundi ATCC 49239, and Haloquadratum walsbyi C23 genomes using NCBI blastn (Altschul et al., 1990) with an e-value cutoff of 1e-06 and NCBI tblastx with an e-value cutoff of 1e-30. Color-coded cas detection results for genes in the Haloquadratum, Halorubrum, and Haloferax cas operons were drawn with the synthetic biology tool Pigeon (Bhatia & Densmore, 2013). A library of cas genes from the three major halobacterial orders (the Haloferacales, Halobacteriales, and Natrialbales) were downloaded from NCBI Genbank (see Tables S11–S13) and aligned against each set of assembled contigs (assembled with each method and from each combined metagenome) to determine the taxonomic affiliation of matching cas genes. BLAST alignments were performed with the same parameters described above. Taxonomy information was extracted from matching cas gene Genbank records with the NCBI Efetch utility to determine the taxonomic affiliations of the assembled cas contigs.

Taxonomic profiling identifies predominant members of saltern microbial communities

We profiled the taxonomic abundances of metagenomes sequenced from Cahuil (C34), Chula Vista (CV6-8, CV12-14 and CV27-30), Isla Cristina (IC21), and Santa Pola Saltern (SS13, SS19, SS33, SS37) sites (where the ID numbers indicate salinity percentages) to better understand the population structures of saltern microbial communities and to compare host abundances amongst sites. We found the low salinity Chula Vista (CV6-8) and SS13 salterns had different microbial communities from the other salterns, containing large proportions of bacterial Actinobacteria, Proteobacteria and Rhodobacteraceae (7.64, 14.28, and 20.69% in SS13 and 6.6, 10.05, and 11.13% in CV6-8, respectively—see Fig. 1A). In addition, we found Haloquadratum was the major taxon identified in the SS19, SS33, and SS37 salterns (88.32, 65.55, and 78.06% abundance, respectively), while Halorubrum and Halobacteriaceae were major components of the Cahuil and IC21 salterns (31.79 and 64.55% in Cahuil; 32.5 and 51.63% in IC21, respectively–see Fig. 1A). Bacterial Salinibacter was identified in the SS19, SS33, and SS37 salterns (2.18, 7.16, and 6.82% relative abundances, respectively). While Haloquadratum was the major identified taxon in the higher salinity Santa Pola salterns, Halobacteriaceae composed the most abundant identified taxon identified in the medium and high salinity Chula Vista salterns (51.46 and 77.02%, respectively). We realize that the Halobacteriaceae component of the saltern communities may also include species in the new order Haloferacales (Gupta, Naushad & Baker, 2015), but this taxon was defined after the MetaPhyler marker library was constructed (Liu et al., 2011a). These taxonomic profiles and clustering results further guided our analyses of host CRISPR arrays, as we used taxonomy to guide our searches for CRISPR spacers.

We further compared the taxonomic profiles through principal component analysis, which found that the first two principal components accounted for 88% of variation amongst the examined taxonomic profiles. PCA grouped the SS19, SS33, and SS37 profiles; low salinity Chula Vista (CV6-8) and SS13 profiles; and medium salinity Chula Vista (CV12-14), high salinity Chula Vista (CV27-30), IC21, and C34 profiles (Fig. 1B). These groups were found to be statistically significant in their between-group differences (ADONIS; F=19.817; p = 0.002). We furthermore found percent salinity to be a significant feature in structuring the similarities amongst metagenomes when fitting this variable to the PCA ordination either as a vector or a response surface (envfit: p = 0.038; ordisurf, knots=9: p = 0.0291).

We also determined the percentage of reads mapping to phylogenetic markers to detect read length-dependent differences in taxonomic profiling among metagenomes (Fig. S2). The percentage mapping was similar among Santa Pola, Isla Cristina, and Cahuil (C34) metagenomes (an average of 0.51% of reads with a standard deviation of 0.0483% among these metagenomes), but lower among the Chula Vista metagenomes (an average of 0.21% of reads with a standard deviation of 0.041% among these metagenomes) due to lower coverage (see Tables 1 and 2).

We have thus built upon prior metagenomics studies to show that Haloquadratum is the majority taxon in many salterns (Oh et al., 2009), particularly those that are environmentally stable with high salinity (Fig. 1). The Santa Pola and Isla Cristina metagenomes have previously been taxonomically profiled through mapping of reads to 16S rRNAs obtained from the Ribosomal Database Project (Fernandez et al., 2013; Fernández et al., 2014b; Ventosa et al., 2014). These 16S rRNA profiling studies revealed that the lowest salinity saltern (SS13) contained a community of bacterial Actinobacteria and Proteobacteria and archaeal Euryarchaeota, while members of the Euryarchaeota increased in abundance with increasing salinity (Fernández et al., 2014c; Fernández et al., 2014a; Ventosa et al., 2015). Compared with the Santa Pola salterns, Cahuil and Isla Cristina salterns may be more dynamic, with more rapidly changing daily and yearly temperatures and higher rainfall (Fernández et al., 2014c). This may lead to lower concentrations of potassium, magnesium, and sulfate in those salterns, thus leading to lower abundances of Haloquadratum, and higher abundances of other haloarchaea (Podell et al., 2014; Fernández et al., 2014c).

Functional profiling reveals abundances of glycerol metabolism genes and their taxonomic affiliations

Glycerol is the main carbon source for saltern microbes (Oren, 1993). We used two functional profiling tools to study the glycerol metabolism genes critical for carbon cycling in these ecosystems. We used the tool ShotMAP to determine abundances of glycerol metabolism gene families (Nayfach et al., 2015) and metAnnotate to study the taxonomic distribution of reads mapping to these gene families (Petrenko et al., 2015). As with MetaPhyler, Cahuil, Isla Cristina, and Santa Pola metagenomic reads mapped to the ShotMAP reference library in higher proportions than the Chula Vista reads (Fig. S3). Using these tools, we examined the abundances and taxonomic affiliations of the glycerol kinase (FGGY_N and FGGY_C), NAD⁺-dependent glycerol-3-phosphate dehydrogenase (NAD_Gly3P_dh_N and NAD_Gly3P_dh_C), glycerol dehydrogenase (iron-containing alcohol dehydrogenase –Fe-ADH and Fe-ADH_2), and dihydroxyacetone kinase (Dak1, Dak1_2, and Dak2) Pfam gene families.

We first determined gene family abundances with ShotMAP and then examined relationships between these abundances and salinity levels. Among the gene families we examined, we noticed both negative and positive correlations between relative gene family abundance and salinity (Fig. 2). Dak1, Dak2, FGGY_N, and Fe-ADH_2 families showed positive correlations with salinity, while Dak1_2, NAD_Gly3P_dh_N, NAD_Gly3P_dh_C, FGGY_C, and Fe-ADH families showed negative correlations with salinity. Only Dak1_2, NAD_Gly3P_dh_N, and NAD_Gly3P_dh_C gene families showed significant correlations with salinity (r =  −0.617, −0.795, and −0.787; p = 0.038, 0.005, and 0.006, respectively). We also noted family abundances calculated for the Chula Vista metagenomes to be larger than those found for any other metagenome. More short reads mapping to these gene families may account for these higher family abundances than those found for the Santa Pola, Isla Cristina, and Cahuil metagenomes.

We next determined the taxonomic affiliations of reads mapping to these gene families using metAnnotate. We generated heat maps of these affiliations hierarchically clustered both with regard to metagenome (Fig. 3A) and gene family (Fig. 3B). As with taxonomic profiles, the taxonomic affiliations of glycerol metabolism genes in the Haloquadratum-rich SS19, SS33, and SS37 metagenomes group together, as do those of the Chula Vista metagenomes, and the SS13, IC21, and Cahuil (C34) metagenomes (Fig. 3A). Across gene family, taxonomic affiliations of the dihydroxyacetone kinase families were closest to each other, as well as those of the alcohol dehydrogenase (Fe-ADH) and glycerol kinase (FGGY) gene families (Fig. 3B). The two NAD⁺-dependent glycerol-3-phosphate dehydrogenase families grouped together, but many reads mapping to these families were found to be unclassified (Fig. 3B). These described gene family and metagenome clusters were found to be statistically significant (ADONIS; F = 7.5468 and 2.8768; p = 0.014 and 0.013, for gene family and metagenome clusters, respectively).

We also analyzed these taxonomic affiliations for each gene family across each individual metagenome (Fig. 3C). Haloquadratum associated with dihydroxyacetone kinase gene families in the SS19, SS33, and SS37 metagenomes while Halorubrum associated with these families in the IC21 and Cahuil metagenomes (Fig. 3C). The two NAD⁺-dependent glycerol-3-phosphate dehydrogenase families, on the other hand, associated mainly to Salinibacter in the Cahuil (C34), high salinity Chula Vista, SS19, SS33, and SS37 metagenomes and to Spiribacter in the SS13 and IC21 metagenomes. Unlike the previous families, alcohol dehydrogenase and glycerol kinase had less clear associations across the metagenomes, with association to Haloquadratum and Halorubrum but also many unclassified reads. In addition, glycerol kinase families associated to Spiribacter (SS13 and IC21) and the alcohol dehydrogenase family to Halobellus (C34, CV27-30, SS19, SS33, and SS37) across several metagenomes.

Our functional studies thus found both positive and negative correlations between glycerol metabolism gene family abundances and salinity (Fig. 2). Prior studies of the Santa Pola and Isla Cristina salterns found that glycerol metabolism genes, such as dihydroxyacetone kinase, glycerol kinase, and NAD⁺-dependent glycerol-3-phosphate dehydrogenase, increased in abundance with salt concentrations (Fernández et al., 2014c; Fernández et al., 2014a). Taxonomic affiliations, however, differed amongst salterns and different gene families (Fig. 3). Dihydroxyacetone kinase most strongly associated with Haloquadratum, except in the IC21 and Cahuil (C34) salterns, where it associated strongly with Halorubrum instead. NAD⁺-dependent glycerol-3-phosphate dehydrogenase, on the other hand, associated with Salinibacter (Fig. 3C). Based on prior understanding of saltern glycerol cycling (Oren, 1993; Oren, 1994; Sher et al., 2004), these results together suggest archaea such as Haloquadratum and Halorubrum consume the glycerol metabolism intermediate dihydroxyacetone in these ecosystems, while bacteria such as Salinibacter may first phosphorylate and then oxidize glycerol to metabolize it (Sher et al., 2004).

Metagenomic CRISPRs detected with two independent methods show variation in CRISPR arrays depending on community composition

We detected CRISPR arrays in metagenomes using both de novo and reference-guided methods. We used Crass for de novo detection (Skennerton, Imelfort & Tyson, 2013) and a custom pipeline (MetaCRAST) for reference-guided detection (see Materials and Methods). We used a combination of Halobacteriales and Haloferacales CRISPR direct repeats downloaded from CRISPRdb as our query for MetaCRAST analyses across all metagenomes.

Overall, more spacers were detected in IC21 and Cahuil (C34) metagenomes compared with the high salinity Chula Vista, SS13, SS19, SS33, and SS37 metagenomes (Fig. 4) both with MetaCRAST (3 errors allowed) and Crass. More spacers were detected with Crass compared to reference-guided methods in the combined Santa Pola, SS13, and SS33 metagenomes (377 vs. 277, 29 vs. 9, and 146 vs. 80 spacers, respectively). For the reference-guided detection analyses, the number of spacers identified increased as the allowable error (insertions, deletions, or substitutions different from the query DR sequence) in alignment was increased from 0 to 3 (Fig. 4). No spacers were detected with either method in the low or medium salinity Chula Vista metagenomes. When we normalized by total number of metagenomic reads, we found that the Cahuil (C34) metagenome had the largest number of spacers detected per million reads (Fig. S4). Nonetheless, it is important to note that CRISPR detection analyses could only be treated as qualitative rather than quantitative because of limited sample size and inability to fully correct for differences in sequencing depth amongst metagenomes.

We then compared the detected direct repeats and spacers among metagenomes to discover common CRISPR array types or spacers (Figs. 5 and 6). Remarkably, common spacers were identified amongst the SS19 and SS37 metagenomes (13 with Crass and 18 with MetaCRAST among these two, respectively) using both detection methods (Figs. 5A and 5B). Amongst the metagenomes with a higher number of detected spacers (IC21 and Cahuil/C34), common spacers were detected with both methods (5 in IC21 and 24 in Cahuil—see Figs. 5C and 5D). Amongst the direct repeats detected with Crass in the combined Chula Vista, combined Santa Pola, and Cahuil metagenomes, there was one direct repeat (DR) sequence common to both Cahuil and combined Santa Pola (Fig. 6). We also searched the combined Santa Pola direct repeats against a full database of direct repeats from CRISPRdb to determine possible taxonomic affiliations. We found that 14 of the 67 detected combined Santa Pola DRs matched DRs in the CRISPRdb database (Table S3). Of these 14 matching DRs, 10 matched haloarchaeal DRs, while 4 matched bacterial DRs.

Our CRISPR detection results showed a smaller number of CRISPR spacers found in salterns where Haloquadratum was identified as the majority taxon compared with salterns where Halobacteriaceae was found to be the majority taxon. Furthermore, comparison of detected metagenomic spacers to those found in the genomes of two different H. walsbyi strains (C23 and DSM 16790) suggests variation in CRISPR array length in the salterns. Strain C23 contains a full set of cas genes and several CRISPR arrays, while strain DSM 16790 lost all its cas genes and all but five spacers (Dyall-Smith et al., 2011). While we could not identify any C23 spacers amongst our metagenomes, we could find at least three DSM 16790 spacers in several metagenomes. A successful mobile element integration event in spite of CRISPR defense is thought to have promoted loss of spacers in Haloquadratum genomes (Dyall-Smith et al., 2011). Taken together, this evidence suggests genomic CRISPR arrays in Haloquadratum species may vary in length.

Metagenomic CRISPRs detected with two independent methods map haloviruses to saltern hosts

We then aligned the spacers detected with both methods against a library of haloviral genomes (Table S1) to determine virus-host interactions in these salterns. Unlike CRISPR spacers detected with Crass in several individual Santa Pola metagenomes, spacers detected in a combination of all Santa Pola and Isla Cristina metagenomes did align to haloviral genomes (Fig. 7). In this Crass analysis of the combined Santa Pola dataset, we detected virus interactions most likely with Haloquadratum (G467), Haloarcula (G42 and G82), and Saccharomonospora (G57) hosts, based on our previous comparison of Crass DRs (listed by group number, such as G467) against CRISPRdb DRs (Fig. 7). For many of the interactions (e.g., G467, G42), we could determine a possible host by such comparison. For several notable groups of interactions, however, we could not determine possible hosts (e.g., G85, G1290).

While Crass only provided a picture of virus-host interactions in the combined Santa Pola metagenome, MetaCRAST detected spacers that mapped to haloviral genomes in each of the individual Santa Pola, Isla Cristina, and Cahuil metagenomes (Fig. 8). Only one interaction was detected in the SS13 dataset (Fig. 8A), but there were strikingly similar groups of virus-host interactions detected in the SS19, SS33, and SS37 metagenomes (Figs. 8B– 8D). While many groups of interactions were similar to those expected from the Haloquadratum walsbyi strain C23 CRISPR array (Fig. S5), no C23 spacers matched any detected metagenomic spacers, while at least three Haloquadratum walsbyi strain DSM 16790 spacers did match spacers detected in SS19, IC21, SS33, and SS37 metagenomes. Aside from Haloquadratum, in the SS33 metagenome, we also detected Haloarcula-virus interactions with eHP-4, eHP-8, and eHP-13 (Fig. 8C). Quite surprisingly, we noticed remarkably different virus-Halobacteriaceae interactions in the IC21 and Cahuil (C34) metagenomes (Figs. 8E and 8F), perhaps due to differences in taxonomic composition at the genus level.

Our CRISPR mapping results not only showed consistent viral interactions with Haloquadratum across salterns but also numerous interactions with less abundant taxa such as Haloferax and Haloarcula. As expected from our comparison of taxonomic and functional profiles, there were consistent sets of interactions amongst SS19, SS33, and SS37 salterns, where Haloquadratum species form the majority of the microbial population (Fig. 8). These interactions suggest similar viral infection patterns across these environments, consistent with “global pond” theories of microbial and viral exchange (Dyall-Smith et al., 2011; Atanasova et al., 2012). Unclear biological requirements for CRISPR spacer incorporation make it difficult to make any conclusion about the abundances of the viruses. Though some models suggest viruses must have reached a critical number within their hosts to induce CRISPR spacer incorporation into host genomes (Garrett, Vestergaard & Shah, 2011), characterization of the Haloferax volcanii CRISPR systems suggests plasmid copy number has little effect on incorporation (Maier et al., 2013).

There were also interactions detectable with Crass that could not map to any known hosts (Fig. 7). Some of these hosts may be poorly characterized members of the Nanohaloarchaea, as recent studies have suggested for eHP-12 and eHP-16 (Martínez-García et al., 2014). Surprisingly, we did not detect any spacers or virus-host interactions with Halorubrum with either CRISPR detection method. Either limitations of our detection tools or a lack of Halorubrum CRISPRs could account for this observation.

Nucleotide usage comparison suggests biological plausibility of virus-host associations

To determine whether our virus-host mappings are biologically plausible, we compared the dinucleotide and trinucleotide frequencies of CRISPR-determined virus-host pairings. Over evolutionary time, viruses’ nucleotide usage steadily approaches that of the host, reflecting adaptation to host codon usage, tRNA availability, intracellular nucleotide pools, and host restriction modification systems (Edwards et al., 2015). We compared these nucleotide usage frequencies both among the predicted hosts and analyzed metagenomes (Fig.  9). Among predicted hosts, Haloquadratum walsbyi, Haloferax volcanii, and Haloarcula marismortui showed the smallest overall differences in dimer and trimer usage, suggesting strong virus-host association, while Haloarcula hispanica and Halorhabdus utahensis showed the highest overall differences (Figs. 9A and 9B). Only one virus associated to Halogeometricum borinquense, and it had a greater difference in trimer frequency compared with dimer frequency relative to the other predicted hosts. Among the metagenomes, the smallest differences were found in the SS33 metagenome; the next smallest differences were found in the Cahuil (C34), SS19, and SS37 metagenomes, while the largest difference was found in the single virus-host pair detected in the SS13 metagenome (Figs. 9C and 9D).

Assembled metagenomic contigs contain halobacterial cas genes

To further confirm our CRISPR detection and mapping results, we searched through metagenomic contigs for the cas genes required for CRISPR activity. Because all microbial metagenomes studied were sequenced with 454 technology, and Newbler is designed to assemble 454 reads, we expected longer contigs to be assembled with Newbler and found this in our assembly results (Table S2). We both searched for contigs that matched several different haloarchaeal cas operons (Fig. 10A) and profiled the taxonomic affiliations of haloarchaeal cas genes that best matched contigs (Fig. 10B). We were able to assemble almost all genes in cas operons in three different metagenomes (Cahuil, combined Chula Vista, and combined Santa Pola) from several different Type I-B and Type I-D cas operons (Fig. 10A).

We also searched a library of cas genes from three different haloarchaeal orders (Halobacteriales, Haloferacales, and Natrialbales) and generated profiles of taxonomic affiliations for the genes that best matched contigs (Fig. 10B). Haloarcula composed a large proportion of Halobacteriales contig affiliations, especially amongst the Cahuil (C34) contigs (71.4% among Newbler contigs; 41.4% among Velvet contigs). On the other hand, Haloquadratum and Halorubrum composed large proportions of the Haloferacales affiliations, especially for the combined Santa Pola dataset (14.5% and 25.8% for Newbler contigs; 23.5% and 26.6% for Velvet contigs). We also calculated the total number of hits in the search against the cas gene libraries (Fig. S6) and profiled the taxonomic affiliations of all cas gene alignments to the contigs (Fig. S7). There were smaller Haloarcula, Halonotius, Halorhabdus, and Haloquadratum proportions in the total hit profiles (Fig. S7) compared with the best hit profiles (Fig. 10B).