In the current survey, we examined structural and functional responses of the bacterioplankton community toward a phytoplankton bloom. Samples for community analysis were taken randomly at different locations and different depths within a P. globosa bloom in the German Bight (**Figure 1**, **Table 1**). Six samples were taken in presence of the phytoplankton bloom (samples 5, 6, 9, 10, 13, and 15) and three in bloom absence (samples 3a, 3b, and 4). One sample was taken near the Weser river outfall (sample 1). Salinity ranged from 30.7 to 32.7 psu. Fluorescence was approximately 0.45 and 2.2 mg/m3 outside and inside the algal bloom, respectively. Temperatures ranged from 8.2 to 11.8◦C. All environmental parameters are listed in **Table 2**. Based on our previous analysis, most measured parameters were significantly linked to algal bloom presence (see Wemheuer et al., 2014). Only the suspended particulate matter content (SPM) and the nitrite concentration exhibited no direct correlation to bloom presence.

### Bloom Presence Affects Bacterial Community Structures

Total and active bacterioplankton community structures were assessed by pyrosequencing-based analysis of the V3–V6 region of the 16S rRNA amplified from environmental DNA and RNA, respectively. A total of 50,125 and 31,982 high-quality bacterial 16S rRNA sequences were obtained across all 10 samples at DNA and RNA level, respectively (**Table 6**). Calculated rarefaction

TABLE 7 | Alpha diversity indices at 97 and 80% genetic similarity derived from the 16S rRNA analysis.

OTU, operational taxonomic unit.

curves (data not shown) as well as diversity indices revealed that the majority of the bacterial community was recovered by the surveying effort (**Table 7**).

Classification of the obtained 16S rRNA sequences revealed that Proteobacteria and Bacteroidetes were the most abundant bacterial phyla across all samples (approximately 78% and 20%, respectively). At higher taxonomic resolution, the majority of the obtained sequences was affiliated to 17 bacterial groups, clades, and genera (**Figure 2**). These groups represented different lineages within the Alpha-, Beta-, and Gammaproteobacteria and the Bacteroidetes. These results are in accordance with our previous study (Wemheuer et al., 2014) and recent investigations of bacterial communities in the North Sea (Alderkamp et al., 2006; Sapp et al., 2007; Teeling et al., 2012). However, in our previous study, the number of Bacteroidetes was rather low which can be attributed to the differences in primer pairs and variable regions of the 16S rRNA gene used in our previous study.

Alphaproteobacteria accounted for 50% of all sequences with a higher abundance at DNA and RNA level (59 and 41%, respectively). The opposite was recorded for the Gammaproteobacteria, which accounted for 16% (DNA) and 31% (RNA), respectively. The increased abundance of the Gammaproteobacteria was mainly attributed to the higher abundances of Pseudospirrilii and the SAR92 clade. Changes in the abundances of the different alphapoteobacterial taxa were mainly attributed to the overall low abundance of the SAR11 clade at RNA level. A low activity of SAR11 is supported by other studies (West et al., 2008; Lamy et al., 2010; Klindworth et al., 2014). For example, Lamy et al. found an overall low abundance and activity of the SAR11 clade in a P. globosa bloom in the eastern English Channel. In another study, Alonso and Pernthaler (2006b) showed that SAR11 is highly abundant but not very active in costal North Sea waters. In addition, West et al. (2008) demonstrated that SAR11 was more abundant at DNA level than at RNA level in the Southern Ocean.

Several groups responded significantly toward algal bloom presence at DNA and/or RNA level, e.g., the abundance of the SAR92 clade was three times higher at RNA level and in bloom presence. This is in accordance with previous studies (Pinhassi et al., 2005; West et al., 2008; Klindworth et al., 2014; Wemheuer et al., 2014). For example, a phytoplankton bloom induced by inorganic nutrient enrichment influenced SAR92 in a mesocosm experiment (Pinhassi et al., 2005). Klindworth et al. (2014) found that members of the Rhodobacteraceae and SAR92 clade exhibited high metabolic activity levels during a bloom succession, which indicates their important role during bloom events. In addition, the 16S cDNA estimates for SAR11 were notably lower in the earlier bloom sample. The authors suggest that members of this clade could not profit from the increasing availability of nutrients in the decaying bloom and thus were outcompeted by other clades. This is in line with our study in

which we found significant lower abundances of SAR11 in the bloom presence and at RNA level.

Members of the two genera Marinoscillum and Polaribacter were significantly more abundant in bloom samples both at DNA and RNA level. Bacteroidetes are widespread in marine systems and play an important role in organic matter degradation (Gómez-Pereira et al., 2010). The higher abundance of this phylum during the phytoplankton bloom was verified by recent findings (Alderkamp et al., 2006; Lamy et al., 2010; Tada et al., 2012; Teeling et al., 2012). The strongest increase in activity during the senescent stage of a P. globosa bloom in the North Sea was observed for Bacteroidetes (Alderkamp et al., 2006) A mesocosm experiment targeting bacterial succession patterns during a diatom bloom revealed that Bacteroidetes had a relatively high growth potential as the bloom peaked (Tada et al., 2012). The authors suggested that the early development contributed to the initial stage of bloom decomposition. Therefore, this phylum seems to benefit from the conditions provided by the algal bloom and might play an important role in the degradation of phytoplankton-derived organic matter. Klindworth et al. (2014) mapped metatranscriptomic data on assembled and taxonomically classified metagenomic data and found that Formosa and Polaribacter acted as major algal polymer degraders. A similar conclusion was drawn in a study of a P. globosa bloom in the eastern English Channel. Here, members of Bacteroidetes group dominated the activities and the abundances during the growth phase of the algae (Lamy et al., 2010).

### Bloom Presence Affects Bacterioplankton Gene Expression

After removal of ribosomal RNA, nearly 45 million Illumina reads remained and were used for environmental gene expression analysis. Generated mRNA datasets were initially mapped on the 15 reference genomes belonging to abundant marine genera and lineages. However, only 10% of the sequences mapped to these reference genomes. Most of these sequences were affiliated to the genome of P. temperata RCA23 (see Supplementary Table S3). This strain was isolated in the German Wadden Sea (Giebel et al., 2011, 2013), and its genome was recently described (Voget et al., 2014).

Mapping mRNA datasets on genomes is a common approach when analyzing metatranscriptomic data (e.g., Gifford et al., 2013). The advantage of this approach is that community functions can be linked to a certain organism. However, most reads are not included in the analysis because reference genomes for many marine lineages are still missing. Another problem is the data normalization when mapping metatranscriptomic data on genomes. In a transcriptomic approach, all sequences derive from a single organism and data can be normalized by the number of reads mapped. However, in a metatranscriptomic approach, the amount of sequences affiliated to an organism is not only linked to its gene expression but also to the gene expression of all other community members. Thus, a decrease in abundance can be caused either by a lower gene expression of the organism or by an increased expression of other community members.

Here, we mapped the mRNA data on reference genomes representing marine lineages, which were abundant in the 16S rRNA gene analysis. These genomes included data from P. temperata RCA23 and HTCC2255, belonging to two Roseobacter lineages, the RCA and NAC11-7 clusters. Other reference genomes used were the genomes of HTCC2207 und MOLA455, both members of the SAR92 clade (Stingl et al., 2007; Courties et al., 2014). Mapping the datasets on assembled metagenomic and metatranscriptomic data resulted in an overall alignment rate of more than 86.22% (see **Table 5**). This overall high coverage is higher than in other studies. For example, Kopf et al. (2015) mapped up to 80% of their metatranscriptomic data on the corresponding metagenome. We were able to map our data on approximately 600,000 genes with almost 800,000 different functions.

Most of the genes affiliated to the two members of the SAR92 clade were upregulated which corresponds to their increasing abundance at 16S rRNA transcipt level. For example, one leucyltRNA synthetase and three isoleucyl-tRNA synthetases affiliated to MOLA455 were significantly upregulated in the bloom (**Figure 3A**). Morover, two leucine-tRNAs were significantly upregulated in the bloom in P. temperata RCA23 (**Figure 3B**). In addition, an isoleucyl-tRNA synthetase affiliated to P. temperata RCA23 (**Figure 3A**) and a isoleucyl-tRNA synthetase and a leucyl-tRNA synthetase of P. temperata RCA23 (**Figure 3B**) were marginal significantly upregulated in the bloom (P < 0.1). This is in line with a study by West et al. (2008). The authors found that the Roseobacter groups NAC11-7 and RCA as well as the SAR92 clade were the most important contributors to leucine incorporation during the peak of a naturally iron-fertilized phytoplankton bloom in the Southern Ocean. This result is confirmed by a study about a P. globosa bloom in the English Channel (Lamy et al., 2010). Here, Bacteroidetes and Gammaproteobacteria were the most abundant and active groups during the growth period of the algae. Gammaproteobacteria and Alphaproteobacteria dominated by the Roseobacter clade accounted for the major part of leucine incorporation after the disappearance of the bloom. In addition, the contributions of different bacterial groups to bulk abundance and leucine incorporation were partly correlated with cellspecific exoproteolytic and exoglucosidic activities and with particulate organic carbon. This indicates some specificity of these bacterial groups with respect to their ecological role in the environment. Interestingly, we identified two betaglucosidases affiliated to HTCC2207 being expressed only in bloom samples (**Figure 3A**). In the study of Teeling et al. (2012), metagenomic and metaproteomic data indicated the presence of distinct sets of carbohydrate-active enzymes (CAZymes) and transporters, which suggested a positive selection for bacteria with the capacity to decompose phytoplankton biomass. Four HTCC2207 indicator genes have been described to contain cadherin domains involved in complex carbohydrate degradation via cell aggregation and direct binding to cellulose, xylan, and related compounds (Gifford et al., 2013). This might explain the increase of the SAR92 clade as observed in the present study. In addition, the role of Gammprotaobacteria during polysaccharide degradation has been recently addressed in a study by Wietz

et al. (2015). The authors showed that alginate and agarose degradation covaried with the abundance of different lineages within the Gammaproteobacteria.

Interestingly, a major fraction of the expressed genes affiliated to SAR92 clade and other bacterial lineages was linked to several genes such as RNA polymerases, heat shock proteins, chaperons, sigma factors, and ribosomal proteins. This is consistent with a study from Klindworth et al. (2014). In this study, the most abundant mRNA transcripts with known functions coded for housekeeping genes including DNA gyrase, elongation factors, and sigma factors. The authors indicated that this reflects differences in nutritional ecological strategies of the dominant bacterial classes.

However, the high number of heat shock and other stressrelated genes overexpressed in bloom samples in members of the SAR92 clade might not necessarily reflect its ecological role as a polymer degrader but might also be caused by a higher stress tolerance toward the rapidly changing conditions during the phytoplankton bloom (Wemheuer et al., 2014). PON, Chl a, and phaeopigments of stations in the bloom area were significantly higher than that outside the bloom area. An increasing pH from 7.9 to 8.7 was observed as a result of CO2 net fixation into the alga during a phytoplankton bloom (Brussaard et al., 1996). Members of SAR92 clade might benefit from bloom conditions due to their high stress tolerance level rather than filling one of the specialized ecological niches formed during a phytoplankton bloom.

### Adaption to Environmental Changes by Higher Genome Plasticity

The overall expression level of P. temperata RCA23 was higher inside the bloom (**Figure 4**). Numerous genes encoding for transposases in its genome were highly overexpressed. Transposases are the most abundant and most ubiquitous genes in nature (Aziz et al., 2010). In addition, investigation of the metatranscriptomic bins revealed the presence of several transposases that were affiliated to P. temperata and overexpressed in the bloom (**Figure 4**). It has been shown that some bacteria expressed transposases under changing environmental conditions to rearrange genome architecture. For example, up to 81 genes encoded for transposases were upregulated in Microcystis aeruginosa relative to the control when grown on urea (Steffen et al., 2014). Genome rearrangements and the resulting genome mosaics have been also found in other members of the Roseobacter clade. The genomes of Octadecabacter arcticus and O. antarcticus are highly different despite their similarity on 16S rRNA gene sequence level and the presence of some unique gene features (Vollmers et al., 2013). This is attributed to genomic rearrangements caused by an unusually high number of transposases in the genomes of both Octadecabacter strains. We assume that the recorded overexpression could result in a higher genome plasticity/heterogeneity of this population and thus might be a possible adaptation strategy of P. temperata to environmental changes. Moreover, as found in other members of the Roseobacter clade, it might be one of the key features of this group explaining its high abundance in marine ecosystems and its ability to adapt to various marine niches. However, comparative genome studies are missing because only one genome of the genus Planktomarina is currently available. Consequently, this issue cannot been fully answered yet.

### Conclusions

Active bacterial communities in the North Sea are dominated by only a few marine groups such as the Roseobacter RCA cluster. Some of these lineages responded significantly toward the P. globosa bloom investigated in this study. For example, the SAR92 clade was three times more abundant at active bacterial community level and in bloom presence. The metatranscriptomic approach revealed that these groups are not dominated by well-studied isolates or type species as only 10% of all metatranscriptomic sequences mapped on reference genomes. Therefore, in situ experiments employing available isolates do not necessarily reflect environmental conditions and, thus, only provide limited information on the ecological role of the studied isolates. However, mapping these reads on assembled metagenomic and metatranscriptomic sequences led to an overall mapping rate of more than 85% demonstrating the power of this combined approach. The functional analysis performed in this study provides insights into gene expression patterns of the abundant community members. The high abundance of the SAR92 clade, which is supposed to be involved in polymer-degradation during and after the bloom, is attributed to a higher stress tolerance indicated by the high number of

heat shock expressed in the bloom. Although the number of field studies targeting the active bacterial community either by metatranscriptomic or metaproteomic approaches has been increased over the past years, the complex dynamics of marine environments are still largely unexplored. This study provides a deep insight into structural and functional responses of the bacterioplankton community toward a phytoplankton bloom. Therefore, it paved the way for a better understanding of the complex dynamics of marine bacteria and their interactions with the surrounding environment.

### Author Contributions

RD and BW conceived and designed the experiments; BW, FW, JH, SV, and FM performed the experiments and analyzed the data; BW, FW, and RD wrote the paper; all authors reviewed, edited, and approved the manuscript.

### References

### Acknowledgments

We thank the crew of the research vessel Heincke for their valuable support during the sampling campaign. We are grateful to Peter Meinicke and Heiko Liesegang for the help during data analysis. This work was funded by the Deutsche Forschungsgemeinschaft (DFG) as part of the collaborative research center TRR51 and the Alfred Wegener Institute under grant number AWI-HE327\_00. Additionally, we acknowledge support by DFG and the Open Access Publication Funds of the Göttingen University.

### Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fmicb. 2015.00805

a representative of a ubiquitous proteorhodopsin-producing group in the ocean. Genome Announc. 2, e01203–e01213. doi: 10.1128/genomea.01203-13

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2015 Wemheuer, Wemheuer, Hollensteiner, Meyer, Voget and Daniel. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Metagenome and Metatranscriptome Revealed a Highly Active and Intensive Sulfur Cycle in an Oil-Immersed Hydrothermal Chimney in Guaymas Basin

### *Ying He1,2, Xiaoyuan Feng1, Jing Fang1, Yu Zhang1,2,3 and Xiang Xiao1,2,3\**

*1 State Key Laboratory of Microbial Metabolism, School of Life Sciences and Biotechnology, Shanghai Jiao Tong University, Shanghai, China, 2 State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, Shanghai, China, 3 Institute of Oceanology, Shanghai Jiao Tong University, Shanghai, China*

#### *Edited by:*

*Roy D. Sleator, Cork Institute of Technology, Ireland*

#### *Reviewed by:*

*Huiluo Cao, The University of Hong Kong, Hong Kong Pat G. Casey, University College Cork, Ireland*

> *\*Correspondence: Xiang Xiao<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Evolutionary and Genomic Microbiology, a section of the journal Frontiers in Microbiology*

*Received: 12 May 2015 Accepted: 26 October 2015 Published: 10 November 2015*

#### *Citation:*

*He Y, Feng X, Fang J, Zhang Y and Xiao X (2015) Metagenome and Metatranscriptome Revealed a Highly Active and Intensive Sulfur Cycle in an Oil-Immersed Hydrothermal Chimney in Guaymas Basin. Front. Microbiol. 6:1236. doi: 10.3389/fmicb.2015.01236* The hydrothermal vent system is a typical chemosynthetic ecosystem in which microorganisms play essential roles in the geobiochemical cycling. Although it has been well-recognized that the inorganic sulfur compounds are abundant and actively converted through chemosynthetic pathways, the sulfur budget in a hydrothermal vent is poorly characterized due to the complexity of microbial sulfur cycling resulting from the numerous parties involved in the processes. In this study, we performed an integrated metagenomic and metatranscriptomic analysis on a chimney sample from Guaymas Basin to achieve a comprehensive study of each sulfur metabolic pathway and its hosting microorganisms and constructed the microbial sulfur cycle that occurs in the site. Our results clearly illustrated the stratified sulfur oxidation and sulfate reduction at the chimney wall. Besides, sulfur metabolizing is closely interacting with carbon cycles, especially the hydrocarbon degradation process in Guaymas Basin. This work supports that the internal sulfur cycling is intensive and the net sulfur budget is low in the hydrothermal ecosystem.

Keywords: hydrothermal vent, metagenomics, metatranscriptomics, sulfur cycle, carbon cycle

## INTRODUCTION

Hydrothermal vents are often discovered in ocean ridges where hydrothermal fluid is emitted after the hydrothermal circulation and alteration of seawater entrained through geothermally heated subseafloor basalt (Von Damm, 1990). The deep-sea hydrothermal vent fluid is commonly characterized by its high temperature, varied salinity, enriched metallic elements, and particularly high contents of reduced chemicals, such as H2, CH4, and H2S (Jannasch and Mottl, 1985). A thermodynamic non-equilibrium is created when the hydrothermal vent fluid encounters sea water that is cold and at a rather high oxidative state, which allows various abiotic and biotic reactions occur. Thus, the hydrothermal vent system is a typical chemosynthetic ecosystem in which microorganisms play essential roles in the generation, consumption, and modification of energy available in the environment (Reysenbach and Shock, 2002).

In the hydrothermal vent ecosystem, almost all types of inorganic sulfur compounds (e.g., S2−, S, S2O2 2−, SO2, S2O3 2−, and SO4 2−) are abundant and actively converted through chemosynthetic pathways to provide energy and thus sustain the microbial population in the ecosystem (Nakagawa et al., 2005). For example, in the Lost City hydrothermal field, the dominant *Thiomicrospira*-like group, which consists of sulfuroxidizing chemolithoautotrophs, was observed in the carbonate chimney (Brazelton and Baross, 2010). In the Lau Basin hydrothermal vent field, sulfur-oxidizing Alphaproteobacteria, Gammaproteobacteria, and Epsilonproteobacteria have been suggested to be dominant in the exterior chimney, whereas putative sulfur-reducing Deltaproteobacteria are dominant in the interior of the chimney (Sylvan et al., 2013). In the Guaymas Basin hydrothermal vent field, sulfate-reducing microorganisms, e.g., Desulfobacterales, have been detected and are hypothesized to be involved in the anaerobic methane-oxidation process (Biddle et al., 2012). Moreover, the sulfur cycling is alternated by the chemical reactions that occur during the emitting and growth of the hydrothermal vent. Reduced sulfur compounds are extremely sensitive to oxidants and easily precipitated with metal ions to form chimney or nodule structures (Orcutt et al., 2011). Moreover, shifts in temperature and fluid composition have been observed during the life span of a hydrothermal vent. For example, at 9◦N East Pacific Rise, Bio9 vent fluids were 368◦C in 1991, increased to an estimated temperature greater than or equal to 388◦C after a second volcanic event in 1992, and thereafter declined over the next similar to 2 years reaching a temperature of 365◦C in December 1993 (Fornari et al., 1998). The hydrogen concentration in the hydrothermal plum in the NE Lau Basin dropped from 14843 nM in 2008 to 4410 nM in 2010 then further to 7 nM in 2012 (Baumberger et al., 2014). As a result, environmental fluctuations may be induced between sulfateand sulfur-reducing archaea and contribute to the diverse roles of these microorganisms in the ecosystem (Teske et al., 2014). Therefore, a better understanding of sulfur cycling is essential for describing the geobiochemistry and providing hints to identify the life status of a hydrothermal vent ecosystem.

Due to the complexity of microbial sulfur cycling resulting from the numerous parties involved in the process, the sulfur budget in a hydrothermal vent is poorly characterized. To date, most studies have focused on the abundance and diversity of sulfur oxidizers and sulfate reducers in environmental samples through a metagenomic approach (Nakagawa et al., 2005). The exception is the study conducted by Anantharaman et al. (2013), who combined metatranscriptomic and metagenomic analyses of a hydrothermal plume sample and demonstrated the novel metabolic potentials of the SUP05 group of uncultured sulfur-oxidizing Gammaproteobacteria. However, this finding is based on the near-complete genomes of two SUP05 populations, and the information is restricted to this particular group of sulfur oxidizers (Anantharaman et al., 2013). The in-depth mining of the metatranscriptomic data remains too scarce to allow construction of the entire sulfur cycle and thus further illustrate the interactions of this process with the biological cycling of C, N, and O elements.

The Guaymas Basin in the Gulf of California is a young marginal rift basin characterized by the active hot venting of reduced sulfur compounds and the rapid deposition of organic-rich sediments. These features make the sulfur cycle in this ecosystem particularly intensive and closely interact with the carbon cycle, including hydrocarbon degradation (Bergmann et al., 2011). Thus, this sampling site is ideal for illustrating all of the possible microbial sulfur metabolic pathways and to evaluate the maximal biomass contribution of sulfur-metabolizing microorganism to the hydrothermal vent ecosystem. In this study, we performed an integrated metagenomic and metatranscriptomic analysis on a chimney sample from Guaymas Basin to achieve a comprehensive study of each sulfur metabolic pathway and its hosting microorganisms and constructed the microbial sulfur cycle that occurs in the site.

### RESULTS

### Composition of the Microbial Community

The composition and function of this microbial community were assessed at both the DNA and RNA levels to estimate the community metabolic potential and activity, respectively. The metagenome and metatranscriptome sequencing resulted in 199,903,215 and 1,885,022,958 bp clean sequences, respectively (**Table 1**). The metagenome raw reads were assembled into 49,055 contigs with an average length of 544 bp. In total, 5,417,253 reads (26.2%) from the metatranscriptome were mapped onto metagenomic contigs for quantification of the gene transcripts. 222 and 690,059 16S rRNA gene fragments were identified from the metagenome and metatranscriptome, respectively. The class-level taxonomic compositions of the metagenome and metatranscriptome revealed obvious differences in the presence and the activity of microbes in this community (**Table 1**). At the DNA level (**Figure 1A**), Archaeoglobi were found to be the most abundant, with 24.0% of the sequences assigned, and followed by Deltaproteobacteria (23.6%) and

Epsilonproteobacteria (11.3%). At the RNA level (**Figure 1B**), the same dominant groups were found: Deltaproteobacteria (31.8%), Archaeoglobi (13.3%), and Epsilonproteobacteria (12.8%). As reported previously (He et al., 2013), 53,034 gene features were predicted and then followed by manual examination and 19,491 gene features (36.8%) were considered to have expressions determined by transcriptomic reads mapping (see Materials and Methods). A total of 8929 (45.3%) and 4628 (23.7%) of all of the expressed genes were assigned (based on the BLAST results as described in Section "Materials and Methods") to Bacteria and Archaea, respectively, and the remaining sequences were not assigned to any category. Among the 13,557 expressed genes with taxonomic information, 2135 (15.7%) were from the highly abundant Archaeoglobi, which is consistent with the results from the 16S rRNA gene analysis. Although the assignment of bacterial genes could not be resolved well at the family level, the dominance of Deltaproteobacteria and Epsilonproteobacteria was still observed. As the archaeal cells typically have fewer copies of the 16S rRNA gene compared with bacterial cells, the proportion of active Archaeoglobi in this community was underestimated. Nevertheless, the predominant active players in this microbial community were Deltaproteobacteria, Archaeoglobi, and Epsilonproteobacteria.

The *de novo* assembly of metagenomic reads and binning by tetranucleotide signatures (Dick et al., 2009) identified three genomic bins (Supplementary Figure S1 and Supplementary Table S1). These three bins (herewith denoted bin20, bin21, and bin22) were assigned based on their phylogenomic marker genes to *Desulfobacteraceae*, Desulfovibrionales and *Archaeoglobus*. The identified genes in the obtained bins ranged from 486 to 1224. The genome completeness was estimated to range from ∼10 to 34%, based on singlecopy gene estimation (Supplementary Table S1). These three genomic bins will improve the taxonomic assignment of the expressed genes and the reconstruction of the metabolic pathways.

### Sulfur Metabolism

The genes involved in the oxidation of reduced sulfur (ORS) are sulfide quinone oxidoreductase (*sqr*), which mediates the oxidation of sulfide (HS−) to elemental sulfur (S0), the Sox enzyme complex (*soxABXYZ*), which is responsible for the oxidation of thiosulfate (S2O3 2−) to elemental sulfur, the reverse dissimilatory sulfite reductase complex (*rdsr*), which is responsible for the oxidation of elemental sulfur to sulfite (SO3 2−), and adenosine 5 -phosphosulfate reductase (*apr*) and sulfate adenylyltransferase (*sat*) for oxidation of sulfite to sulfate (SO4 2−; Anantharaman et al., 2013). Conversely, the genes associated with the dissimilatory sulfate reduction (DSR) pathway (Fritz et al., 2002) are *sat*, *apr,* and sulfite reductase (*dsr*). The repertoire of genes associated with the ORS and DSR pathways were found to be expressed in this community (**Table 2**). Both *apr* and *dsr* were found at high expression levels in bin21 and bin22, confirming their active presence in SRB and *Archaeoglobus*. The *sqr* gene, key gene in the ORS pathway, is found present and active in Epsilonproteobacteria, of which the most highly expressed representative was classified into *Sulfurimonas* (**Figure 2**) that is one of the most abundant sulfur-oxidizing bacteria found in hydrothermal vent chimneys (Cao et al., 2014). The *sox* genes were not identified in either the metagenome or metatranscriptome (**Table 2**). In Epsilonproteobacteria, the proposed microorganism in the present study to perform the ORS pathway, *sat* gene was found to exhibit high and medium expression levels (**Table 2**). However, either *aprAB* or *dsrAB* was identified in the metagenome or metatranscriptome. This finding may be due to the fact that the 454-based metagenomes are still with low coverage and unable to present all the important functional genes. In

Deltaproteobacteria and Archaeoglobales, which were proposed to conduct the DSR pathway in this study, *aprAB* and *dsrA* genes were found to be highly expressed in both of these two taxonomic groups, whereas *sat* and *dsrB* genes were found only in Deltaproteobacteria. The phylogenies of *aprA* and *dsrA* further confirmed their assignment to Deltaproteobacteria (Supplementary Figures S2A,B). In a previous study, the *aprA* with the highest abundance was assigned to the genus *Desulfobulbus* (Cao et al., 2014). In our study, the *aprA* gene with the highest expression was assigned to *Desulfovibrio* (Supplementary Figure S2A). To summarize, the taxonomic assignment and expression of key genes in the sulfur cycle suggest that both the ORS and DSR pathways are highly active in this oil-immersed microbial community, and the energy generated by the sulfur metabolism supports the dominant and active group (**Figure 3**).

Because there are no metatranscriptome published for any hydrothermal vent chimneys, we compared the expression patterns of the sulfur-metabolizing genes in this metatranscriptome to those in the available metatranscriptome of a plum sample that was also collected from Guaymas Basin (Lesniewski et al., 2012). As shown in **Figure 4**, sulfur metabolizing (including oxidation and reduction) genes were among the most abundant genes found in the metatranscriptome, and a significant difference (*p*value *<* 0.001) in the expression profiles of sulfur metabolizing genes was observed between the chimney and the plume metatranscriptome. Therefore, the sulfur-metabolizing genes were highly abundant and expressed in this GB chimney sample, and displayed significantly higher expression pattern than

those of a hydrothermal vent plume sample from Guaymas Basin.

### Carbon Metabolism

In this study, the complete WL pathway was identified in Archaeoglobales with high expression levels (Supplementary Table S2). The CBB cycle was not identified. The genes involved in the complete rTCA cycle were found to be actively present in both Deltaproteobacteria and Epsilonproteobacteria that dominated this chimney microbial community (**Table 3**). The key gene in the rTCA cycle, ATP-citrate lyase (*acl*), identified in this study to exhibit the highest expression was from Epsilonproteobacteria and exhibited the highest similarity to *Sulfurovum*, a novel sulfur-, nitrate-, and thiosulfate-reducing and strictly anaerobic chemolithoautotroph bacterium isolated from a deep-sea hydrothermal vent chimney at the Central Indian Ridge (Mino et al., 2014). In this study, the key enzyme for the utilization of acetate, acetyl-CoA synthetase (*acd/acs*), was found to be expressed and was assigned to sulfate-reducing bacteria (SRB; bin21 as shown in **Table 3**). In addition, the rTCA cycle and WL pathway were found to be the main pathways for carbon fixation by the dominant Bacteria and Archaea, respectively. This result suggests that, in combination with sulfur metabolism, autotrophic carbon fixation may play an important role in the survival and dominance of these species in the community. Moreover, as shown in Supplementary Table S3, genes involved in the flagellar assembly process were found to be actively present in Desulfovibrionales (bin21). The active role of the flagellar system in SRB may facilitate the movement toward electron donors and nutrients that occurs under the highly fluctuating conditions resulting from eruptions of hydrothermal vents. SRB have been reported to have the potential to anaerobically oxidize diverse hydrocarbons, such as alkanes, in Guaymas Basin sediments and chimney samples (Rueter et al., 1994). In this study, the activity and expression level of the presumably key gene in fumarate addition, a process through which alkanes are added to the double bond of fumarate based on the activity of alkylsuccinate synthase (*ass*), was checked. The *ass* genes were found to be highly active in this community, as determined through their expression level, and the most highly expressed hits were from *Desulfoglaeba alkanexedens* (Agrawal and Gieg, 2013), a typical sulfate-reducing and alkane-oxidizing bacterium (Supplementary Table S4). Moreover, the enzymes required for the degradation of a variety of organic compounds, such as hydrocarbons, fatty acids, chitins and proteins, have been detected in both the metagenome and metatranscriptome (Supplementary Table S5). Despite their important roles in carbon and global sulfur cycle, the energy metabolism of SRB remains poorly understood. After taxonomic assignment (see Materials and Methods), cyctochrome c (*cytC*), formate dehydrogenase (*fdh*), F-type ATPase (*atp*), NADH-quinone oxidoreductase (*nuo*), electron transport complex protein (*rnf*) and hydrogenases, such as Ni/Fe-hydrogenase I (*hyaAB*) and hydrogenase nickel incorporation and accessory protein (*hypA* and *hypB*), were found with expressions and assigned to SRB (Supplementary Table S6). The presence of hydrogenases and *fdh*

fragments per kilobase of transcript per million fragments mapped) are displayed and discussed. The processes conducted by Archaea are shown on the left, whereas those conducted by bacteria are presented on the right. Detailed information of these genes is displayed in Tables 2–4.

abundance of the gene transcripts was normalized to the length of the gene fragment and the total number of all of the transcripts.

may suggest that H2 or formate and play important roles in the flow of electrons during sulfate reduction. As shown above, the sulfur cycle in this community was particularly intensive and closely interacted with the carbon cycle, including carbon fixation and hydrocarbon degradation, to sustain the primary production in this ecosystem.

#### TABLE 3 | Genes identified in the rTCA pathway in Delta- and Epsilonproteobacteria species.

∗*The taxonomy assignments were determined by two methods, as described in Section "Materials and Methods." The binning index is explained in Supplementary Table S1. #FPKM is based on the maximal expression value of the annotated genes.*

### Nitrogen Metabolism

The key genes involved in the nitrogen metabolism were found, and some of these were found to be actively expressed (**Table 4**). Many Bacteria and Archaea have the potential to perform denitrification (Philippot, 2002), and numerous organic and inorganic compounds can be used as electron donors for denitrification. The genes involved in denitrification, including *nar* (nitrate reductase), *nap* (nitrate reductase), *nir* (nitrite reductase), *nor* (nitric oxide reductase) and *nosZ*, were found to be present in the metagenome. The *narG* gene was assigned to *Beggiatoa*, a nitrate-respiring and sulfide-oxidizing bacterium that has been found to dominate microbial mats in hydrothermal sediments in the Guaymas Basin (Winkel et al., 2014). *narJ* was found to be expressed in Alteromonadales, whereas *napA* and *napB* were found to be expressed in Epsilonproteobacteria. To summarize, a complete set of denitrification genes were found in the bacterial community of the chimney, though some of them were found at low expression levels (**Table 4**). Based on this observation, we propose that nitrogen denitrification present in this community is most likely mediated by Gammaproteobacteria and Epsilonproteobacteria, with electrons generated by the ORS pathway.

### DISCUSSION

Since the discovery of the deep-sea hydrothermal ecosystem in 1977, it has been proposed that hydrogen sulfide-oxidizing chemoautotrophs may potentially sustain the primary production in these ecosystems (Kvenvolden et al., 1995), where hydrogen sulfide or sulfide is primarily supplied via the high temperatures of seawater-rock interactions in the subseafloor hydrothermal reaction zones (Jannasch and Mottl, 1985). The chemical and microbial oxidation and reduction reactions of sulfur compounds probably establish the overall sulfur metabolism in the ecosystem (Yamamoto and Takai, 2011). There is no doubt that the sulfur cycle is one of the most important microbial chemosynthetic pathways in the microbial habitats of hydrothermal vents, but few studies have attempted to characterize the process, particularly at the function and activity levels. To date, the mechanism through which a microbial community in hydrothermal fields can be fueled by sulfate metabolism remains unclear. In particular, metagenomic approaches have not been widely applied in studies of energy generation by the microbial sulfur cycle in hydrothermal systems. In this study, a combined metagenomic

∗*The taxonomy assignments were determined by two methods, as described in Section "Materials and Methods." The binning index is explained in Supplementary Table S1. #FPKM is based on the maximal expression value of the annotated genes.*

and metatranscriptomic study of a chimney in the Guaymas Basin provides insight into the complete sulfur cycle based on the results from not only the genomic but also the expression analysis, the combination of which has not been previously used for the analysis of a deep-sea hydrothermal vent chimney sample.

The accumulation of hydrogen sulfides at the outer chimney promoted the coupling of sulfide oxidation to the electron acceptors present in the nearby marine water, including oxygen and nitrate, as supported by the retrieval of the functional and expressed genes described herein (**Tables 2–4** and **Figure 3**). These findings suggest that the coupling between sulfur oxidation and denitrification may fuel some N-metabolizing microorganisms at the sulfide-enriched outer chimney. As proposed in this study, the microorganisms involved in this process were Epsilonproteobacteria as the sulfur-oxidizing bacteria, and Gammaproteobacteria and Epsilonproteobacteria as potential denitrifiers. The other sulfur-metabolizing group, namely sulfate-reducing prokaryotes, may use hydrogen and/or dissolved organic matter as electron donors, as hydrogenases and key genes for the degradation of organic compounds have been identified in this study (Supplementary Tables S5 and S6).

Carbon fixation pathways other than the Calvin–Benson– Bassham (CBB) cycle have been found to exhibit a notable contribution to carbon fixation, mostly at deep-sea hydrothermal vents (Campbell and Cary, 2004). The rTCA cycle was found to be highly expressed in the dominant Delta- and Epsilonproteobacteria. The key enzyme for the utilization of acetate was also identified to be expressed in this study (**Table 3**). Generally, the rTCA cycle appears to be dominant in habitats with a temperature ranging from 20 to 90◦C, whereas the CBB cycle and the Wood-Ljungdahl (WL) pathway may be the principal pathways at temperatures lower than 20◦C and greater than 90◦C, respectively (Hugler and Sievert, 2011). In the present sample, the CBB cycle was not found present, which is consistent with the fact that this sample was collected from a high-temperature condition (He et al., 2013). In addition, the enzymes for the degradation of a variety of organic compounds, such as hydrocarbons, fatty acids, chitins and proteins, have been detected at both DNA and RNA level (Supplementary Table S5). Together, all of these organic compounds may be the carbon source for this microbial community.

In this scenario, both autotrophic and heterotrophic SRB could inhabit the inner chimney (**Figure 3**), where sulfate reduction is coupled to carbon fixation and hydrocarbon oxidation. Based on the expression levels of key genes in rTCA (**Table 3**) and alkane degradation (Supplementary Table S4), hydrocarbon degradation might contribute substantially to the linking of S and C cycle at inner layer chimney. In another word, heterotrophic SRB, commonly found at vent systems, may be the major player in coordinating and influencing the S and C cycle. Compared the expression of key genes in sulfur metabolizing and the rest processes (**Figure 4**), the reduced sulfur would be quickly and intensively oxidized to fuel the community, where sulfate-reducing microbes were found dominated. The composition of the sulfate-reducing community was determined by the way that microbes perform carbon metabolism. In our sample, heterotrophic SRB was found prevalent with their capabilities in hydrocarbon degradation. This finding may improve our understanding on the structure, function, and interaction within microbial community in hydrothermal vent.

Meta-omics based approaches have the advantages in studying the entire microbial community without pure cultures or prior knowledge on the sample. Functional omics approaches, such as transcriptome and proteome, could further confirm the metabolic potential at the active level. More efforts will be spent on quantification and comparison of these function omics datasets. Together with *in situ* carbon stable isotope measurement, and lipid type and diversity analysis, the activity, rate and interaction of key process in a given environmental condition could be accessed and estimated.

### MATERIALS AND METHODS

### Sample Collection and Processing

The sample 4558-6 under investigation was collected from the outer layer of a black-smoker chimney in the Guaymas Basin and was previously described through a metagenome-based study (He et al., 2013). The sample was fixed with RNAlater (Sigma-Aldrich, Munich, Germany) and stored at −80◦C prior to DNA and RNA extraction. DNA isolation was conducted as described previously (Wang et al., 2013). Metagenome pyrosequencing was performed using a 454 Life Sciences GS FLX system with a practical limit of 400 bp. RNA was isolated with a RNA isolation kit (Omega Bio-Tek, Doraville, GA, USA) following the user's manual provided by the manufacturer. RNA samples were treated with DNAse (Thermo) for 45 min at 37◦C, and then used as a template for PCR to detect undigested DNA. The mRNA fraction was enriched through the enzymatic digestion of rRNA molecules (mRNA-ONLY Prokaryotic mRNA Isolation kit, Epicentre Biotechnologies, Madison, WI, USA) followed by the subtractive hybridization of rRNA with capture oligonucleotides (Ambion MICROBExpress kit, Life Technologies, Gaithersburg, MD, USA). The mRNA isolates were first amplified (MessageAmp II-Bacteria kit, Ambion, Life Technologies) and then reversely transcribed into complementary DNA. Afterward, the cDNA was directly sequenced using the Illumina (BGI-Shenzhen, China) Hiseq2000 platform (2∗90 bp pair-end) for metatranscriptome analysis.

### Metagenome Assembly and Annotation

The reads obtained through metagenome sequencing were assembled and annotated as previously described (He et al., 2013). Briefly, low quality sequencing reads were trimmed in Geneious 6.04 (Biomatters Ltd.) and technical replicates were removed with cd-hit (at 96% sequence identity; Fu et al., 2012). After removing short reads (*<*100 bp), the remaining reads were assembled with Velvet (Zerbino and Birney, 2008). Coding regions of the metagenomic assembly were predicted using FragGeneScan (Rho et al., 2010) and then BLASTed (Altschul et al., 1997; 1e−5) against an NCBI non-redundant (NR) protein database. The 16S rRNA genes were picked using Sortmerna and BLASTed against GreenGene database (e-value *<* 1e−5) respectively. For functional annotation, sequences with matches to the COG (Tatusov et al., 2003), Pfam (Finn et al., 2014), and KEGG (Ogata et al., 1999) databases were retrieved to establish the functional categories and reconstruct the metabolic pathways. The genes of interest, such as transposases, were subjected to manual checkup, and spurious annotations (putative, like-, similar to) were excluded from further analysis.

### Taxonomic Assignment

Two different methods were applied to assess the taxonomic information. First, the assembled metagenomic sequences was binned using the tetranucleotide frequencies in emergent self-organizing maps (ESOMs; Dick et al., 2009) with a window size of 8 kbp, a sliding window size of 4 kbp, and the minimum fragment size of 2 kbp. Complete genomic sequences of 20 species were used as references (designated as bin1–20), these microorganism were listed as following: *Acinetobacter pittii* ANC 4052, *Alteromonas macleodii* str. 'Deep ecotype,' *Candidatus Pelagibacter ubique* HTCC1062, uncultured marine crenarchaeote E37-7F, Marine group II euryarchaeote SCGC AAA288-C18, Marine Group II euryarchaeote SCGC AB-629-J06, uncultured marine group II euryarchaeote (marine metagenome), Marine Group III euryarchaeote SCGC AAA007-O11, Marine Group III euryarchaeote SCGC AAA288-E19, *Marinobacter nanhaiticus* D15-8W, *Methylobacter tundripaludum* SV96, *Methylophaga aminisulfidivorans* MP, *Methylotenera mobilis* JLW8, *Nitrosopumilus maritimus* SCM1, *Candidatus Nitrospira defluvii*, *Planctopirus limnophila* DSM 3776, *Pseudomonas denitrificans* ATCC 13867, *Candidatus Ruthia magnifica* str. Cm (Calyptogena magnifica), SAR324 cluster bacterium SCGC AAA240-J09 and SAR86 cluster bacterium SAR86E. After binning, the completeness and taxonomic classification of the genomes within bins were then estimated by counting and BLASTing universal single-copy genes as previously described (Rinke et al., 2013). Alternatively, each predicted sequence feature in the metagenome and metatranscriptome was assigned to a certain taxon if at least 75% of the BLAST hits of this query were from that specific taxon. A BLAST search of all of the reads against the non-redundant protein database in NR was performed. All of the hits obtained from the BLAST searches were retained, and their taxonomic affiliations were determined using MEGAN (Huson et al., 2007) with bit-score values of 100. The taxonomic compositions of each predicted gene feature was then visualized using MEGAN.

### Metatranscriptome Mapping and Transcript Quantification

The raw shotgun sequencing metatranscriptomic reads obtained by Illumina pair-end sequencing were dereplicated (100% identity over 100% lengths) and trimmed using sickle1 . The dereplicated, trimmed, and paired-end Illumina reads were then mapped to the metagenome using Bowtie (Langmead and Salzberg, 2012) with the default parameters. The unique mapped reads were selected, and FPKM (expected fragments per kilobase of transcript per million fragments mapped) was used to estimate the expression level of each gene using a script downloaded from GitHub2 .

### Estimation of the Completeness of Genomic Bins

The complete genome sizes of the genomic bins were estimated based on an analysis of conserved single-copy genes (CSCGs) as described by Lloyd et al. (2013). In total, we were able to collect 162 and 139 universal CSCGs for the archaea and bacteria genomes, as in the previous study (Rinke et al., 2013). The ratios between the numbers of CSCGs present in the metagenome and the number of total CSCGs were then used to estimate the size of each genome bin.

1https://github*.*com/najoshi/sickle

2https://github*.*com/minillinim/sam2FPKG

### Comparative Analysis

The expression patterns of the sulfur-metabolizing genes in this metatranscriptome were compared to those in the metatranscriptome of a plum sample from Guaymas Basin (Lesniewski et al., 2012). Comparisons between two metatranscriptomes were conducted using the Mann–Whitney *U*-test. The gene expression profiles were compared between two samples using the normalized rank from 0 to 1 in each respective sample as the input. A difference was considered significant if the *p*-value was lower than 0.001.

### Construction of a Phylogenetic Tree

The predicted sequence features were checked across multiple annotation databases and then aligned with ClustalW (Larkin et al., 2007), and any gaps were removed manually. To construct functional gene phylogenies, the aligned sequences were analyzed by maximum likelihood-based FastTree (Price et al., 2010) using the Jones–Taylor–Thornton (JTT) with CAT approximation.

### Metabolic Pathway Identification

The gene products were searched for similarity against the KEGG database. A match was counted if the similarity search resulted in an expectation e-value below 1e−5. All of the occurring KO (KEGG Orthology) numbers were mapped against the KEGG pathway functional hierarchies and the COG database. For genes

### REFERENCES

with multiple hits, only the genes with the highest expression value (FPKM) are displayed in the figures and tables and further discussed in the text.

### Data Availability

The metatranscriptome sequences are available on NCBI as SRX1008212. The assembled sequence was uploaded to IMG with a project ID Ga0072503.

### ACKNOWLEDGMENTS

We thank Anna-Louise Reysenbach for providing the chance to attend the expedition, and all the crew members from AT-26 cruise. This work was supported by National High Technology Research and Development Program of China (Grant No. 2012AA092103), China Ocean Mineral Resources R & D Association (Grant No.DY125-22-04 and DY125- 15-T-04).

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal*.*frontiersin*.*org/article/10*.*3389/fmicb*.* 2015*.*01236

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 He, Feng, Fang, Zhang and Xiao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Degradation Network Reconstruction in Uric Acid and Ammonium Amendments in Oil-Degrading Marine Microcosms Guided by Metagenomic Data

*Rafael Bargiela1, Christoph Gertler2†, Mirko Magagnini3, Francesca Mapelli4, Jianwei Chen5, Daniele Daffonchio4,6, Peter N. Golyshin2\* and Manuel Ferrer1\**

#### *Edited by:*

*Eamonn P. Culligan, Cork Institute of Technology, Ireland*

### *Reviewed by:*

*Romy Chakraborty, Lawrence Berkeley National Lab, USA Efthymios Ladoukakis, National Technical University of Athens, Greece*

> *\*Correspondence: Manuel Ferrer<EMAIL_ADDRESS>Peter N. Golyshin<EMAIL_ADDRESS>
#### *†Present address:*

*Christoph Gertler, Friedrich Loeffler Institute– Federal Research Institute for Animal Health, Institute for Novel and Emerging Infectious Diseases, 17493 Greifswald, Germany*

#### *Specialty section:*

*This article was submitted to Evolutionary and Genomic Microbiology, a section of the journal Frontiers in Microbiology*

*Received: 06 July 2015 Accepted: 30 October 2015 Published: 24 November 2015*

#### *Citation:*

*Bargiela R, Gertler C, Magagnini M, Mapelli F, Chen J, Daffonchio D, Golyshin PN and Ferrer M (2015) Degradation Network Reconstruction in Uric Acid and Ammonium Amendments in Oil-Degrading Marine Microcosms Guided by Metagenomic Data. Front. Microbiol. 6:1270. doi: 10.3389/fmicb.2015.01270*

*1 Systems Biotechnology, Department of Biocatalysis, Institute of Catalysis, Consejo Superior de Investigaciones Científicas, Madrid, Spain, 2 School of Biological Sciences, Bangor University, Bangor, UK, 3 EcoTechSystems Ltd., Ancona, Italy, 4 Department of Food, Environmental and Nutritional Sciences, University of Milan, Milan, Italy, 5 Beijing Genomics Institute, Shenzhen, China, 6 Biological and Environmental Science and Engineering Division, King Abdullah University of Science and Technology, Thuwal, Saudi Arabia*

Biostimulation with different nitrogen sources is often regarded as a strategy of choice in combating oil spills in marine environments. Such environments are typically depleted in nitrogen, therefore limiting the balanced microbial utilization of carbon-rich petroleum constituents. It is fundamental, yet only scarcely accounted for, to analyze the catabolic consequences of application of biostimulants. Here, we examined such alterations in enrichment microcosms using sediments from chronically crude oil-contaminated marine sediment at Ancona harbor (Italy) amended with natural fertilizer, uric acid (UA), or ammonium (AMM). We applied the web-based AromaDeg resource using as query Illumina HiSeq meta-sequences (UA: 27,893 open reading frames; AMM: 32,180) to identify potential catabolic differences. A total of 45 (for UA) and 65 (AMM) gene sequences encoding key catabolic enzymes matched AromaDeg, and their participation in aromatic degradation reactions could be unambiguously suggested. Genomic signatures for the degradation of aromatics such as 2-chlorobenzoate, indole-3-acetate, biphenyl, gentisate, quinoline and phenanthrene were common for both microcosms. However, those for the degradation of orcinol, ibuprofen, phenylpropionate, homoprotocatechuate and benzene (in UA) and 4-aminobenzenesulfonate, *p*-cumate, dibenzofuran and phthalate (in AMM), were selectively enriched. Experimental validation was conducted and good agreement with predictions was observed. This suggests certain discrepancies in action of these biostimulants on the genomic content of the initial microbial community for the catabolism of petroleum constituents or aromatics pollutants. In both cases, the emerging microbial communities were phylogenetically highly similar and were composed by very same proteobacterial families. However, examination of taxonomic assignments further revealed different catabolic pathway organization at the organismal level, which should be considered for designing oil spill mitigation strategies in the sea.

Keywords: ammonium, biostimulation, crude oil degradation, enrichment, Mediterranean Sea, metagenomics, microcosm, uric acid

## INTRODUCTION

Oil pollution still is a global problem (Yakimov et al., 2007; Bargiela et al., 2015). At present, in many sea regions containment and recovery of oil using booms and skimmers is the method of choice for oil spill first responders (Walther, 2014). Especially in the open sea, the use of dispersants in combination with biostimulation and bioaugmentation agents based on non-toxic, natural low cost formulations, is encouraged, although the majority of tests have been performed at lab-scale (Das and Chandran, 2010; Nikolopoulou and Kalogerakis, 2010; Alvarez et al., 2011; Nikolopoulou et al., 2013). In marine systems, the low concentration of nitrogen, phosphorous, and oxygen, together with their low bioavailability are main factors limiting the degradation of carbon-rich hydrophobic compounds (Howarth and Marino, 2006; Venosa et al., 2010; Ly et al., 2014). Attempts have been made to use different nitrogen sources to promote the growth and selection of different microbial strains with greater catabolic capacity for combating oil spills compared to natural attenuation (Teramoto et al., 2009; Venosa et al., 2010). However, crude oil biodegradation requires about 0.04 g of nitrogen per gram of oil (Atlas, 1981) which makes the choice of nitrogen source pivotal for the whole treatment. Recent data highlighted the possible link between N cycling processes and hydrocarbon degradation in marine sediments (Scott et al., 2014). Therefore, it is essential to select appropriated N-containing biostimulants.

The sources of nitrogen for the degradation tests – mostly performed at lab-scale and in minor occasions at field-scale – included nitrate, ammonium (AMM), urea, uric acid (UA), amino acids and the hydrophobic substance lecithin (Garcia-Blanco et al., 2007; Li et al., 2007; Martínez-Pascual et al., 2010; Venosa et al., 2010; Nikolopoulou et al., 2013; Mohseni-Bandpi et al., 2014). Slow-release nitrogen (AMM-based) fertilizers have also been successfully used for growth stimulation in microbial oil remediation (Miyasaka et al., 2006; Teramoto et al., 2009; Reis et al., 2013). However, AMM has been proved ineffective in treatment of real oil spill due to co-precipitation with phosphates in seawater. In a recent study, we have shown that biodegradable natural fertilizers like UA can be used as costefficient biostimulant for enhancing bacterial growth in polluted sediments (Gertler et al., 2015). Each nitrogen source has its advantages and disadvantages, yet overall results have shown that the microbial populations were initially different from those found in the absence of biostimulants and that the degradation efficiency generally increased. It is therefore critical to establish how the whole microbial biodegradation network is affected and whether different pollutants are preferentially degraded as a consequence of amendments of biostimulants.

In an early work using the recently developed AromaDeg analysis (Duarte et al., 2014) and a meta-network graphical approach, we reconstructed the catabolic networks associated to microbial communities in a number of chronically polluted sites (Bargiela et al., 2015). The approach focuses on the usage of metagenomic data, which directly leads to a network that included catabolic reactions associated to genes encoding enzymes annotated in the genomes of the community organisms. We found key catabolic variations associated to changes in community structure and environmental constraints (Bargiela et al., 2015). In this work, this approach was applied to draft the catabolic networks of two different enrichment microcosms set up with sediments from chronically crude oil-contaminated marine sediments from Ancona harbor (Italy) and the natural fertilizer UA or AMM as nitrogen sources (Gertler et al., 2015). Ancona harbor is very close to the urban area and hosts a multi-purpose port receiving cruise liners, passenger ferries, commercial liners and fishing boats. A minor part of the related airborne pollutants is due to the vessels calling at the port while the main contribution comes from road traffic and other human activities. Furthermore, sediments in Ancona harbor are heavily contaminated due to its role as a major ferry terminal and industrial port on the Adriatic Sea. We hypothesize that the microbial community shifts previously observed after addition of UA and AMM (Gertler et al., 2015) may have an influence in the selection of certain catabolic pathways. Potential protein-coding genes (≥20 amino acids long) obtained by direct Illumina HiSeq sequencing of DNA material of the corresponding microcosms (Gertler et al., 2015) constituted the input information in our study.

### MATERIALS AND METHODS

### Study Site, Microcosm Set-up and Sequence Accession Numbers

The starting point of this study were the meta-sequences previously obtained by direct sequencing from two microcosm sets created using sediment samples from the harbor of Ancona (Italy; 43◦37 N, 13◦30 15E), as described previously (Gertler et al., 2015). Both microcosm setups were identical in size, composition, incubation, sampling regime and nutrient concentration with exception of the type of nitrogen source applied. Either AMM or UA were supplied in equimolar amounts of nitrogen. Briefly, one-liter Erlenmeyer flasks (duplicates) were filled with 150 g of sand (Sigma–Aldrich, St. Louis, MO, USA), sterilized and spiked with 10 mL of sterile filtered Arabian light crude oil. One gram of sediment from the sampling site was mixed into the oil-spiked sand as the inoculum. Three hundred milliliters of modified ONR7a medium (Dyksterhouse et al., 1995) (omitting AMM chloride and disodium hydrogen phosphate) was added. We added 5 mL of Arabian light crude oil, which based upon average literature values for density and molecular weight equals about 300 mM of C (Wang et al., 2003), 5 mM of NH4Cl and 0.5 mM of Na2HPO4 resulting in a molar N/P ratio of approximately 10:1. For UA treatment microcosm, 0.21 g (1.25 mmol = 5 mmol N) of UA was provided as nitrogen source while the AMM treatment microcosms were each supplied with 2.5 mL of a 2 M AMM chloride solution (5 mmol; pH 7.8). Both treatments also contained 2.5 mL of a 0.2 M disodium hydrogen phosphate solution (0.5 mmol; pH 7.8). Excess amounts of crude oil were added to compensate for the 35% carbon losses due to evaporation of volatile hydrocarbons over the course of the experiment. Including losses due to evaporation, the C/N/P ratio was approximately 400:10:1. Control treatments were set up: (i) a negative control contained only sterile sand and ONR7a; (ii) two further controls contained sand, ONR7a, crude oil and either UA or AMM chloride solution but no sediment sample; and (iii) one control contained oil, sterile sand, ONR7a medium and a sediment sample, but no additional nitrogen source or phosphorus source was provided. No significant growth was detected under tested control conditions. Under the given assay conditions, the utilization of UA as carbon source is minimal, as the amount of carbon introduced by UA into the microcosms was disproportionately low in contrast to the residual carbon in the sediment and the carbon introduced in form of oil. Briefly, we added 300 mmols of carbon in form of oil and only 6.25 mmols of carbon in form of UA. In addition, the molar ratio C/N in the system (between 10:1 and 40:1, depending UA or AMM was added) implies there was excess of carbon in the medium and thus the growth was limited by N.

The resulting microbial communities from microcosms were destructively sampled after 21 days of incubation at 20◦C, the isolated DNA subjected to the paired-end sequencing (Illumina HiSeq 2000) at Beijing Genomics Institute (BGI; China), and gene calling performed as described (Gertler et al., 2015). Taxonomic affiliations of potential protein-coding genes were predicted as described previously (Guazzaroni et al., 2013; Bargiela et al., 2015).

The meta-sequences are available at the National Center for Biotechnology Information (NCBI) with the IDs PRJNA222664 [for MGS-ANC(UA)] and PRJNA222663 [for MGS-ANC(AMM)]. The Whole Genome Shotgun projects are also available at DDBJ/EMBL/GenBank under the accession numbers AZIH00000000 [for MGS-ANC(UA)] and AZIK00000000 [for MGS-ANC(AMM)]. All original non-chimeric 16S small subunit rRNA hypervariable tag 454 sequences were archived at the EBI European Read Archive under accession number PRJEB5322. Note that the samples were named based on the code 'MGS', which refers to MetaGenome Source, followed by a short name indicating the origin of the sample and the nitrogen source, as follows: MGS-ANC(AMM) (the harbor of Ancona and AMM as nitrogen source); MGS-ANC(UA) (the harbor of Ancona and UA as nitrogen source).

### Biodegradation Network Reconstruction: Scripts and Commands for Graphics

The web-based AromaDeg resource (Duarte et al., 2014) was used for catabolic network reconstruction. AromaDeg is a web-based resource with an up-to-date and manually curated database that includes an associated query system which exploits phylogenomic analysis of the degradation of aromatic compounds. This database addresses systematic errors produced by standard methods of protein function prediction by improving the accuracy of functional classification of key genes, particularly those encoding proteins of aromatic compounds' degradation. In brief, each query sequence from a genome or metagenome [MGS-ANC(AMM) and MGS-ANC(UA), in this study] that matches a given protein family of AromaDeg is associated with an experimentally validated catabolic enzyme performing an aromatic compound degradation reaction. Individual reactions, and thus the corresponding substrate pollutants and intermediate degradation products, can be linked to reconstruct catabolic networks. We have recently designed an in-house script allowing the automatic reconstruction of such networks in a graphical format, which was used in present work. The script allows visualization and comparison of the abundance levels of genes encoding catabolic enzymes assigned to distinct degradation reactions as well as substrates or intermediates possibly degraded by distinct microbial communities. The complete workflow, including the scripts and commands used for catabolic network reconstruction has recently been reported (Bargiela et al., 2015).

Note that the sequence material used in the present investigation for biodegradation network reconstruction was based upon single biological microcosm replicate to preserve maximum coverage and sequencing depth as well as for other technical reasons, as described previously (Gertler et al., 2015). For each of the metagenome datasets the rarefaction curves of the observed species were estimated to analyze the species sampling coverage, and found that the rarefaction curves indicate closeness to saturation in each of the samples (Gertler et al., 2015). Therefore, with a single run of paired-end Illumina sequencing we determined populations that really represent the actual state of the microbial community in the microcosms and that biases were not introduced due to differences in microbial coverage. Whether or not more replicates may introduce some differences in the present study was not examined. However, because of the low standard deviation in the cultures (also checked for the representativeness of the microcosm by 16S small subunit rRNA hypervariable tag 454 sequences fingerprinting; Gertler et al., 2015) and the fact that sampled 16S rRNA diversity indicated closeness to saturation, we considered that the presented data are valid. Note that experimental validations (see Experimental Validations of Predicted Biodegradation Capacities) were performed in triplicates (with appropriated standard deviations), on the basis of which metagenome-based predictions were confirmed. Therefore, we considered that the differences at the taxonomic, gene content levels and catabolic capacities herein presented are most likely due to actual biological variability and are not random.

### Experimental Validations of Predicted Biodegradation Capacities

The ability of each of the microcosms to grow on pollutants expected to be degraded, was confirmed as follows. First, UA and AMM microcosms (in triplicates) were obtained as described above but omitting Arabian light crude oil; instead, a mix of pollutants containing naphthalene, 2,3-dihydroxybiphenyl, benzene, *p*-cumate, orcinol, 2-chlorobenzoate, phthalate and phenylpropionate, all from Fluka-Aldrich-Sigma Chemical Co. (St. Louis, MO, USA), was added at a final concentration of 2 ppm each. These pollutants were selected on the basis of existing analytical methods to quantify their concentrations (Bargiela et al., 2015). Control cultures without the addition of sediments but with chemicals and cultures plus sediments but without the addition of chemicals were set up.

The extent of degradation in test and control samples was quantified as follows. Briefly, bacterial cells (from 300 ml culture) were separated by centrifugation at 13,000 *g* at room temperature for 10 min. After supernatant separation, bacterial pellet was used for methanol extraction by adding 1.2 mL of cold (−80◦C) highperformance liquid chromatography (HPLC)-grade methanol. The samples were then vortex-mixed (for 10 s) and sonicated for 30 s (in a SonicatorR 3000; Misonix) at 15 W in an ice cooler (−20◦C). This protocol was repeated twice more with a 5-min storage at −20◦C between each cycle, and the final pellet was removed following centrifugation at 12,000 *g* for 10 min at 4◦C. Methanol solution was stored at −80◦C in 20 mL penicillin vials until they were analyzed by mass spectrometry and different and complementary separation techniques, namely liquid chromatography electrospray ionization quadrupole timeof-flight mass spectrometry (LC-ESI-QTOF-MS) in positive and negative mode, and gas chromatography-mass spectrometry (GC-MS), as described previously (Bargiela et al., 2015). The abundance levels of mass signatures of tested pollutants and key degradation intermediates, namely, salicylate, gentisate, catechol, benzoate and protocatechuate, were used as indicator of the presence of the corresponding enzymes encoded by catabolic genes.

### RESULTS AND DISCUSSION

### Bacterial Community Structures in Microcosms

A graphical approach recently described (Bargiela et al., 2015) was applied to draft the catabolic networks of two different oil-degrading marine microcosms. They were obtained from Ancona harbor sediments which were applied in a series of two enrichment microcosms, where AMM or UA were supplied to introduce equivalent amounts of nitrogen. Using partial 16S rRNA gene sequences obtained in the non-assembled Illumina reads through a metagenomic approach, it was firstly found a relatively high degree of similarity in the emerging communities (Gertler et al., 2015). Proteobacteria were the most abundant (AMM: 74.5%; UA: 74.2%, total sequences), in agreement with the fact that this bacterial group is the most abundant in other chronically crude oil-contaminated marine sediments within the Mediterranean Sea (Bargiela et al., 2015). Noticeably, all proteobacterial families were found in both microcosms (for details see **Table 1**). However, differences in the abundance of some community members could be observed on the basis of corresponding read frequency. As an example, the percentage of members of the *Rhodobacteraceae* and *Enterobacteriaceae* was elevated in microcosms supplied with AMM (18.2% AMM vs. 0.8% in UA and 5.6% in AMM vs. 3.2% in UA, correspondingly). Conversely, lower percentages of members of the *Alteromonadaceae* (9.6%/19.2%), *Halomonadaceae* (5.6%/7.8%), *Moraxellaceae* (0.5%/7.9%) and *Flavobacteriaceae* (1.8%/5.7%) could be detected in the AMM-supplemented microcosm in comparison to UA-based microcosms. At a genus level, 55 out of 57 identified proteobacterial taxa were common in both communities. However, enrichments containing AMM were characterized by higher percentages (referred to total reads) of Alphaproteobacteria, such as *Roseovarius* sp. (1.4% in AMM

*Results are based on the analysis of the partial 16S ribosomal RNA (rRNA) gene sequences extracted from non-assembled DNA sequences obtained by pairedend Illumina HiSeq 2000 sequencing.*

1*Only lineages with abundance of reads >1% are shown. Data from Gertler et al. (2015).*

vs. 0.1% in UA), *Ruegeria* spp. (1.1%/0.1%) and *Sulfitobacter* sp. (1.5%/0.1%), and some Gammaproteobacteria such as *Vibrio* sp. (2.4%/1.1%). In stark contrast to this, the UA-based enrichments showed significantly elevated percentages of members of the Firmicutes (7.4% in UA enrichments/4.9% in AMM enrichments) and Gammaproteobacteria, such as *Aeromonas* spp. (1.5%/0.5%) and *Pseudoalteromonas* sp. (1.9%/0.4%). Highly elevated percentages in UA enrichments were observed for the genera *Acinetobacter* (0.9% in UA enrichments/0.1% in AMM enrichments), *Halomonas* (6.1%/4.2%), *Marinobacter* (16.8%/6.3%) and *Psychrobacter* (6.8%/0.2%). A direct comparison of percentages of potentially oil degrading microbial genera in both microcosms showed a higher percentage of *Acinetobacter* sp. (0.9%/0.1%), *Idiomarina* sp. (0.8%/0.3%), *Oleiphilus* sp. (0.2%/0.03%) and *Marinobacter* sp. (16.8%/6.3%) but lower percentages of *Alcanivorax* sp. (3.9%/4.9%) and *Thalassolituus* sp. (0.04/0.8%) in the UA treatments (Gertler et al., 2015).

### Biodegradation Networks

As we were interested in obtaining networks that emphasized the catabolic differences within both microcosms, we selected a metagenomic approach to query the presumptive degradation capacities associated to both microcosms. The identification depends heavily on gene abundance and, despite the fact that a substantial faction of less abundant DNA in metagenomes remains undiscovered, the identified catabolic genes are assumed to represent the dominant presumptive pathways in each system. A rarefaction curve of the observed species for both samples to analyze species sampling coverage indicated closeness to saturation in each of the two microcosms (Gertler et al., 2015). In combination with the fact that both samples were sequenced to a similar extent (24,752,834 bp for AMM and 19,364,101 bp for UA; Gertler et al., 2015), this suggests that biases during the comparative analysis within the metagenomes were not introduced due to differences in microbial and sequence coverage.

Using as a query the 27,893 (for UA) and 32,180 (for AMM) potential protein-coding genes (for ≥20 amino acids-long polypeptides) (Gertler et al., 2015), we identified respectively a total of 45 (or 0.16% relative abundance in UA referred to the total number of protein-coding genes) and 65 (or 0.20% relative abundance in AMM) genes encoding catabolic enzymes with matches in AromaDeg (Duarte et al., 2014). This suggests that the biostimulants did not have much influence on the relative abundance of catabolic genes. However, significant differences can be observed when examining the diversity of genes encoding catabolic enzymes assigned to different families (**Figure 1**). The amount of genes encoding Rieske non-heme iron oxygenases and extradiol dioxygenases (EXDO) of the cupin superfamily increased 2- and 4-fold, respectively, and proved more abundant in the AMM microcosm in comparison to those in the UA microcosm (**Figure 1**).

The differences in family shifts may have an influence on degradation capacities provided by microorganisms in AMM and UA microcosms. To assess this, the presumptive aromatic degradation reactions and the substrate pollutants or intermediates possibly degraded by each of the two communities were predicted, and the corresponding degradation networks constructed (**Figure 2**). For that we used the AromaDeg web system that allows identifying catabolic genes and appropriated scripts and commands for graphics (for details see Biodegradation Network Reconstruction: Scripts and Commands for Graphics). Unambiguous reaction specificities could be detected for 35 (in UA) and 48 (in AMM) catabolic genes and were considered in the degradation network (**Figure 2**). However, no clear specificities could be assigned to 4 (in UA) and 11 (in AMM) Rieske oxygenases and 12 (six in UA and six in AMM) EXDO, which subsequently were not considered in the network. As shown in **Figure 2**, on the basis of the presence of genes encoding catabolic genes involved in particular transformations, the potential degradation of nine intermediates involved in the degradation of six key pollutants (2 chlorobenzoate, indole-3-acetate, biphenyl, gentisate, quinoline

FIGURE 1 | Number and diversity of sequences of gene families encoding key catabolic enzymes involved in the degradation of aromatic pollutants. Catabolic genes were identified as follows. Briefly, predicted open reading frames (ORFs) in the metagenomic DNA sequences were filtered by sequence homology (*>*50%) and minimum alignment length (*>* 50 amino acids) according to their similarity to the AromaDeg sequences of key aromatic catabolic gene families (and sub-families) involved in the degradation of aromatic pollutants (Duarte et al., 2014). After a manual check, a final list of gene sequences encoding enzymes potentially involved in degradation was prepared.

substrates possibly degraded by the communities are summarized. Common and microcosm-specific initial pollutants or intermediates for which presumptive degradation signatures were identified are specifically indicated in the Venn diagram. Solid lines represent single step reactions while dotted lines represent degradation steps where multiple reactions are involved (for details see Bargiela et al., 2015). Codes for proteins encoded by genes as follows: Abs, 4-aminobenzenesulfonate 3,4-dioxygenase; Bph, biphenyl dioxygenase; Bzn, benzene dioxygenase; Bzt, benzoate dioxygenase; Cat, catechol 2,3-dioxygenase; 2CB, 2-chlorobenzoate dioxygenase; Cum, *p*-cumate dioxygenase; Dhb, 2,3-Dihydroxybiphenyl dioxygenase; Dpp, 2,3-dihydroxyphenylpropionate dioxygenase; Gen, gentisate dioxygenase; Hna, 1-hydroxy-2-naphthoate dioxygenase; Hpc, homoprotocatechuate 2,3-dioxygenase; Ibu, ibuprofen-CoA dioxygenase; Ind, Rieske oxygenase involved in indole acetic acid degradation; Odm, 2-oxo-1,2-dihydroxyquinoline monooxygenase; Orc, orcinol hydroxylase; Pca, protocatechuate 3,4-dioxygenase; Pht, phthalate 4,5-dioxygenase; Thb, 2,2 ,3-trihydroxybiphenyl dioxygenase.

and phenanthrene) was found to be common for both microcosms. They include the transformation of biphenyl by Bph, 2,3-dihydroxybiphenyl by Dhb, benzoate by Bzt, indole-3-acetate by Ind, catechol by Cat, gentisate by Gen, 2-oxo-1,2-dihydroquinoline by Odm, 1-hydroxy-2-naphthoate by Hna, and 2-chlorobenzoate by 2-chlorobenzoate dioxygenase (2CB). Within them, genes encoding Cat were most abundant in both communities (UA: 17; AMM: 14), in agreement with the fact that catechol is the central intermediate for most cyclic aerobic hydrocarbons degradation (Pérez-Pantoja et al., 2009; Vilchez-Vargas et al., 2013). Gentisate and benzoate/2-chlorobenzoate may be most likely preferentially degraded by microorganisms in the AMM microcosm (10 Gen and 4 Bzt/2CB) in comparison to the UA microcosm (1 Gen and 1 Bzt/2CB). Genomic signatures for the degradation of orcinol (or 3,5-dihydroxytoluene) by Orc, phenylpropionate by Dpp, homoprotocatechuate by Hpc, and benzene by Bzn, were only found in the UA microcosm. The potential degradation of ibuprofen by Ibu, although not being a constituent of the crude oil but possibly originated from bilge water from the cruise lines or urban run-off, was also identified in UA microcosm. In stark contrast, the degradation of 4-aminobenzene-sulfonate by Abs, *p*-cumate by Cum, dibenzofuran by Thb, phthalate by Pht and protocatechuate by Pca, was characteristic for the AMM microcosm.

Note that within all pollutants predicted as being potentially degraded by bacteria inhibiting Ancona port (**Figure 2**), independently whether they are enriched with AMM or UA, only the potential degradation of Bargiela et al. Catabolic Alterations in Enrichment Microcosms

ibuprofen and 4-aminobenzene-sulfonate was not found associated to bacteria from other chronically crude oilcontaminated sites in oil-polluted sites along the coastlines of the Mediterranean Sea (Bargiela et al., 2015). This suggests that the pollution type and pollutant diversity in Ancona port, which receives chemicals such as alkyl benzene sulfonate detergents and drugs coming from human activities (Martínez-Pascual et al., 2010; Paíga et al., 2013), may have supported the presence of ibuprofen- and sulfonate benzene-growing bacteria. Such bacteria may be further stimulated by either the addition of UA or AMM, respectively.

### Experimental Analysis of Catabolic Capacities in AMM and UA Microcosms

Experimental validation assays were conducted to prove the extent of agreement with metagenomic-based predictions. For that, AMM and UA enrichment cultures were set up in triplicates as described in Section "Experimental Validations of Predicted Biodegradation Capacities," in which instead of Arabian light crude oil as the carbon source (used for the initial microcosms), naphthalene, 2,3-dihydroxybiphenyl, benzene, *p*-cumate, orcinol, 2-chlorobenzoate, phthalate and phenylpropionate (2 ppm each) were used. The capacity to degrade other pollutants predicted as potential substrates such as ibuprofen, phenanthrene, dibenzofuran, indole-3-acetic acid, 4-aminobenzene-sulfonate and quinoline, could not be experimentally proved because no analytical procedures could be designed for their analysis in the pollutant mix.

Samples were taken at 21 days of incubation at 20◦C. Fingerprinting by LC-ESI-QTOF-MS and GC-MS was used to confirm the degradation of the initial substrates as well as the existence of degradation intermediates in both cultures. A careful inspection of the mass signatures confirmed the lowering in the abundance level of naphthalene, 2,3-dihydroxybiphenyl, and 2-chlorobenzoate, and the presence of catechol, salicylate, gentisate, and benzoate in both microcosms (**Figure 3**). This demonstrates that the naphthalene-to-salicylate-togentisate, 2,3-dihydroxybiphenyl-to-benzoate-to-catechol, and 2-chlorobenzoate-to-catechol degradation pathways occurred or were active in both microcosms. Note that the lower abundance level of gentisate in AMM microcosm may correlate with the 10-fold overabundance of genes encoding Gen enzymes in AMM as compared to UA; this may decrease the pool of gentisate in the microcosm when growing in naphthalene. We further found a decreased level of *p*-cumate only associated to the AMM enrichment. Phthalate degradation mostly associated to the AMM microcosm, as confirmed by the higher extend of phthalate degradation by meaning of its residual percentage at the end of the assay (21 ± 1.8% in AMM vs. 92.2 ± 2.3% in UA) and the 22.2-fold higher abundance of protocatechuate in AMM compared to UA assays. In addition, decreased level of orcinol, benzene and phenylpropionate associated only to UA enrichments (**Figure 3**). Accordingly, the benzeneto-catechol, orcinol-, and phenylpropionate-degradation pathways occurred or were active in the UA microcosm,

while *p*-cumate degradation mostly occurred in the AMM enrichments.

The identification of degrading capacities on microcosms depends heavily on enrichment conditions (including cultivation time frame) and bacteria and protein abundance. While these drawbacks are known, the experimental data presented above (**Figure 3**) fully confirmed our sequence-based predictions (**Figure 2**) for the degradation of all eight pollutants tested in each of the two amendments. This suggests that the differences herein predicted in UA and AMM microcosms (**Figure 2**) are due to real biological differences and not random. Uncertainty remains only for phthalate degradation in UA microcosm: experimental analysis demonstrated the slight degradation of this chemical (**Figure 3**), which was not predicted by sequence analysis (**Figure 2**).

### Phylogenetic Identities of Catabolic Genes

We further attempted to analyze the contributions of particular sets of microbes to the entire reconstructed catabolic network, where multiple proteins from multiple organisms may contribute to organic pollutants' decomposition.

As the community structure of the two enrichment cultures was well-characterized (Gertler et al., 2015), the taxonomic affiliations of the catabolic genes identified could be unambiguously established at the family and phylum level. For that, we used tools recently published that provide a high level of confidence (Guazzaroni et al., 2013; Bargiela et al., 2015). **Figure 4** shows the contribution of members assigned to the different bacterial families and phyla in both microcosms to pollutant degradation. They included populations closely related to members of *Aeromonadaceae, Alcanivoracaceae, Alteromonadaceae, Halomonadaceae, Oceanospirillaceae, Piscirickettsiaceae, Pseudomonadaceae, Rhodobacteraceae, Vibrionaceae,* and *Xhantomonadaceae*, as well as to a lesser extent for the phyla Actinobacteria and Firmicutes. These comprise bacterial groups well known for their oil biodegrading capabilities (Yakimov et al., 2007; Jin et al., 2012; Guazzaroni et al., 2013). A further careful examination of the data presented

names/codes are identical to those presented in Figure 2.

in **Figure 4** clearly leads to the occurrence of a different pathway organization at organism level for the catabolism of 18 different pollutants predicted to be degraded.

As can be seen in **Figure 4**, members of *Alcanivoracaceae, Alteromonadaceae,* and *Rhodobacteraceae* were the major contributors to the networks. They contribute, in combination, to the degradation of 16 out of 18 pollutants predicted in the catabolic network, including dibenzofuran, phenanthrene, indolacetic acid, biphenyl, *p*-cumate, 2-chlorobenzoate, phenylpropionate, aminobenzenesulfonate and gentisate. This is in agreement with the fact that they were among the most abundant members in the established microcosms based on 16S rRNA (**Table 1** and **Figure 5**). Interestingly, *Pseudomonadaceae* which was the second most abundant microbial clade at the level of 16S rRNA in both microcosms (**Table 1** and **Figure 5**), did not contribute to the degradation network in AMM but it does in the UA microcosm (**Figure 4**), where it supports the biphenyl-to-benzoate and homoprotocatechuate degradation.

As shown in **Figure 4**, among the common degradation capacities, a number of observations can be made. First, the degradation of indole acetate by Ind was supported by members of *Alteromonadaceae* in AMM and *Rhodobacteraceae* in UA, which suggests a catabolic replacement. This was also observed for the degradation of biphenyl and benzene (by Bph/Bzn), most likely supported by members of the *Pseudomonadaceae* in UA but *Rhodobacteraceae* in AMM. We identified members of five proteobacterial families (*Aeromonadaceae, Alcanivoracaceae, Alteromonadaceae, Oceanospirillaceae*, and *Rhodobacteraceae*) and of the Firmicutes phylum as key groups for the degradation of gentisate (by Gen) in AMM. By contrast, only members of *Alteromonadaceae*were predicted to support gentisate catabolism in UA. In agreement with this it has been found that AMM promotes the growth of such multiple marine bacteria with the ability to utilize naphthalene (the precursor of gentisate) as a sole carbon in enrichment cultures (Hedlund et al., 1999). Also, the increased abundance of bacteria of the Firmicutes phylum has been demonstrated during bio-stimulation with ammonia (Guazzaroni et al., 2013). The naphthalene-to-gentisate catabolism by bacteria of the family *Alteromonadaceae* has also been found during microcosm assays using seawater and sediment samples obtained after an oil spill along the Korean shoreline without AMM addition (Jin et al., 2012); this agrees with the enrichment of gentisate catabolism by bacteria of this family in UA microcosm.

Multiple bacteria also contributed to the degradation of catechol (by Cat), with members of *Alcanivoracaceae, Alteromonadaceae,* and *Rhodobacteraceae* being common in both treatments. These bacterial groups are known for their capacity to degrade aromatics and haloaromatics to catechol, which can be further catabolised (Brusa et al., 2001; Antunes et al., 2011). Members of the Actinobacteria phyum and *Oceanospirillaceae* family contributed to catechol catabolism exclusively in the UA microcosm, whereas those of *Vibrionaceae* family did so in the AMM treatment. Note that, in accordance with the fact that *cat* genes are the most abundantly present (**Figure 2**) in both microcosms, the number of bacterial groups

involved in its catabolism was also the highest (8 in total; **Figure 4**). Therefore, a number of bacterial groups within the microcosms exhibit also partial catabolism redundancy.

Interestingly, we noticed that bacteria from the *Halomonadaceae* family contributed also to degradation of aromatics, particularly, 2-chlorobenzoate (through 2CB) and biphenyl (through Dhb) in the UA microcosm (**Figure 4**). This suggests that halomonads not only participate in the conversion of UA to AMM, which further stimulated growth of hydrocarbonoclastic bacteria (Gertler et al., 2015), but also play specific roles in degradation as herein suggested. This agrees with the fact that bacteria from the genus *Halomonas* are capable of degrading chlorobenzoates (de la Haba et al., 2011) and aromatics compounds such as benzoate and catechol (Piubeli et al., 2012), that are intermediate products of biphenyl and 2-chlorobenzoate degradation.

### CONCLUSION

Here, we report that different biostimulants applied in chronically polluted sediments have caused significant alteration in degradation capacities, while having no major effect on the taxonomic composition of microbial communities at the level of the family or higher. Experimental validation was conducted for at least eight of the predicted catabolic capacities, and good agreement with metagenomics-based predictions was observed. On the other hand, the metagenomics-guided metabolic reconstruction allowed us to refine the assignment of roles of community members in the utilization of multiple substrates and found different pathway organization at organism level. For example, while biphenyl degradation by Bph, DhB, and Bzt enzymes seems to be carried out by bacteria of *Pseudomonadaceae*, *Halomonadaceae,* and *Rhodobacteraceae* in UA, those of *Alteromonadaceae, Oceanospirillaceae, Picirickettsiaecae,* and Firmicutes may be involved in an alternative pathway in AMM. This demonstrates that different microbial members within microcosms obtained with different nitrogen sources may exhibit partial functional redundancy, and thus, may have a high level of common catabolic capacities. The present investigation provides an estimation of such common and distinct degrading capacities. Indeed, herein we found that 50% of the predicted degradation capacities were common for microorganisms in AMM and UA microcosms (**Figure 2**). However, according to the microbial biodegradation networks herein reconstructed, we also found that the two different biostimulants investigated, UA and AMM, have also changed substrate utilization capacities and preferences, which must be considered for the design of petroleum bioremediation techniques. This was demonstrated by showing that UA enriched for bacteria with the capability of degrading pollutants otherwise not degraded, or possibly degraded at low level, by those stimulated by the addition of AMM, and vice versa.

Therefore, the results of this study show that smart formulations based on the application of multiple nitrogen sources, rather than commonly used single sources (mostly AMM), for example, may increase the efficiency of the biological removal of the widest diversity of aromatic pollutants and could be essential to support effective biodegradation strategies in response to an oil spill incident or in response to chronical pollution. Thus, as herein demonstrated, the utilization of both AMM and UA in conjunction will have a double aim. In one side, AMM may most likely enhance the bio-stimulation of bacterial populations capable of degrading heavy oil components such as naphthalene, phenathrene and dibenzofuran, as well as sulfonated-benzenes and substituted benzoate derivatives such as p-cumate (**Figure 2**). In other side, UA will promote the growth of bacteria most active against benzene, orcinol-, ibuprofen- and phenyl-propionate (**Figure 2**). This will provoke a significant increase in multiple aromatics consumption in polluted areas. Having said that, this work seems to introduce a promising way for future oil-based contamination handling techniques. In this context, it would be very interesting to test the overall cleaning capacity (if any) on a real oil-contaminated marine sample. For that, also another point will be to use the combination of the UA and AMM, which was herein not presented in microcosm assays. It would be interesting to see their combinatory effect in the overall degradation capacity and taxonomic distribution of the microbial niche depending also on their ratio, so to find optimal nitrogen-containing formulations in real scenarios.

We would like to stress the attention to the fact that similarities regarding microbial community composition in the AMM microcosm from Ancona port with those reported in enrichments from surface water bodies at other Mediterranean sites, were found (Gertler et al., 2012). However, a similar comparison with the results from UA microcosm cannot be established because the limited information available. In fact, to the best of our knowledge, there have been only three studies that thoroughly investigated the use of UA in bioremediation trials (Koren et al., 2003; Knezevich et al., 2006; Nikolopoulou et al., 2013). Those studies, however, did not use UA in comparison to other nitrogen sources such as AMM, both in respect to their effect in microcosm population structures and catabolic preferences. Accordingly, herein we reported first evidences linking UA to catabolic preferences at the bacterial level, in comparison to the commonly use nitrogen source AMM.

### ACKNOWLEDGMENTS

This research was supported by the European Community Projects MAGICPAH (FP7-KBBE-2009-245226), ULIXES (FP7- KBBE-2010-266473) and KILL-SPILL (FP7-KBBE-2012-312139). We thank EU Horizon 2020 Program for the support of the Project INMARE H2020-BG-2014-2634486. This work was further funded by grants BIO2011-25012, PCIN-2014-107 and BIO2014-54494-R from the Spanish Ministry of Economy and Competitiveness. The authors gratefully acknowledge the financial support provided by the European Regional Development Fund (ERDF). The present investigation was also funded by the Spanish Ministry of Economy and Competitiveness within the ERA NET IB2, grant number ERA-IB-14-030. FM was supported by Università degli Studi di Milano, European Social Fund (FSE) and Regione Lombardia (contract "Dote Ricerca").

## REFERENCES

Yakimov, M. M., Timmis, K. N., and Golyshin, P. N. (2007). Obligate oil-degrading marine bacteria. *Curr. Opin. Biotechnol.* 18, 257–266. doi: 10.1016/j.copbio.2007.04.006

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Bargiela, Gertler, Magagnini, Mapelli, Chen, Daffonchio, Golyshin and Ferrer. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Novel circular single-stranded DNA viruses identified in marine invertebrates reveal high sequence diversity and consistent predicted intrinsic disorder patterns within putative structural proteins

*Karyna Rosario, Ryan O. Schenck, Rachel C. Harbeitner, Stephanie N. Lawler and Mya Breitbart\**

*College of Marine Science, University of South Florida, St. Petersburg, FL, USA*

#### *Edited by:*

*Eamonn P. Culligan, University College Cork, Ireland*

#### *Reviewed by:*

*Kenneth Stedman, Portland State University, USA Purificacion Lopez-Garcia, Centre National de la Recherche Scientifique, France*

#### *\*Correspondence:*

*Mya Breitbart, College of Marine Science, University of South Florida, 140 7th Avenue South, St. Petersburg, FL 33701, USA<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Evolutionary and Genomic Microbiology, a section of the journal Frontiers in Microbiology*

*Received: 27 April 2015 Accepted: 23 June 2015 Published: 10 July 2015*

#### *Citation:*

*Rosario K, Schenck RO, Harbeitner RC, Lawler SN and Breitbart M (2015) Novel circular single-stranded DNA viruses identified in marine invertebrates reveal high sequence diversity and consistent predicted intrinsic disorder patterns within putative structural proteins. Front. Microbiol. 6:696. doi: 10.3389/fmicb.2015.00696* Viral metagenomics has recently revealed the ubiquitous and diverse nature of singlestranded DNA (ssDNA) viruses that encode a conserved replication initiator protein (Rep) in the marine environment. Although eukaryotic circular Rep-encoding ssDNA (CRESS-DNA) viruses were originally thought to only infect plants and vertebrates, recent studies have identified these viruses in a number of invertebrates. To further explore CRESS-DNA viruses in the marine environment, this study surveyed CRESS-DNA viruses in various marine invertebrate species. A total of 27 novel CRESS-DNA genomes, with Reps that share less than 60.1% identity with previously reported viruses, were recovered from 21 invertebrate species, mainly crustaceans. Phylogenetic analysis based on the Rep revealed a novel clade of CRESS-DNA viruses that included approximately one third of the marine invertebrate associated viruses identified here and whose members may represent a novel family. Investigation of putative capsid proteins (Cap) encoded within the eukaryotic CRESS-DNA viral genomes from this study and those in GenBank demonstrated conserved patterns of predicted intrinsically disordered regions (IDRs), which can be used to complement similarity-based searches to identify divergent structural proteins within novel genomes. Overall, this study expands our knowledge of CRESS-DNA viruses associated with invertebrates and explores a new tool to evaluate divergent structural proteins encoded by these viruses.

Keywords: single-stranded DNA virus, CRESS-DNA virus, circular DNA virus, intrinsically disordered proteins (IDPs), intrinsically disordered regions (IDRs), marine invertebrate, crustaceans

### Introduction

Viral metagenomics, or shotgun sequencing of total nucleic acids from purified virus particles, enables examination of viral communities without prior knowledge of the viruses present, thus resulting in an unprecedented view of viral diversity (Breitbart et al., 2002; Edwards and Rohwer, 2005; Angly et al., 2006). This technique has uncovered many novel viral types and extended the environmental distribution of known viral groups (Delwart, 2007; Rosario and Breitbart, 2011). In particular, the incorporation of rolling circle amplification (RCA) into viral metagenomic studies has unearthed a high diversity and wide distribution of eukaryotic viruses with circular, single-stranded DNA (ssDNA) genomes that encode a conserved replication initiator protein (Rep; Delwart and Li, 2012; Rosario et al., 2012a). Before the metagenomics era, eukaryotic circular Rep-encoding ssDNA (CRESS-DNA) viruses were only known in agricultural and medical fields since they are known plant (*Geminiviridae* and *Nanoviridae*) and vertebrate (*Circoviridae*) pathogens. However, over the past decade metagenomic approaches have revealed the ubiquitous nature of eukaryotic CRESS-DNA viruses, with reports from various environments, including deep-sea vents (Yoshida et al., 2013), Antarctic lakes and ponds (López-Bueno et al., 2009; Zawar-Reza et al., 2014), wastewater (Rosario et al., 2009b; Roux et al., 2013; Kraberger et al., 2015; Phan et al., 2015), freshwater lakes (Roux et al., 2012, 2013), oceans (Rosario et al., 2009a; Labonte and Suttle, 2013; Roux et al., 2013), hot springs (Diemer and Stedman, 2012), the near-surface atmosphere (Whon et al., 2012; Roux et al., 2013), and soils (Kim et al., 2008; Reavy et al., 2015).  Notably, CRESS-DNA viruses similar to circoviruses, which were previously thought to only infect vertebrates, have now been identified in a myriad of invertebrates, including insects (Ng et al., 2011; Rosario et al., 2011, 2012b; Dayaram et al., 2013; Padilla-Rodriguez et al., 2013; Pham et al., 2013a,b; Garigliany et al., 2015), crustaceans (Dunlap et al., 2013; Hewson et al., 2013a,b; Ng et al., 2013; Pham et al., 2014), cnidarians (Soffer et al., 2014), and gastropods (Dayaram et al., 2015a), suggesting that CRESS-DNA viruses may be prevalent amongst unexplored taxa.

Well-studied viruses from the *Circoviridae*, *Nanoviridae*, and *Geminiviridae* families demonstrate the rapid evolutionary potential of CRESS-DNA viruses due to high nucleotide substitution rates (Duffy et al., 2008; Duffy and Holmes, 2009) as well as mechanistic predispositions to recombination (Lefeuvre et al., 2009; Martin et al., 2011). These characteristics, combined with the high level of recently reported diversity, highlight the need to continually revisit taxonomic classification of this viral group to add new species, genera and/or families. However, this task is complicated by the fact that many of the CRESS-DNA virus genomes exhibit novel genome architectures, only share similarities to the highly conserved Rep of known viruses, and have similarities to viruses belonging to multiple different taxonomic groups (Rosario et al., 2012a; Roux et al., 2013). In addition, the definitive hosts for many of these CRESS-DNA viruses remain unknown, hindering their classification according to traditional standards.

CRESS-DNA viruses are characterized by small genomes (∼1.7–3 kb) that contain 2–6 protein-encoding genes. The smallest monopartite CRESS-DNA viruses, members of the *Circoviridae* family, exhibit only two major open reading frames (ORFs), which encode a Rep and a capsid protein (Cap). Many of the novel eukaryotic CRESS-DNA viral genomes obtained from environmental samples or individual organisms through either metagenomic sequencing or degenerate PCR (herein referred to as "metagenomic CRESS-DNA viruses") exhibit similarities to circoviruses and have been referred to as 'circo-like' viruses. Although many of the metagenomic circo-like virus genomes are highly divergent, these surveys have uncovered a novel CRESS-DNA viral group, the proposed Cyclovirus genus (Li et al., 2010a). Cycloviruses, which form a sister group to the *Circovirus* genus within the family *Circoviridae*, have been identified from both vertebrates (Li et al., 2010a; Smits et al., 2013; Tan Le et al., 2013; Garigliany et al., 2014; Zhang et al., 2014) and invertebrates (Rosario et al., 2011, 2012b; Dayaram et al., 2013, 2014, 2015b; Padilla-Rodriguez et al., 2013).

Similarities to circoviruses are mainly based on the Rep whereas the second major ORF in novel circo-like metagenomic CRESS-DNA viruses generally does not have any significant matches in the database but is assumed to encode for a structural protein based on the genomic architecture of known circoviruses. In lieu of significant matches to known structural proteins in the GenBank database, it is important to investigate putative novel Caps in CRESS-DNA viruses to provide evidence regarding theirstructural function. A potential avenue to identify conserved patterns in highly divergent structural proteins, such as those observed in novel metagenomic CRESS-DNA viruses, is to investigate the presence of predicted intrinsically disordered regions (IDRs). IDRs are regions within a protein that lack a rigid or fixed (i.e., ordered) structure, allowing a protein to exist in different states depending on the substrate with which it is interacting (Dunker et al., 2001; Brown et al., 2011). Research examining IDRs within viral proteomes has revealed that smaller viral genomes, such as those of CRESS-DNA viruses, contain a higher proportion of predicted disordered residues than larger viruses (Xue et al., 2012, 2014; Pushker et al., 2013). Therefore it has been suggested that small viruses may exploit IDRs to encode multifunctional proteins (Xue et al., 2012, 2014; Pushker et al., 2013). Since structural proteins in several viral families commonly contain IDRs (Chen et al., 2006; Goh et al., 2008a,b; Chang et al., 2009; Jensen et al., 2011), the presence of similar patterns of predicted disorder amongst unidentified CRESS-DNA proteins may provide one line of evidence for these proteins representing putative Caps.

To contribute to efforts exploring the diversity of CRESS-DNA viruses in invertebrates, this study investigated various marine invertebrate species for the presence of these viruses. A total of 27 novel CRESS-DNA genomes were recovered from 21 invertebrate species, expanding the known diversity of CRESS-DNA viruses associated with marine organisms and providing the first evidence of viruses associated with some under-sampled taxa. The well-conserved Rep of CRESS-DNA viruses was used to explore the relationships between these novel viruses and previously reported eukaryotic CRESS-DNA viruses in GenBank, including metagenomic CRESS-DNA viruses. In addition, the non-Rep-encoding ORFs (i.e., putative Caps) within these genomes were investigated for IDRs. Disorder prediction methods suggest that CRESS-DNA viral Caps exhibit conserved patterns of predicted disorder, which can be used to complement similarity-based searches to identify structural proteins within novel CRESS-DNA viral genomes.

## Materials and Methods

### Sample Processing and Genome Discovery

CRESS-DNA viruses were investigated in a variety of marine invertebrate species that were collected as samples of opportunity (**Table 1** and Supplementary Table S1). Specimens were identified with the highest degree of taxonomic resolution possible based on morphology. Whole organisms or tissue sections were serially rinsed three times using sterile SM Buffer [0.1 M NaCl, 50 mM Tris-HCl (pH 7.5), 10 mM MgSO4]. Viral particles were partially purified from each specimen prior to DNA extraction. For this purpose, samples were homogenized in one of two ways depending on the size of the specimen. Smaller organisms or dissected tissues that could be placed in a 1.5 ml microcentrifuge tube were homogenized in 1 ml of sterile SM Buffer through bead-beating using 1.0 mm sterile glass beads in a bead beater (Biospec Products). Homogenates were then centrifuged at 6000 × *g* for 6 min. Larger organisms or tissues of dissected organisms, such as muscle or gonads, were placed in a gentleMACSTM M tube (Miltenyl Biotec) containing 3 ml of sterile SM buffer. Samples were then homogenized using a gentleMACS dissociator (Miltenyl Biotec) followed by centrifugation at 6000 × *g* for 9 min. The supernatant from both homogenization methods was filtered through a 0.45 µm Sterivex filter (Millipore) and nucleic acids were extracted from 200 µl of filtrate using the QIAmp MinElute Virus Spin Kit (Qiagen).

DNA extracts were amplified through RCA using the illustra TempliPhi Amplification kit (GE Healthcare) to enrich for small circular templates (Kim et al., 2008; Kim and Bae, 2011). RCA-amplified DNA was digested with a suite of FastDigest restriction enzymes (Life Technologies; BamHI, EcoRV, PdmI, HindIII, KpnI, PstI, XhoI, SmaI, BgiII, EcoRI, XbaI, and NcoI) following manufacturer's instructions in separate reactions to obtain complete, unit-length genomes for downstream cloning and sequencing. Restriction enzyme digested products were resolved on an agarose gel and bands ranging in size from 1000 to 4000 bp were excised and cleaned using the Zymoclean Gel DNA Recovery Kit (Zymo Research). Products resulting from blunt-cutting enzyme digestions were cloned using the CloneJET PCR Cloning kit (Life Technologies), whereas products containing sticky ends were cloned using pGEM-3Zf(+) vectors (Promega) pre-digested with the appropriate enzyme. All clones were commercially Sanger sequenced using vector primers and genomes exhibiting significant similarities to eukaryotic CRESS-DNA viruses were completed through primer walking.

### Genome Annotation

Genomes were assembled using Sequencher 4.1.4 (Gene Codes Corporation). Putative ORFs *>*100 amino acids were identified and annotated using SeqBuilder version 11.2.1 (Lasergene). Partial genes or genes that seemed interrupted were analyzed for potential introns using GENSCAN (Burge and Karlin, 1997). The potential origin of replication (*ori*) for each genome was identified by locating a canonical nonanucleotide motif (NANTATTAC; Rosario et al., 2012a) and confirming predicted stem-loop structures using Mfold with constraints applied to prevent hairpin formation within the nonanucleotide motif and a folding temperature set at 17◦C (Zuker, 2003). Final annotated genomes have been deposited to GenBank with accession numbers KR528543–KR528569.

### Database Sequences and Sequence Analysis

To conduct sequence comparisons, members of the *Circovirus* genus, as well as complete eukaryotic CRESS-DNA viral genomes obtained from environmental samples or individual organisms through either metagenomic sequencing or degenerate PCR (herein referred to as "metagenomic CRESS-DNA viruses") were retrieved from GenBank. Since the Rep is the only conserved protein among CRESS-DNA viruses (Ilyina and Koonin, 1992; Rosario et al., 2012a) this protein was used to compare the different genomes. Rep pairwise identities were calculated using SDT v1.2 (Muhire et al., 2014) and summarized using heat maps generated in R (R Core Team, 2014). A maximum likelihood (ML) phylogenetic tree based on Rep amino acid sequences was also constructed. For this purpose, alignments were performed in MEGA 6.06 (Tamura et al., 2013) using the MUSCLE algorithm (Edgar, 2004) and manually edited. Sequences were inspected for the presence of conserved amino acid motifs that have been shown to play a role in rolling circle replication (RCR) of eukaryotic CRESS-DNA viruses, including three RCR and three superfamily 3 (SF3) helicase motifs (Gorbalenya et al., 1990; Ilyina and Koonin, 1992; Gorbalenya and Koonin, 1993; Rosario et al., 2012a). Although all the recently reported CRESS-DNA viruses are included in the heatmap, only sequences exhibiting all six motifs are included in the phylogenetic analysis. In addition, divergent regions that were poorly aligned, as shown by a high percentage of gaps, were removed from the alignment (Supplementary Data Sheet 1). Since the *Nanoviridae* and *Geminiviridae* are also CRESS-DNA viral families that are evolutionarily related to the *Circoviridae* (Ilyina and Koonin, 1992; Rosario et al., 2012a), select representatives of these families were included in the phylogenetic analysis. The ML phylogenetic tree was inferred using PHYML (Guindon et al., 2010) implementing the best substitution model (rtRev+I+G+F; Dimmic et al., 2002) according to ProtTest (Abascal et al., 2005). Branch support was assessed using the approximate likelihood ratio test (aLRT) SH-like method (Anisimova and Gascuel, 2006).

### Intrinsically Disordered Region (IDR) Analysis of Putative Capsid Proteins

To determine if the non-Rep-encoding ORFs from the CRESS-DNA viral genomes presented here (*n* = 25), circoviruses (*n* = 15), and metagenomic CRESS-DNA viruses (*n* = 259; including 37 cycloviruses) represent putative Caps, these proteins were evaluated for IDRs. Disordered protein regions were predicted using the DisProt VL3 disorder predictor (Obradovic et al., 2003; Sickmeier et al., 2007). This artificial neural network utilizes an ensemble of feed forward neural networks with 20 attributes (18 amino acid frequencies, average flexibility,

TABLE 1 | CRESS-DNA genomes identified in this study, the organism they were obtained from, and genome details (acronym, genome length, nonanucleotide motif, genome type, and ORFs identified).

1*Genome names contain abbreviation aCV for associated circular virus or aCG for associated circular genome. ID within parentheses corresponds to ID used throughout the paper.*

2*Nonanucleotide motif sequences that were not identified within a stem-loop structure are denoted with an asterisk (*∗*).*

3*Non-Rep encoding ORFs were identified as putative capsid proteins based on BLAST results. However, many non-Rep-encoding ORFs did not exhibit any significant matches (marked with an asterisk*∗*).*

and sequence complexity; Obradovic et al., 2003). Disorder disposition scores above a 0.5 threshold indicate intrinsic disorder. Counts and statistical analysis for the fraction of disorder- and order-promoting amino acid residues was conducted using R with the "seqinr" package (Charif and Lobry, 2007).

### Results

A total of 27 CRESS-DNA genomes were recovered from 21 marine invertebrates (**Table 1**). Most of the recovered genomes (66.7%) were identified from *Crustacea*, mainly from the order *Decapoda*. Recovered genomes ranged in size from 1063 to 2469 nt and exhibited a variety of genome architectures. Of the 27 genomes identified, 23 exhibited a common putative *ori* marked by a conserved nonanucleotide motif (NANTATTAC) at the apex of a predicted stem-loop structure (**Table 1**). The remaining four genomes lacked a stem-loop structure (*n* = 2) or a stem-loop structure and a nonanucleotide motif (*n* = 2). Genomes lacking the canonical nonanucleotide motif could not be assigned to any genome type; therefore only 25 genomes were assigned to genomic architecture types previously described by Rosario et al. (2012a) (**Figure 1**). The predominant genomic architecture observed was Type I (*n* = 13), which is typical of members of the *Circovirus* genus. However, other genomic architectures were observed including Types II (*n* = 5), III (*n* = 1), IV (*n* = 1), V (*n* = 3), and VII (*n* = 2) (**Figure 1**). It is important to note that genomes exhibiting a Type VII genome architecture only exhibit a single major ORF encoding a Rep. This type of architecture is observed in genomic components of multipartite viruses from the *Nanoviridae* family and satellite DNA molecules that require helper viruses for encapsidation (Gronenborn, 2004; Briddon and Stanley, 2006). Therefore genomes exhibiting only a single major ORF may

represent partial genomes of multipartite viruses or non-viral mobile genetic elements such as plasmids (Rosario et al., 2012a).

The majority of the CRESS-DNA viruses detected in marine invertebrates were most similar to viral sequences identified through metagenomic surveys of marine samples (Supplementary Table S1). However, one of genomes, *Lytechinus variegatus* variable sea urchin associated circular virus\_I0021, was most similar to plant viruses from the *Geminiviridae* family. Most of the viral genomes had database similarities for the Rep; except for *Sicyonia brevirostris* brown rock shrimp associated circular virus\_I0722, which only had similarities for the putative Cap (Supplementary Table S1).

Pairwise identities indicate that the CRESS-DNA viruses detected in marine invertebrates share less than 60.1% sequence identity (average sequence identity = 26.04%) with previously identified Reps from CRESS-DNA viruses in GenBank, indicating that these viruses represent novel species (**Figure 2**). Twenty-one of the 27 recovered Reps contained all six conserved RCR and helicase motifs (see Materials and Methods) and were used for phylogenetic analysis. Analysis of these Reps with representative

CRESS-DNA viral Reps from GenBank, including available metagenomic CRESS-DNA viral Reps, show that most of the sequences from marine invertebrate associated viruses detected here are more closely related to circo-like viruses recovered through metagenomic surveys of the marine environment than

to previously defined CRESS-DNA viral groups (**Figure 3**). Eleven of the 21 Reps from marine invertebrate associated viruses do not form distinct clusters with each other or any known sequences (**Figure 3**). However, ten of the Reps form a well-supported clade that also includes sequences detected

FIGURE 3 | Multifurcation maximum likelihood phylogenetic reconstruction based on the Reps of CRESS-DNA genomes recovered here, metagenomic CRESS-DNA viruses, cycloviruses, circoviruses, and representative members of the *Nanoviridae and Geminiviridae* families. Reps obtained from CRESS-DNA genomes obtained in this study are highlighted in blue font. Branches are colored for the different CRESS-DNA viral groups including the Marine

Clade 1 (red), circoviruses (purple), cycloviruses (pink), nanoviruses (orange), and geminiviruses (green). Representative nanoviruses (*n* = 4) and geminiviruses (*n* = 15) have been condensed into their family names. Reps from genomes exhibiting a single ORF are highlighted using an asterisk (∗). Branches with less than 60% aLRT branch support have been collapsed. Description of acronyms used can be found in Supplementary Table S4.

in the Gulf of Mexico (GOM00443; JX904231.1), Straight of Georgia (JX904106.1), McMurdo Ice Shelf (YP\_009047125.1; YP\_009047137.1), and a semi-enclosed shallow estuary (Avon-Heathcote Estuary associated circular virus 24; AJP36460.1). Pairwise identity scores indicate that all members of this clade, named Marine Clade 1 for the purposes of this study, share more than 32.7% identity, with an average pairwise identity score of 47.2% (**Figure 2**). Members of the Marine Clade 1 seem to be more closely related to members of the *Nanoviridae* (31.95% average pairwise identity) than any other known CRESS-DNA viral group; however, members of this clade exhibit different genomic architectures compared to these plant viruses. CRESS-DNA viral genomes from the Marine Clade 1 encode two major ORFs in an ambisense organization (i.e., Type I architecture), which is similar to members of the *Circoviridae*, rather than the single ORF, Type VII genome organization observed in genomic components from the *Nanoviridae*.

### Capsid Analysis

Only half of the CRESS-DNA viral genomes described here contained an ORF that had significant BLASTX matches (e-value *<* 0.001; amino acid identities ranging from 26–54%) to proteins annotated as putative Caps in GenBank (**Table 1**). Furthermore, most of the matches in the database were to putative CRESS-DNA viral Caps detected through metagenomic surveys, which are not supported by biochemical data and have not necessarily been well curated. Therefore, alternative methods were explored to investigate non-Rep-encoding ORFs (i.e., putative Caps) found in CRESS-DNA viral genomes.

The majority of metagenomic CRESS-DNA viruses reported from marine invertebrates in this study and in GenBank are most similar to previously described circoviruses. Therefore, the predicted IDP profiles of well-characterized members of the *Circovirus* genus were examined in an effort to identify conserved patterns in structural proteins encoded by these viruses. These circovirus IDP profiles were then compared against profiles observed in cycloviruses (the proposed sister group to the circoviruses, which exhibit conserved features and share high identities with circoviruses) and other metagenomic CRESS-DNA viruses.

The DisProt VL3 disorder prediction analysis revealed that Caps encoded by members of the *Circovirus* genus (*n* = 15) exhibit one of two protein disorder profiles, distinguished here as Type A or Type B, based on the first 125 amino acids of these proteins (**Figure 4A**). Type A Caps exhibit IDP profiles that are predicted to have the highest degree of disorder closest to the N-terminus (i.e., amino acid residues 1–50) before the profile tapers to a structured region with variable predicted disorder. Type A Caps exhibit significant enrichment for amino acid residues that promote disorder (R, K, E, P, S, Q, and A) within the first 50 residues relative to amino acid residues 51– 125 (ANOVA with *post hoc* Tukey's HSD; *p <* 0.05) and a depletion of order promoting amino acid residues (W, C, F, I, Y, V, L, and N) within the first 25 residues relative to amino acid residues 26–125 (ANOVA with *post hoc* Tukey's HSD; *p <* 0.05; **Figure 4B**). On the other hand, Type B Caps exhibit IDP profiles that peak in predicted disorder between amino acid residues 26–75. Type B Caps show an enrichment of disorder promoting residues between residue positions 26 through 75, whereas there is a depletion of predicted order promoting residues in this region compared to residues 1–25 and 76–125 (**Figure 4B**). Beyond 125 amino acids, IDP profiles exhibited more structured regions for both Types A and B Caps, with no distinguishable predicted disorder pattern (**Figure 4A**).

The overwhelming majority of Caps from the *Circovirus* genus (86.7%) exhibited Type A IDP profiles; however, two avian circoviruses, Finch circovirus (YP\_803551.1) and Beak and feather disease virus (NP\_047277.1), had Type B IDP profiles (**Table 2** and Supplementary Table S5). Similarly, 97.3% of cyclovirus putative Caps (*n* = 37) exhibited Type A IDP profiles. Comparison of IDP profiles showed that a majority of metagenomic CRESS-DNA viruses also contained patterns of increased predicted disorder at the N-terminus of the putative Cap, consistent with the *Circoviridae*. Interestingly, Type B IDP profiles were more prevalent among putative Caps from metagenomic CRESS-DNA viral genomes in GenBank (10.8%; *n* = 222) and the novel genomes reported in this study (56%; *n* = 25). Notably, 7 of the 10 viruses found in the Marine Clade 1 described here exhibit Type B Caps. Among the total 299 CRESS-DNA genome sequences analyzed, most putative Caps exhibit Type A IDP profiles (69.9%), followed by Type B (13%). Notably, most of the putative Caps lacking a significant match in the database exhibited one these profiles.

### Discussion

Metagenomic studies have revealed a prodigious amount of diversity in eukaryotic CRESS-DNA viruses in the marine environment (Rosario et al., 2009a; Rosario and Breitbart, 2011; Labonte and Suttle, 2013; McDaniel et al., 2014). However, few studies have isolated these viruses directly from organisms. Building upon recent studies suggesting that CRESS-DNA viruses are associated with marine invertebrates (Dunlap et al., 2013; Hewson et al., 2013a,b; Ng et al., 2013; Pham et al., 2014; Soffer et al., 2014; Dayaram et al., 2015a), this study investigated a variety of marine invertebrates, including under sampled taxa, for the presence of these viruses. Viral genomes presented here were primarily recovered from *Crustacea,* suggesting that this subphylum harbors a rich diversity of CRESS-DNA viruses. This is consistent with previous research that identified CRESS-DNA viruses in copepods (Dunlap et al., 2013), which are the most abundant members of mesozooplankton (Kleppel et al., 1996), as well as different species of shrimp (Ng et al., 2013; Pham et al., 2014), which comprise some of the world's most important food sources (Goss et al., 2000; Paezosuna, 2003). In addition, this is the first study to report viruses associated with marine snails, anemones, sea squirts, and several crab species. Although a definitive host for these viruses cannot be assigned with the present data, this study reveals the need for further examination of viruses associated with common marine invertebrates and experiments to determine their potential impact, if any, on the ecology of these organisms. The grouping of the invertebrateassociated CRESS-DNA viruses reported here with metagenomic

FIGURE 4 | (A) Representative IDP prediction profiles for Type A and Type B capsid proteins (Caps) from the Disprot VL3 predictor. Type A and Type B IDP prediction profiles are based on the Porcine circovirus 2 Cap (NP\_937957.1) and the Beak and feather disease virus Cap (NP\_047277.1), respectively. The grey shaded area represents the amino acid residue interval used in (B). (B) Graphs showing the fraction of predicted disordered (red bars) and ordered (blue bars) residues within discrete amino acid intervals for Type A and Type B Caps identified from all CRESS-DNA viral genomes

analyzed in this study. Significantly different amino acid intervals for each Cap type are distinguished using letters ("A", "B", "C", "D" for statistics based on percentage of predicted disordered residues) or numbers ("1", "2", "3", "4" for statistics based on percentage of predicted ordered residues; ANOVA with *post hoc* Tukey's HSD; *p <* 0.05). Note that the percentage of predicted disordered and ordered residues does not add to 100% due to the presence of residues that are not considered either disordered or ordered (i.e., H, M, T, and D).

CRESS-DNA viruses implies that marine invertebrates may serve as hosts for many of the sequences obtained from marine environments.

The marine invertebrate associated CRESS-DNA viruses identified here are only distantly related to known members of the *Circoviridae* and may represent novel groups. Approximately one third of the novel sequences reported here belong to the Marine Clade 1, whose members share an average pairwise identity of 47.2%. Members of this viral clade share an average pairwise