**VIRUS ECOLOGY AND DISTURBANCES: IMPACT OF ENVIRONMENTAL DISRUPTION ON THE VIRUSES OF MICROORGANISMS**

**Topic Editors Heather K. Allen and Stephen T. Abedon**

#### *FRONTIERS COPYRIGHT STATEMENT*

© Copyright 2007-2015 Frontiers Media SA. All rights reserved.

All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.

The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.

Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.

Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.

As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.

All copyright, and all rights therein, are protected by national and international copyright laws.

The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use.

**ISSN** 1664-8714 **ISBN** 978-2-88919-448-3 **DOI** 10.3389/978-2-88919-448-3

## *ABOUT FRONTIERS*

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

## *FRONTIERS JOURNAL SERIES*

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing.

All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

## *DEDICATION TO QUALITY*

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view.

By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

## *WHAT ARE FRONTIERS RESEARCH TOPICS?*

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area!

Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

## **VIRUS ECOLOGY AND DISTURBANCES: IMPACT OF ENVIRONMENTAL DISRUPTION ON THE VIRUSES OF MICROORGANISMS**

Topic Editors: **Heather K. Allen,** National Animal Disease Center,USA

**Stephen T. Abedon,** The Ohio State University, USA

A no-doubt microbe-laden formation found within Excelsior Geyser Crater in Yellowstone National Park. The scale of the formation is roughly five meters in diameter. Photograph by S.T. Abedon.

Viruses infect numerous microorganisms including, predominantly, Bacteria (bacteriophages or phages) but also Archaea, Protists, and Fungi. They are the most abundant and ubiquitous biological entities on Earth and are important drivers of ecosystem functioning. Little is known, however, about the vast majority of these viruses of microorganisms, or VoMs. Modern techniques such as metagenomics have enabled the discovery and description of more presumptive VoMs than ever before, but also have exposed gaps in our understanding of VoM ecology. Exploring the ecology of these viruses – which is how they interact with host organisms, the abiotic environment, larger organisms, and even other viruses across a variety

of environments and conditions – is the next frontier. Integration of a growing molecular understanding of VoMs with ecological studies will expand our knowledge of ecosystem dynamics.

Ecology can be studied at multiple levels including individual organisms, populations, communities, whole ecosystems, and the entire biosphere. Ecology additionally can consider normal, equilibrium conditions or instead perturbations. Perturbations are of particular interest because measuring the effect of disturbances on VoM-associated communities provides important windows into how VoMs contribute to ecosystem dynamics. These disturbances in turn can be studied through in vitro, in vivo, and in situ experimentation, measuring responses by VoM-associated communities to changes in nutrient availability, stress, physical disruption, seasonality, etc., and could apply to studies at all ecological levels. These are considered here across diverse systems and environments.

# Table of Contents


## Virus ecology and disturbances: impact of environmental disruption on the viruses of microorganisms

#### *Heather K. Allen1 \* and Stephen T. Abedon2 \**

*<sup>1</sup> National Animal Disease Center, Agricultural Research Service, United States Department of Agriculture, Ames, IA, USA*

*<sup>2</sup> Department of Microbiology, The Ohio State University, Mansfield, OH, USA*

*\*Correspondence: heather.allen@ars.usda.gov; abedon.1@osu.edu*

#### *Edited by:*

*John R. Battista, Louisiana State University, USA*

**Keywords: bacteriophages, environmental disturbance, phage ecology, aquatic microbiology, phage therapy, metaviromes, evolution, microarrays**

Viruses of microorganisms (VoMs)—the encapsidated, acelluluar parasites of archaea, bacteria, and microbial eukaryotes (Hyman and Abedon, 2012)—can be both responsive to and the source of change in their environments. Defining the diverse VoMs in an environment has become easier with the advent of various new technologies, so analyzing VoM populations or communities involved with disturbances could be the next step to drive research forward. Contributions assembled in this Research Topic explore the effects of disturbances on one type of VoM, the bacteriophages (phages), as well as disturbances caused by phages. The collection presented in this ebook opens with an editorial describing the importance of change to life on Earth, and by defining the multiple levels at which change can occur: genetic, organismal, and environmental (Allen and Abedon, 2013). We have arranged and introduced the manuscripts in this collection based on which of these levels they address.

## **GENETIC CHANGE**

Organismal responses to certain disturbances can be tightly linked to evolution. Two papers in this Research Topic contribute to evolutionary theory using two entirely different approaches: RNA virus populations *in vitro* and viral metagenomes *in silico*. Goldhill et al. studied the effect of repeated heat stress on an RNA virus to test the relationship between the maintenance of phenotype despite ongoing genetic mutations, which is called robustness, and the heritable ability to adapt to environmental disturbances, or evolvability (Goldhill et al., 2014). Theirs is one of the first studies to investigate these questions in the laboratory, using the φ6 RNA phage of *Pseudomonas syringae* pv. *phaseolicola*. The contribution by Mizuno et al. examines marine viral metagenomes for metaviromic islands, which are gaps in the metagenome represented by significantly fewer sequences than neighboring regions of high coverage (Mizuno et al., 2014). They show that the metaviromic islands frequently correspond to genes involved in genome packaging, molecular recognition, and host recognition, the latter of which is suggestive of providing high diversity within lineages to maximize predatory opportunities.

## **ORGANISMAL CHANGE**

Next are two papers that use examples from agriculture to consider the relationship between phages and disturbances at the organismal level. Bearson et al. show that subinhibitory concentrations of the agricultural antibiotic carbadox induce prophages *in vitro* in the foodborne pathogen *Salmonella* Typhimurium (Bearson et al., 2014). Additional data demonstrate phage-mediated gene transfer by phages induced by this antibiotic. The piece by Meaden and Koskella is a review on phage therapy, with an emphasis on agricultural examples (Meaden and Koskella, 2013). Phage therapy is a potential alternative to antibiotic therapy, and the authors thoughtfully explore the potential risks and unknowns associated with phage therapy.

## **ENVIRONMENTAL CHANGE**

The effect of environmental change on viruses is explored here in aquatic ecosystems. Storms are a powerful macro-scale disturbance to any system, and Williamson et al. examines the effect of stormwater runoff and the viruses therein on freshwater microbial communities (Williamson et al., 2014). Unique to this study is consideration of both the planktonic and particle-associated viral communities. Moving to a marine system, Labonté and Suttle undertake a deep analysis of the *Gokushovirinae* family of single-stranded DNA viruses to assess their diversity on a nearly global scale (Labonté and Suttle, 2013). The lineages studied showed environment-dependent biogeographic distributions.

A special type of environmental change can occur in hostassociated ecosystems, where interaction with hosts adds an additional layer of complexity. The work by Payet et al. continues the aquatic theme with an analysis of viral diversity in coral reef aquatic ecosystems (Payet et al., 2014). They undertook an extended spatiotemporal analysis of viral abundance and lytic activity in the South Pacific Ocean, with the results revealing viral dynamism linked to microbial communities and environmental factors. Next, Yamada presents a hypothesis on the regulation of virulence by filamentous phages of *Ralstonia*, which is the causative agent of bacterial wilt disease in certain agricultural plants (Yamada, 2013). The author suggests a phage-encoded open reading frame that is responsible for mediating virulence. Finally, Santos et al. review the potential for microarrays to dissect viral-microbe interactions, which could be applied to samples from any environment (Santos et al., 2014). The use of various probes, including glycans, peptides, and nucleic acids, provides a powerful platform for investigating a wide range of questions in viral ecology.

## **CONCLUDING REMARKS**

The articles in this Research Topic uniquely and individually address one of the tiniest biological entities on Earth, the viruses of microorganisms, and the enormous role they play in both causing and effecting environmental change. This collection also provides insights into new techniques and approaches that will be valuable to furthering research on viral ecology.

## **ACKNOWLEDGMENT**

United States Department of Agriculture is an equal opportunity provider and employer.

## **REFERENCES**


**Conflict of Interest Statement:** Heather K. Allen declares that the editing of this volume was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. Stephen T. Abedon undertakes paid consulting for phage therapy-associated businesses but was not involved in editing the phage therapy-associated chapter in this volume (Meaden and Koskella, 2013).

*Received: 30 October 2014; accepted: 25 November 2014; published online: 12 December 2014.*

*Citation: Allen HK and Abedon ST (2014) Virus ecology and disturbances: impact of environmental disruption on the viruses of microorganisms. Front. Microbiol. 5:700. doi: 10.3389/fmicb.2014.00700*

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

*Copyright © 2014 Allen and Abedon. 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.*

## That's disturbing! An exploration of the bacteriophage biology of change

#### *Heather K. Allen1 \* and Stephen T. Abedon2 \**

*<sup>1</sup> Agricultural Research Service, United States Department of Agriculture, National Animal Disease Center, Ames, IA, USA*

*<sup>2</sup> Department of Microbiology, The Ohio State University, Mansfield, OH, USA*

*\*Correspondence: heather.allen@ars.usda.gov; abedon.1@osu.edu*

*Edited by:*

*John R. Battista, Louisiana State University and A & M College, USA*

**Keywords: virus ecology, ecological disturbance, phage-bacteria interactions**

The nature of life is change. Organisms change developmentally, morphologically, and physiologically, and they also modify their environments in the process. Such change can be gradual, abrupt, or even imperceptible components of steady states. Temporally, change can range from deterministic and regular processes to stochastic and unusual events. Unusual changes can result from relatively frequent occurrences, such as particularly devastating storms or instead highly uncommon events such as volcanic eruptions, ice ages, and asteroid impacts.

The scale of environmental change can be local or global, seasonal or climatic. The change can be caused by organisms, including as incidental consequences of their activities, such as our own anthropogenic influences on ecosystems and climate. Responses of organisms to change ranges from behavioral to physiological to developmental, and can occur over ecological as well as evolutionary time scales. The study of ecologically relevant physiological responses is described as physiological ecology, or ecophysiology, while evolutionary responses represent Darwinian evolution.

Among effectors of change are parasites, including viruses. In considering the viruses of bacteria—bacteriophages or phages change to hosts can vary from devastating lytic infections to simple genetic modification via lysogeny. In between are phages that are released from bacteria chronically, with productively infected bacteria continuing to replicate. Phage-induced change also can range from seemingly cosmetic chromatin rearrangements, as triggered by phage T2 during infection of *Escherichia coli* (Murray et al., 1950) to lysogenic conversion as can result in phage-encoded augmentation of bacterial pathogenesis (Addy et al., 2012). Bacteria also can change in response to phage infection, becoming immune to infecting phages by producing virus-specific interfering RNAs that are associated with CRISPR/Cas systems (Richter et al., 2012).

Change within the context of phages themselves is often more subtle. Phage virions can diffuse or be moved between microenvironments or ecosystems, resulting in changes in abiotic conditions. These environmental changes can cause physiological changes to their adsorption abilities (Conley and Wood, 1975). Phage physiology also changes dramatically as phages transition from their virion or free state to that of infecting bacteria, and then again from their phage-infecting form back to the free virion state. These transitions sometimes include a state of quiescence called pseudolysogeny. These pseudolysogenic states can be brought about by environmental conditions, such as bacterial starvation, and might help to promote phage survival (Ripp and Miller, 1997).

Examination of populations, communities, and entire ecosystems reveals that phages play integral roles in both eliciting and responding to disturbances (**Figure 1**). Phage biology and the relative impact of phages on bacterial populations can change, particularly as phage densities increase from those causing lower versus higher multiplicities of infection (Abedon, 2012). Environmental change in turn can impact viral population densities, including in terms of antibiotic induction of prophages (Allen et al., 2011). Bacterial fitness can change not just with phage quantity but also with phage quality, with greater bacterial fitness costs potentially associated with bacterial evolutionary responses to predation by multiple versus individual phage types (Koskella et al., 2012).

Selection on phages acts primarily on host acquisition, on rates of progeny production, and on survival until host acquisition again becomes a possibility. Suggesting an existence of tradeoffs associated with the optimization of these organismal properties, phage infection strategies may change in their effectiveness, pleiotropically, as phages change from infecting one bacterial strain to another (Duffy et al., 2006). Changes in the abundances of phages and other viruses in complex communities can occur seasonally in estuarian habitats (Winget et al., 2011) or in halophilic viral communities in response to environmental stressors (Santos et al., 2011). Phage abundance also can vary as a function of bacterial abundance, and in turn the cost to bacteria of phage sensitivity can increase as a function of phage abundance. Bacterial existence at high densities thus can result in phage-induced catastrophic changes in bacterial densities, a phenomenon that has been dubbed, "Killing the winner" (Rodriguez-Valera et al., 2009).

On the level of ecosystems, phages can be key contributors to the mineralization of nutrients as they solubilize host bacteria via lysis. As such they contribute to the primary ecological process of soils, that of decomposition and decay. In aquatic environments, phages potentially impact global carbon cycling by short circuiting the movement of carbon and energy to heterotrophic bacteria rather than from cyanobacteria to consumer eukaryotes

(Wilhelm and Suttle, 1999). In addition, phages can be added deliberately to environments to motivate change, as seen with phage-mediated biocontrol or phage therapy.

Change thus represents an ongoing and intrinsic aspect of phage biology, with phages both affecting and effecting organismal- population-, community-, ecosystem-, and even global environmental change. In this Research Topic we consider

## **REFERENCES**


environment-sensing device. *Proc. Natl. Acad. Sci. U.S.A.* 72, 3701–3705. doi: 10.1073/pnas.72. 9.3701


especially the impact of environmental change on communities and ecosystems as those changes may be propagated through phages, gene transfer agents, and viruses of other microorganisms.

## **ACKNOWLEDGMENTS**

USDA is an equal opportunity provider and employer.


virioplankton production within an estuarine ecosystem. *Proc. Natl. Acad. Sci. U.S.A.* 108, 11506–11511. doi: 10.1073/pnas.1101907108

*Received: 05 September 2013; accepted: 13 September 2013; published online: 14 November 2013.*

*Citation: Allen HK and Abedon ST (2013) That's disturbing! An exploration of the bacteriophage biology of change. Front. Microbiol. 4:295. doi: 10.3389/ fmicb.2013.00295*

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

*Copyright © 2013 Allen and Abedon. 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.*

#### Evolvability and robustness in populations of RNA virus -6

## *Daniel Goldhill , Angela Lee†, Elizabeth S. C. P. Williams and Paul E. Turner\**

*Department of Ecology and Evolutionary Biology, Yale University, New Haven, CT, USA*

#### *Edited by:*

*Stephen T. Abedon, The Ohio State University, USA*

#### *Reviewed by:*

*Richard Heineman, Kutztown University, USA John Dennehy, Queens College of the City University of New York, USA*

#### *\*Correspondence:*

*Paul E. Turner, Osborn Memorial Laboratory, Department of Ecology and Evolutionary Biology, Yale University, 165 Prospect St., PO Box 208106, New Haven, CT 06520-8106, USA e-mail: paul.turner@yale.edu*

#### *†Present address:*

*Angela Lee Hospital Luis Vernaza, Guayaquil, Ecuador*

Microbes can respond quickly to environmental disturbances through adaptation. However, processes determining the constraints on this adaptation are not well understood. One process that could affect the rate of adaptation to environmental perturbations is genetic robustness, the ability to maintain phenotype despite mutation. Genetic robustness has been theoretically linked to evolvability but rarely tested empirically using evolving populations. We used populations of the RNA bacteriophage φ6 previously characterized as differing in robustness, and passaged them through a repeated environmental disturbance: periodic 45◦C heat shock. The robust populations evolved faster to withstand the disturbance, relative to the less robust (brittle) populations. The robust populations also achieved relatively greater thermotolerance by the end of the experimental evolution. Sequencing revealed that thermotolerance occurred via a key mutation in gene P5 (viral lysis protein), previously shown to be associated with heat shock survival in the virus. Whereas this identical mutation fixed in all of the independently evolving robust populations, it was absent in some brittle populations, which instead fixed a less beneficial mutation. We concluded that robust populations adapted faster to the environmental change, and more easily accessed mutations of large benefit. Our study shows that genetic robustness can play a role in determining the relative ability for microbes to adapt to changing environments.

#### **Keywords: genetic robustness, evolvability, thermotolerance, bacteriophage, experimental evolution**

## **INTRODUCTION**

Viruses are often capable of very rapid molecular evolution, allowing adaptation to new hosts, and other novel challenges (Duffy et al., 2008; Wasik and Turner, 2013). Although viruses can quickly adapt in response to environmental changes, possible constraints on viral adaptation are seldom studied (Burch and Chao, 2000; Turner and Elena, 2000). One trait that may affect the rate of adaptation in evolving populations is genetic robustness: the capacity to maintain phenotype despite perturbation from underlying mutations (de Visser et al., 2003; Wagner, 2008; Draghi et al., 2010; Masel and Trotter, 2010). Because robustness buffers mutational effects, it may seem to be the antithesis of evolvability; however, robustness and evolvability may instead positively correlate (McBride et al., 2008; Draghi et al., 2010). For example, relatively robust proteins have greater structural stability; the effects of mutations that destabilize surface residues do not disrupt the core fold in robust proteins, providing an advantage for evolution of novel secondary functions (Bloom et al., 2006; Tokuriki et al., 2008; Tokuriki and Tawfik, 2009). More generally, theory shows that genetically robust populations contain large neutral networks of genotypes that span broad genotypic space, affording greater access to novel phenotypes following mutation, and allowing relatively robust populations to be more evolvable (Wilke and Adami, 2003; Draghi et al., 2010; Wagner, 2011). This link between evolvability and robustness is experimentally shown for RNA secondary structures and proteins, but rarely for evolving biological populations (Bloom et al., 2006; Elena and Sanjuán, 2008; Draghi et al., 2010).

Empirical studies demonstrating the relationship between robustness and evolvability at the level of populations are difficult for two main reasons: firstly, it can be difficult to identify or construct populations that differ in robustness (though see Montville et al., 2005; Sanjuán et al., 2007; Coleman et al., 2008) and secondly, the timespan necessary to conduct evolution experiments is often very long (Blount et al., 2008; Kawecki et al., 2012). RNA viruses offer an experimental system that can overcome these problems (Elena, 2012). These viruses have a high mutation rate and are easily cultured in the laboratory so the effects of differences in genetic robustness can be studied in a short timespan (Montville et al., 2005; Sanjuán et al., 2007). Furthermore, the small genome sizes of RNA viruses offer the possibility of identifying the specific genetic architectures leading to robustness and how they affect evolvability. Reverse genetic techniques allow large-scale manipulations of the genetic code of viruses, such as switching codons whilst leaving the amino acid sequence intact. Codon switching can increase the percentage of mutations that are non-synonymous, thereby decreasing the number of neutral neighbors and reducing the robustness of a virus population (Lauring et al., 2012).

An alternative way to manipulate robustness in viruses is by varying coinfection level (Montville et al., 2005; Gao and Feldman, 2009). Coinfection allows for complementation, where the effects of harmful mutations are buffered because viruses with deleterious or inactive proteins can be complemented by a coinfecting virus with the beneficial or active protein (Froissart et al., 2004; Aaskov et al., 2006; Gao and Feldman, 2009). Thus, passaging virus populations under high levels of coinfection (and hence, complementation) should reduce selection to maintain robustness at the level of an individual virus, because the environment (coinfection) provides the mutational buffering. Montville et al. confirmed this idea by evolving three populations of the dsRNA bacteriophage φ6 for 300 generations at high vs. low multiplicity of infection (MOI), the ratio of infecting viruses to cells (Montville et al., 2005). When clones from the populations were used to found lineages subjected to mutation accumulation (successive bottlenecking that causes mutations to fix via drift), it was revealed that the high-MOI-evolved viruses showed greater variance in the fitness effects of accumulated mutations (reduced robustness) and these populations were termed "brittle." In contrast, lineages founded by clones isolated from low-MOI-evolved populations showed lesser variance in fitness effects of accumulated mutations, defining these populations as "robust." Consistent with these findings, Dennehy et al. later showed that greater sequence diversity existed in the low-MOI populations compared to their high-MOI counterparts (Dennehy et al., 2013). This observation suggested that the more robust, low MOI populations had greater genetic variation because they contained a larger neutral network of genotypes.

To test whether robustness imparted an evolvability advantage in phage φ6, McBride et al. used clones from the robust and brittle populations to found lineages that were subjected to a novel environment: strong selection pressure (high mortality) caused by periodic exposure to 45◦C heat shock (McBride et al., 2008). After passaging the populations through 10 rounds of selection (5 min heat shock, interspersed by 5 generations of growth under normal conditions), improved thermotolerance was observed. However, lineages founded by clones from the robust populations showed a significantly greater increase in thermotolerance, on average, than those founded by clones from the brittle populations. This result indicated that robustness could enhance evolvability, at least in the particular novel environment. However, we note that other virus studies have not always found a positive relationship between evolvability and robustness (Cuevas et al., 2009; Tokuriki et al., 2009).

If viral populations exist in neutral networks, it may be problematic to test the link between robustness and evolvability using populations founded by individual clones (McBride et al., 2008). That is, a clone may not have enough time over the course of a short-term selection experiment to explore the entire neutral network, reducing any differences between small and large networks. Furthermore, studies have shown that rare variants within a population can disproportionately affect evolutionary trajectories (Blount et al., 2008). The current study addresses these caveats by harnessing the same study system to examine whether robustness aids adaptation to the novel environment, but employing genetically-variable population samples (rather than clones) to found test lineages. Similar to McBride et al. (2008), we measured thermotolerance as the key phenotypic change to assess whether robustness enhances evolvability. However, in the current study we monitored phenotypic changes over time (adaptive trajectories) due to the specific prediction that larger neutral networks should foster better access to key mutations that speed adaptation. In addition (unlike the former study), we measured molecular changes via consensus genome sequencing of experimentally evolved populations, to identify whether different beneficial mutations fixed in robust vs. brittle populations.

## **MATERIALS AND METHODS STRAINS AND CULTURE CONDITIONS**

A single colony of *Pseudomonas syringae pv. phaseolicola* (American Type Culture collection #21781) was taken from a source plate each day and cultured overnight at 25◦C in 10 ml of Luria Broth (LB). The six phage populations of φ6 used in this study were previously described (Montville et al., 2005). Three populations were previously evolved for 300 generations under low multiplicity of infection (MOI = 0.002) and characterized as robust: L1–L3 (strains #PT578-PT580); whereas, three otherwise identical populations were evolved at high multiplicity (MOI = 5) and deemed brittle: H1–H3 (strains #PT581-PT583) (**Figure 1**). Viruses were grown at 25◦C in 3 ml of top agar (0.7% agar) containing 200 ul of an overnight bacterial culture, overlaid on LB agar plates (1.5% agar). Lysates were made by removing the top agar, centrifuging in 3 ml of LB and filtering (0.22 μm filter, Millipore) to remove bacteria. Viral lysates were stored in glycerol/LB mixture (2:3 by volume) at −20◦C.

## **SERIAL PASSAGE WITH HEAT SHOCK**

To heat shock a virus population, 50 ul of a diluted viral lysate was mixed with 3 ml of top agar and placed in a heating block at 45◦C for 5 min; then, 200 ul of overnight bacterial culture was added to the lysate and immediately plated as described above. Different dilutions of the virus lysate were heat-shocked and plated to ensure a resulting plate containing <sup>∼</sup>10<sup>4</sup> plaques, which constituted a controlled bottleneck size of the evolving virus population. Lower dilutions with countable numbers of plaques confirmed the bottleneck size was <sup>∼</sup>10<sup>4</sup> and non-heatshocked control dilutions were also plated alongside for comparison, to ensure that heat shock was causing virus mortality. After 24 h, the 10<sup>4</sup> plaques created a "lacy lawn" (highly overlapping plaques), which yielded an extremely high-titer lysate. The passage was repeated by using the fresh lysate and naïve (noncoevolving) bacteria. Each population was subjected to this daily passage for 10 total days, which was equivalent to 50 generations (5 generations per day) (Turner and Chao, 1998) with 5 min heat shock imposed every fifth generation (**Figure 1**). A sample of each lysate (evolving population) was frozen at each daily passage.

## **SEQUENCING**

High-titer lysates were created of the populations to be sequenced and RNA extracted using QIAamp Viral RNA minikits (Qiagen). RNA was converted into cDNA using SuperScript II (Invitrogen) and then used as template for PCR (primers available on request.) Sanger Sequencing was performed by the Yale Science Hill DNA Analysis Facility. Sequences were manually inspected and analyzed using CLC DNA Workbench 6. Sequences are available in Genbank under the following accession numbers: KF996287- KF996304.

#### **SURVIVABILITY ASSAY AND THERMAL NICHE ANALYSIS**

75 ul of diluted viral lysate was placed in a PCR tube and heated for 5 min at a test temperature (42.5–47.5◦C) in a pre-heated Eppendorf Thermocycler. A 50 ul aliquot from the heat shocked lysate was plated, alongside an otherwise identical control plate ("mock" heat shock: 5 min incubation at 25◦C). Survivability fraction was calculated as the number of plaques formed following heat shock divided by those formed on the control plate.

## **RESULTS**

## **THERMOTOLERANCE OF INITIAL POPULATIONS**

We tested whether the low-MOI evolved "robust" populations (L1, L2, L3) and high-MOI evolved "brittle" populations (H1, H2, H3) initially differed in survival following heat shock. To do so, each population was subjected to replicated (*n* = 4 or 5) 45◦C heat shock assays in top agar, as well as assays in a thermocycler at two different temperatures: 42.5 and 45◦C. Results (**Figure 2**) showed no statistical differences in the initial thermotolerance of robust and brittle populations [Two-Way ANOVA, *F*(1, <sup>12</sup>) = 0.25, *p* = 0.63]. The observed low survival confirmed that periodic 45◦C heat shock would create a strong selective pressure, and that initial survival did not differ according to prior ecological history (low vs. high MOI experimental evolution).

## **FASTER THERMOTOLERANCE ADAPTATION IN ROBUST POPULATIONS**

bars show standard error and data points are offset for clarity.

To test how robustness affected evolution of thermotolerance, we evolved the six populations in a novel heat shock environment using a constant bottleneck size of 10,000 pfu (plaque forming units). As population size is correlated with speed of evolution (Szendro and Franke, 2013), it was important to maintain equivalent population size among the robust and brittle experimental populations. The bottleneck size of 10,000 individuals was large enough to permit extensive genetic diversity, whilst allowing ∼5 generations of evolution each day as the populations grew to their maximum size that preceded the bottleneck (Turner and Chao, 1998). Previous experiments showed rapid evolution in phage φ6 populations when cultured at sizes comparable to our study (Burch and Chao, 1999). Montville et al. (2005) showed no fitness differences between robust and brittle populations used to initiate the current study, which suggests that these populations experienced comparable generation numbers during the prior experimental evolution, and that no pre-existing fitness bias favored an evolvability advantage in our experiment.

Following 10 days of evolution, all of the populations showed increased thermotolerance compared to their founding (day 0) population (2-tailed paired *<sup>t</sup>*-test, *<sup>T</sup>* <sup>=</sup> <sup>13</sup>.5, *d.f.* <sup>=</sup> 5, *<sup>p</sup>* <sup>&</sup>lt; <sup>4</sup>∗10−5, see **Figure 3**). As described in Dessau et al. (2012), phage φ6 strains with improved thermotolerance often show a "bull's-eye" plaque phenotype when grown on agar under normal conditions of 25◦C (see **Figure 3**). Three of the six evolved populations (i.e., L1, L3, H2) showed apparent fixation of the bull's-eye phenotype, which was the only morphotype observed on day 10 plate dilutions of these populations. In contrast, two of the evolved populations (i.e., L2, H1) were polymorphic for bull's-eye plaques, with L2 showing ∼50% bull's-eyes and H1 producing ∼10% bull's-eyes. Last, a single evolved population (H3) showed almost entirely clear (wildtype phenotype) plaques at the end of the experiment.

To test whether the rate of phenotypic adaptation differed between robust and brittle populations, samples from each population at intermediate days were tested for thermotolerance. Results (**Figure 3**) showed that robust population L1 had increased in thermotolerance after 3 days, whereas all other populations had not. After 6 days, all of the robust populations showed an increase in thermotolerance, compared to a minimal such increase in the brittle populations; this group-wise difference was statistically significant (2-tailed *t*-test, *T* = 3.74, *d.f.* = 4, *p* = 0.02). The bull's-eye phenotype allowed a visual estimate of the penetrance of thermotolerance mutations. Bull's-eye plaques were first seen in population L1 after the second passage and in populations L2 and L3 shortly thereafter. The bull's eye phenotype was apparently fixed in population L1 by the sixth passage. The bull's-eye phenotype reflected the observed thermotolerance data in population L1: a small increase in thermotolerance by day three coincided with low frequencies of observed bull's eyes, while the higher increase in thermotolerance by day 6 coincided with bull's eyes as the only visible phenotype. To confirm that there were no bull's-eyes present in the original population, we heat shocked the day 0 L1 population and plated multiple dilutions to visualize individual plaques; here we observed no bull's-eyes after screening ∼10,000 plaques.

#### **MOLECULAR EVOLUTION**

The genome size of phage φ6 is ∼13 kb and consists of three dsRNA segments (Small, Medium and Large). To investigate whether there were differences among the evolved populations in mutations leading to thermotolerance, we sequenced the endpoint (day 10) evolved populations. The consensus sequences revealed several polymorphisms but few fixed protein changes across the six populations (**Figure 4**). The only locus that showed changes shared among multiple populations was the gene for protein P5, the viral lysin. One non-synonymous change in P5, G2238T (a valine to phenylalanine mutation), was fixed in three populations (L1, L3, H2) and polymorphic in two other populations (L2 and H1). This mutation was previously shown to increase thermotolerance by stabilizing the enzyme under elevated heat, while simultaneously causing a bull's-eye plaque phenotype, indicating reduced virus growth under 25◦C conditions (Dessau et al., 2012). One other substitution in P5, A1857G (Lysine to Glutamic Acid), was also shared across two populations (H1, H3). Population L2 had a pre-existing polymorphism at position 2274 which remained in the evolved population. In addition, L2 showed two other non-synonymous polymorphisms in P5: G2229C and A2254G. Although sequencing showed that the G2238T mutation was not present in H3, bull's-eye plaques were observed at low levels (∼1%) in population H3. We chose four bull's-eye plaques from population H3 on days 9 and 10 for sequencing, to test whether the G2238T mutation was present. However, the G2238T mutation was not found in any of the chosen plaques. Instead, we observed two other mutations near the end of gene P5 in population H3 that presumably caused the bull's-eye phenotype.

#### **EVOLVED CHANGES IN THERMAL NICHE**

We tested the survivability of the evolved populations at three different temperatures to approximate a thermal reaction norm for each population (**Figure 5**). We found no significant difference at 42.5◦C, as all populations showed very high survival (2 tailed *t*-test, *T* = 1.36, *d.f.* = 4, *p* = 0.18). However, at 45 and 47.5◦C, the robust populations were significantly more thermotolerant than the brittle populations [Two Way ANOVA, *F*(1, <sup>8</sup>) = 10.27, *p* = 0.01]; there were no significant differences between the robust populations at these temperatures [Two Way ANOVA, *F*(2, <sup>12</sup>) = 2.18, *p* = 0.15].

As key mutations were shared between populations, we wanted to test whether the genetic background affected the phenotype and in particular, how the G2238T mutation affected phenotype because it was shared between both brittle and robust populations. The three populations in which the G2238T mutation was fixed- L1, L3 and H2- showed no significant difference in survival at any of the three temperatures suggesting that the mutation had a similar effect on all three genetic backgrounds. Brittle populations, which did not have the G2238T mutation (H3) or had the mutation in a low percentage of the population (H1) shared an alternative mutation (A1857G) and showed lower thermotolerance at 45 and 47.5◦C. The overall significantly lower survival at high temperatures of the brittle populations was caused by populations H1 and H3 and not by population H2, which did not differ from the robust populations.

#### **DISCUSSION**

We subjected three robust and three brittle populations of RNA phage φ6 to a novel environment—periodic heat shock selection for 10 days—after which all the populations were observed to increase in thermotolerance. However, the robust populations evolved thermotolerance earlier in the experiment, and were significantly advantaged in thermotolerance relative to most of the

plotted for day 10 populations. Squares and circles represent fixed mutations and polymorphisms respectively. Amino-acid changes are shown for non-synonymous mutations (filled symbols) and nucleotide changes are shown for all mutations in P5. Open symbols represent polymorphisms present in the starting population. Positions 2126 and 2242 on the small segment and 439, 1613, 3201, and 3561 on the medium segment are fixed mutations resulting from losses of polymorphism.

brittle populations. Thus, we concluded that robustness tended to promote evolvability, when populations of phage φ6 experienced the heat-shock selection conditions imposed in our study.

Theory states that populations containing a larger neutral network in genotype space should be more evolvable, because they can access a greater range of mutations from different points on the network (Wagner, 2008). This increased evolvability could occur either through faster access to key mutations- as there is no need for a permissive primary mutation before a secondary mutation evolves- or through access to a greater variety of mutations with large phenotypic effects. Both explanations relate to the findings in our experiment. All of the robust populations (L1, L2, L3) evolved higher thermotolerance by day 6 than the brittle (H1, H2, H3) populations (**Figure 3**). The robust populations were able to access a key thermotolerance mutation (G2238T) earlier than the brittle populations, confirmed by their greater increase in thermotolerance by day 6 and by the phenotypic appearance of bull's-eye plaques that are associated with the G2238T substitution (Dessau et al., 2012). It is very unlikely that the G2238T mutation was present in any of the founding populations, as the plaquing phenotype was easily distinguishable and was not observed when founding populations were screened in survival assays (i.e., where genotypes producing bull's eye plaques are strongly positively selected). Rather, bull's-eye plaques were not observed until after day 2 (generation 10) in L1 and even later in other populations.

The only mutations shared by multiple populations were in gene P5 on the Small segment. P5 encodes the viral lysin, a lytic transglycosylase enzyme used by phage φ6 to penetrate the bacterial cell wall during entry and exit of the host cell (Mindich and

day 10 population were heat shocked in a thermocycler at three different temperatures and then plated (*n* = 3). A control was also plated and percent survival was calculated. Error bars show standard error and data points are offset for clarity.

Lehman, 1979; Dessau et al., 2012). The lytic enzyme resides in between the protein shell (P8) surrounding the nucleocapsid and the lipid envelope, but little is known about interactions between P5 and other structural elements of phage φ6.

The G2238T mutation in P5 was seen in every population except for H3 and was present at low frequency in H1; in contrast, these two populations had a different mutation in P5: A1857G. It seems likely that this mutation also increases thermotolerance because it occurred following selection in two independent populations, it is a non-synonymous change from a positively to a negatively charged amino acid, and it is the only evolved change in gene P5 of population H3, which showed increased thermotolerance but no bull's-eye phenotypes. The A1857G substitution seemed to be a less-beneficial P5 mutation, because it increased thermotolerance to a lesser degree than the G2238T mutation (cf. **Figures 4**, **5**).

The large number of polymorphisms in population L2 suggested the presence of at least two subpopulations. Intriguingly, L2 had a pre-existing polymorphism in gene P5 at position 2274, which remained in the evolved population. This polymorphism might have a negative epistatic interaction with the G2238T mutation, perhaps explaining why G2238T did not fix in this population. L2 had two other polymorphisms in gene P5 (G2229C and A2254G) that may have evolved because the G2238T mutation was selected against in a subpopulation. This will be the subject of further experiments.

All of the robust populations accessed a more beneficial mutation than population H3, in the evolutionary time allowed. Initially, population H3 had no mutations in the amino acid sequence of gene P5 that could have restricted the evolution of the G2238T mutation, whereas H1 and L2 both had mutations close to this locus, such that negative epistasis may explain why the mutation did not fix in either population. Population H3 showed a synonymous mutation in P5; although this change could not affect the protein structure, it could still be negatively epistatic by altering RNA structure. Brittle populations occupy smaller genotypic networks so it is not surprising that there was variation between populations in which mutations fixed: population H2 was in a network that allowed the evolution of G2238T, whereas population H3 may have been in a network that prevented the evolution of this mutation or promoted a less beneficial mutation. When we sequenced bull's-eye plaques to search for the G2238T mutation in H3, we did not find the mutation. Its absence suggests this highly beneficial mutation could not arise or was selected against in H3. When the G2238T mutation evolved, there was no evidence that the phenotypic effects differed between robust and brittle populations depending on the genetic background, as populations L1, L3, and H2 did not significantly differ in survival across a range of temperatures (**Figure 5**).

Our results complement and extend those of McBride et al. who showed that lineages founded by robust clones evolved greater thermotolerance than those founded by clones from brittle populations (McBride et al., 2008). We showed that the magnitude change in thermotolerance reported by McBride et al. was likely due to either slower evolution in brittle clones or failure in some clones to find a key mutation. Our results also showed that some brittle populations were able to find the same mutation and catch up to the robust populations, as after 10 days there was no significant difference in phenotypic effect. However, some populations appeared unable to evolve a key mutation altogether. This highlighted the importance of starting with a population as opposed to clones as it demonstrated that in some brittle populations, the loss of intrinsic robustness meant the population had shifted to a less evolvable network. That is, we interpret that 300 generations of evolution at high MOI resulted in the loss of intrinsic (genotypic) robustness, because these high-MOI-evolved populations experienced frequent complementation among coinfecting viruses, an "environmental effect" which should have relaxed selection to maintain individual-level robustness. The resulting evolved brittleness of these populations would cause a concomitant shift to a brittle, less evolvable network. In the current study, initiating experimental evolution with population samples would allow us to test whether such brittle populations were indeed on a less evolvable network, whereas starting the experimental evolution with a single clone could only demonstrate that a single point on the network was less evolvable. Starting with populations also allowed us to test the possibility that rare mutations present in the initial populations could impact the evolutionary trajectory. However, we found no evidence of this in our experiment as all the mutations in P5 occurred *de novo* and were not initially present. The initial diversity within the populations of the current study was necessarily limited by the size of the serial-passage bottleneck (*N* ≈500) they experienced in the study immediately preceding ours (Montville et al., 2005). Although we did not explicitly measure this initial diversity prior to starting the current study, the variation must have exceeded that of a single clone; this differing manipulation of initial variation highlights the key difference between our study and that of McBride et al. (2008), when examining relative evolvability of the populations for increased thermotolerance. Beyond this conservative difference in starting conditions, we note that even higher initial diversity in the founding populations could have led to greater differences between our results and those of McBride et al.

Our results differ from those involving a different RNA virus, Vesicular Stomatitis Virus (VSV), in which no link was found between robustness and evolvability (Cuevas et al., 2009). However, unlike our study, the VSV study founded lineages using clones, rather than populations. It is possible that starting from clones removed initial variation and did not allow the full exploration of the neutral network. Cuevas et al. suggested that in robust clones, beneficial mutations might have reduced phenotypic effect, whereas in our study there was no evidence that beneficial mutations had less phenotypic effect in robust populations. Further work will be necessary to resolve conflicting results for different types of RNA viruses and to assess the generalizability of the link between robustness and evolvability.

The mechanism allowing robustness to affect evolvability in our study remains unclear. The robust populations were shown to harbor greater genetic diversity than the brittle populations (Dennehy et al., 2013), which is consistent with the definition of robustness (i.e., expectation that more genotypes of equal fitness can exist in a robust population than in a brittle one). This greater diversity of genetic backgrounds in a robust population affords increased possibility for *de novo* evolution of positive epistasis with new mutations, and could partly explain the relatively faster rate of thermotolerance evolution in robust populations. There were no easily identifiable permissive mutations in the robust populations that were required before the evolution of thermotolerance, as there were no other mutations seen in multiple populations. However, there were candidate mutations in the brittle populations (as well as a subpopulation of L2) that could have restricted the evolution of the most beneficial mutation G2238T. Although the phenotypic effects were similar between the robust and one of the brittle populations at day 10, the increased thermotolerance of the robust populations by day 6 implies that, under direct competition, the robust populations would have a significant advantage.

One potential problem with our experiment is that genetic robustness and environmental robustness have been shown to be linked, often called plastogenetic congruence, though data is sparse at the level of organisms (Ancel and Fontana, 2000; Meyers et al., 2004; Novella et al., 2013). This means that thermotolerance, a form of environmental robustness, might not be an appropriate trait. We controlled for this by testing initial survival at both 42.5 and 45◦C (**Figure 2**). We found no evidence that the robust populations were pre-adapted to be more thermotolerant as there was no survival difference between robust and brittle populations at 42.5 and 45◦C. However, these measurements were taken for the entire population and there may have been individual clones, which were more thermotolerant (Ogbunugafor et al., 2009). Future work could test if our evolvability results are generalizable to other traits.

It is important to study robustness and evolvability using biological populations, as the results may differ when examining more complex systems than experimenting on individual proteins *in vitro*. Our study is one of the first to experimentally examine the effects of robustness on evolvability using populations that differ in robustness, but future work should extend to other viruses. Robustness of viruses might be important in determining many important viral traits such as the ability to switch hosts and the degree to which a virus is pathogenic (Ogbunugafor et al., 2010; Lauring et al., 2012; Remold, 2012). Robustness may also be a determinant of the success of antiviral therapies such as lethal mutagenesis (Bull et al., 2007). We know little about robustness of viruses in natural systems. If *de novo* mutations are required to respond to environmental changes as opposed to selection on standing variation, relatively robust populations of viruses are likely to have an advantage. A better understanding of robustness will hopefully elucidate how organisms respond to environmental challenges and evolve novel functions.

## **AUTHOR CONTRIBUTIONS**

Daniel Goldhill and Paul E. Turner designed the study, Daniel Goldhill and Angela Lee conducted the experimental evolution, Daniel Goldhill and Elizabeth S. C. P. Williams conducted thermotolerance assays, Daniel Goldhill and Paul E. Turner analyzed the data and prepared the manuscript. All authors read and approved the final manuscript.

## **ACKNOWLEDGMENTS**

We would like to thank Nadya Morales, Jason Shapiro and other members of the Turner Lab for their help. Funding was provided by NSF grant DGE-1122492 to Daniel Goldhill and NSF grant DEB-1021243 to Paul E. Turner. In addition, Daniel Goldhill was supported by the NIH Genetics Predoctoral Research Training Grant (T32-GM007499).

## **REFERENCES**


**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.

*Received: 08 November 2013; accepted: 19 January 2014; published online: 05 February 2014.*

*Citation: Goldhill D, Lee A, Williams ESCP and Turner PE (2014) Evolvability and robustness in populations of RNA virus* -*6. Front. Microbiol. 5:35. doi: 10.3389/fmicb. 2014.00035*

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

*Copyright © 2014 Goldhill, Lee, Williams and Turner. 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.*

## Evidence for metaviromic islands in marine phages

## *Carolina Megumi Mizuno , Rohit Ghai and Francisco Rodriguez-Valera\**

*Evolutionary Genomics Group, Departamento de Producción Vegetal y Microbiología, Universidad Miguel Hernández, Alicante, Spain*

#### *Edited by:*

*Heather K. Allen, National Animal Disease Center, USA*

#### *Reviewed by:*

*Amparo Latorre, University of Valencia, Spain Karin Holmfeldt, Linnaeus University, Sweden*

#### *\*Correspondence:*

*Francisco Rodriguez-Valera, Evolutionary Genomics Group, Departamento de Producción Vegetal y Microbiología, Universidad Miguel Hernández, San Juan de Alicante, 03550 Alicante, Spain e-mail: frvalera@umh.es*

Metagenomic islands (MGIs) have been defined as genomic regions in prokaryotic genomes that under-recruit from metagenomes where most of the same genome recruits at close to 100% identity over most of its length. The presence of MGIs in prokaryotes has been associated to the diversity of concurrent lineages that vary at this level to disperse the predatory pressure of phages that, reciprocally, maintain high clonal diversity in the population and improve ecosystem performance. This was proposed as a Constant-Diversity (C-D) model. Here we have investigated the regions of phage genomes under-recruiting in a metavirome constructed with a sample from the same habitat where they were retrieved. Some of the genes found to under-recruit are involved in host recognition as would be expected from the C-D model. Furthermore, the recruitment of intragenic regions known to be involved in molecular recognition also had a significant under-recruitment compared to the rest of the gene. However, other genes apparently disconnected from the recognition process under-recruited often, specifically the terminases involved in packaging of the phage genome in the capsid and a few others. In addition, some highly related phage genomes (at nucleotide sequence level) had no metaviromic islands (MVIs). We speculate that the latter might be generalist phages with broad infection range that do not require clone specific lineages.

**Keywords: population genomics, phages, metaviromes, host recognition, marine phages, virome**

## **INTRODUCTION**

One of the major conundrums of microbiology is the population biology of bacteria and their phages (Rodriguez-Valera, 2002). Pure culture, the work-horse of bacteriology, imposes a strong prejudice about how bacteria really are in nature (Hugenholtz et al., 1998). Owing to their asexual reproduction the population in a pure culture is essentially clonal, composed of nearly identical cells. Likewise, in an infectious disease, typically the pathogen is a clone and repeated isolation of the same clone is at the root of modern epidemiology (Stenderup and Orskov, 1983). However, we have very little knowledge of what the situation is in natural environments such as soil or waters. One approach to study the population biology of a microbe is to isolate and compare many strains from the same sample (Tettelin et al., 2005; Jacobs-Sera et al., 2012). However, this approach is fraught with challenges. The isolation process can be heavily biased (Rodriguez-Valera, 2004) and the same clone may be retrieved many times. To get a realistic picture many strains might be required, and the cost of sequencing and analyzing their genomes may be disproportionate to the information gain.

On the other hand, high throughput metagenomics provides a composite of the individual cells present in a sample that can be gathered using the genome of one isolate for retrieval of the homologous regions in the clonal lineages present in the population. This virtual experiment has already been done repeatedly for cellular microbes revealing very interesting patterns (Coleman et al., 2006; Legault et al., 2006; Cuadros-Orellana et al., 2007; Rodriguez-Valera et al., 2009). Typically, a large part of the strain genome is covered at a very high identity. Actually, several 100% identical fragments are found for most of the core genome. The retrieval of identical fragments could have been expected, particularly when the isolate comes from the same sample as the metagenome. However, normally there is also a cloud of reads with similarities within the range from 100 to ca. 95% (Caro-Quintero and Konstantinidis, 2012). This similarity cloud reflects the natural variation of the clonal frames that coexist within the population (López-Pérez et al., 2013). Even less expected was the appearance of selected tracts of the reference genome that recruit much less and sometimes not at all. These regions have been named metagenomic islands (MGIs) and overlap mostly (but not completely) with the flexible genome (defined when comparing different strain genomes) (Coleman et al., 2006; Cuadros-Orellana et al., 2007; Pašic et al., 2009 ´ ). Overall, these data can be explained by a high degree of concurrent genomic diversity in prokaryotic populations (Gonzaga et al., 2012; López-Pérez et al., 2013). The finding of such high diversity is contradictory to the classical paradigm derived from laboratory enrichments indicating that clonal sweeps happen regularly in prokaryotic populations narrowing the genetic diversity (Atwood et al., 1951; Cohan and Koeppel, 2008). Particularly, it raises the question of how selection events that favor a single high-fitness clone (clonal sweeps) are prevented from decreasing clonal diversity. Considering these issues and the high diversity of concurrent clonal lineages found in the environments analyzed, a model was proposed in which high clonal diversity was maintained by the predatory pressure exercised by bacteriophages (Rodriguez-Valera et al., 2009). The constant-diversity model (C-D) posits that several clonal lineages coexist, enriching the genetic wealth of the population, and these are kept under control by the stabilizing influence of a similar diversity of viruses that prey on the population. It is actually a revival of the classic "kill-thewinner" model proposed years ago (Thingstad and Lignell, 1997; Thingstad, 2000) to explain the plankton paradox, but this time applied to clonal frames within a single species, rather than to different species with similar niches.

One of the predictions of C-D is that diverse populations of phages preying on the same species coexist in a single environmental niche. Actually, the clonal diversity of phage populations has been also approached by isolating multiple phages and indicated that indeed many lineages coexist at a single time and place, for example, roseophages (Angly et al., 2009), cyanophages (Labrie et al., 2013), and *Alteromonas* phages (Garcia-Heredia et al., 2013) among others. However, drawbacks similar to those encountered in prokaryotic pure culture methods limit the weight of the evidence. Rodriguez-Brito et al. studied the diversity of phages in individual saltern ponds by a viral metagenomic approach and found evidence for a high diversity that was maintained through seasons and years, even though within a short time scale, rank switching could be detected (Rodriguez-Brito et al., 2010). Along similar lines, a high conserved diversity of phages has been detected in a hypersaline lake (Emerson et al., 2013). However, given the scarcity of phage genomes isolated from the same environments, a phage equivalent of the MGIs detected in bacteria and archaea has been hitherto missing. There is only the report of Garcia-Heredia et al, indicating that some host recognition genes (glucanases) under-recruited in metaviromes from hypersaline waters (Garcia-Heredia et al., 2012).

For prokaryotic genomes, the presence of MGIs has been interpreted as indicative of the presence of different clonal lineages with variations at some gene clusters diluting the recruiting efficiency of these regions. These gene clusters are diverse but include an abundance of outer cellular structures that are candidate recognition target of phages, such as the O-chain polysaccharide synthesis genes, that were always among the least recruiting in the genome. We wondered if the reciprocal was also true. That is, that the regions that were different among the viral lineages were involved in host recognition as some previous evidence seemed to indicate (for example, Angly et al., 2009). In this work, taking advantage of the availability of a large set of phage contigs from the Mediterranean Deep Chlorophyll Maximum (MedDCM), by cloning metagenomic DNA in fosmids and a metavirome from the same location (Mizuno et al., 2013), we examine patterns of variability and the presence of the viral equivalent of MGIs or metaviromic islands (MVIs) in the most abundant phages. The results show a remarkable degree of clonal diversity of concurrent phages, even more so than the prokaryotic host and affecting different genes, with some (but not all) involved in host recognition.

## **MATERIALS AND METHODS**

## **SAMPLING, FOSMID, AND METAVIROME SEQUENCING**

The sampling, fosmid library construction, assembly and recovery of phage genomes has been described previously (Ghai et al., 2010; Mizuno et al., 2013). Briefly, a fosmid library of ∼13000 clones was constructed from the DNA from a 0.22μm filter from a sample from the Mediterranean DCM. The sampling date was October 15, 2007, depth 50 m and location was off the coast of Alicante, Spain (38◦4 6.64N 0◦13 55.18W). DNA from ∼6000 metagenomic fosmids was sequenced in 24 batches using Illumina. Each batch was assembled independently. 1148 phage derived contigs were identified using sequence based approaches and classified into groups based on high sequence identity. Complete genomes of phages were termed as Complete Genome Representatives (CGRs) and incomplete contigs related to CGRs were termed Complete Genome Fragments (CGFs). All 1148 contigs assembled were submitted to DDBJ and are available using the accession numbers AP013358-AP014505.

Sampling and sequencing of metavirome MedDCM-Vir has also been described previously (Mizuno et al., 2013). Briefly, the metaviromic biomass originates from another sample of 20 L from the Mediterranean DCM (sampling date August 29, 2011, depth 65 m, from the same location as above). After filtration through 0.22μm, phages were concentrated using tangential flow filtration (TFF) with a 30 kD polyethersulfone membrane and multiple amplification displacement (MDA) was performed. The resulting DNA was sequenced in one third of an Illumina lane. The metavirome has been deposited in NCBI SRA with the Bioproject number PRJNA210529.

### **SELECTION OF PHAGE CONTIGS AND MVI DETECTION**

From a total of 1148 phage contigs described previously (Mizuno et al., 2013), we have focused on 208 that were identified as complete genomes (referred to as CGRs) and those highly related to them, albeit incomplete (524 Complete Genome Fragments, CGFs). The CGRs, by comparison to reference phage genomes were further classified into 21 groups (G1, G2 etc.) that are akin to high level taxa (such as family) (Mizuno et al., 2013). In order to reliably identify MVIs we first identified the most abundant CGRs and CGFs in a MDA-amplified metavirome obtained from the same habitat and location 4 years later (Mizuno et al., 2013). Fragment recruitment was done at very high nucleotide identity (*>*98% over at least 50 nucleotides) and only those with a minimum median base coverage of at least 5× and a value of reads per Kb of sequence per Gb of metavirome (RPKG) *>*1 were used for further analysis (total of 224 phage contigs from which 51 are CGRs and 173 CGFs). In these 224 selected contigs, regions that had a coverage of less than 20% of the median coverage of the entire contig and were *>*500 bp in length were considered as MVIs (Data S1).

#### **PROJECTING FRAGMENT RECRUITMENT ON PROTEIN STRUCTURE**

The base coverages of each C1q gene were projected onto the amino acid sequence of the same using the average of the read coverage for each base for each codon (that is, data for each amino acid position was computed by averaging the read coverages of all three bases in the codon). Protein structure prediction for the C1q domain was performed using the I-TASSER server (Zhang, 2008). The results of the mapping of the read coverages onto the predicted protein structure were visualized using PyMOL version 1.2r1 (Delano, 2002).

## **RESULTS AND DISCUSSION IDENTIFICATION OF MVIs**

We have focussed our analyses on 224 abundant contigs (RPKG *>*1) of which 51 were CGRs and 173 were CGFs. The results of the CGR recruitment analysis against the metavirome are summarized in **Figure 1**. These phage genomes and the metavirome originate from the same location, but 4 years apart. All but one of the predicted hosts for these abundant phages were Alphaproteobacteria, e.g., *Ca.* Pelagibacter (SAR11) and *Ca.* Puniceispirillum (SAR116) (Mizuno et al., 2013). This was not unexpected, as these cellular taxa are among the most abundant in the marine habitat (Delong et al., 2006; Ghai et al., 2010; Oh et al., 2010). However, there are also very abundant phage genomes for which the host is as yet unidentified (for example, those belonging to group G16) (**Figure 1A**, Data S1). We used a very high recruitment threshold (98% identity over at least 50 bp), designed to recruit only reads that belong to the same phage clone (accounting for an error rate of ca. 1%). It is a well proven fact that MDA, used to amplify the DNA of the metavirome, introduces biases. However, we believe the under-recruiting regions (by the large difference in coverage required) actually reflect reduced frequencies of the genes within the population.

We identified a total of 165 contigs (40 of them CGRs) out of the 224 that contained at least one MVI (See Methods, Data S1, S2). Among the CGRs where MVIs were identified, the size and recruitment levels of each MVI was variable, ranging from very small, barely covering a protein domain in a gene (close to our established lower island size limit of 500 bp) to up to 4.6 kb, covering sometimes three or four genes (**Figure 1**, Data S1). The fraction of the genome represented by MVIs was on average ca. 13%, (ranging from 1.53 to 29%). **Figure 1** shows the coverage plots of the most recruiting individual CGRs from groups G15, G16, G17, and G19. The CGRs belonging to G15 and G17 (likely infecting microbes from the "Pelagibacterales" clade) have a similar coverage profile, with the islands representing around 10% of the phage genomes and a maximum island length of 1.5 kb (**Figures 1B,D**). A very different pattern with much more uneven recruitment was found for the G16 CGR representative with islands covering around 20% of the phage genome. Besides, the 3.7 kb largest island in this genome had almost no coverage (**Figure 1C**), suggesting that the whole fragment is absent in the metavirome.

This presence of under-recruiting islands in phage genomes is reminiscent to what was discovered years ago for cellular genomes (both bacteria and archaea) in which similar regions were found even when the metagenomes were from the same habitat from which the strains providing the genomes were isolated (Cuadros-Orellana et al., 2007). There are two possible explanations for this remarkable phenomenon. One would be that during the time span between the isolation of the strain and the sample taken to generate the metagenome the microbe has changed and most of the clones present in the metagenome are different in the underrecruiting region, in effect a Red-Queen (R-Q) arms race between phage and host accounting for the variability (Stern and Sorek, 2011). The other alternative is that prokaryotic species live in populations that are heterogeneous and diverse, containing several clonal lineages that share stretches of their genomes (more or less coincidental with their core genomes), and other regions that vary from one clone to another, and that correspond to the flexible genome. Moreover, there is evidence indicating that clones survive without significant changes for periods of decades in natural environments (López-Pérez et al., 2013). In addition, although there are fewer available examples to examine, even when the strain genome and the metagenome are from the same sample the under-recruiting islands are still detectable (Gonzaga et al., 2012). Some of us proposed that to maintain this diversity of coexisting clonal lineages, a kill-the-winner dynamics has to occur, involving phages that equalize the prokaryotic populations, preventing any clone from sweeping the others out of the population (C-D) (Rodriguez-Valera et al., 2009). To carry out this role, phage population would have to be also polyclonal with multiple concurrent lineages. This would also explain the situation depicted by **Figure 1**.

We also found 11 CGRs (ca. 20% of the abundant ones) in which not a single island could be identified (**Figures 1A,E** and Data S2). All except one of the MVI free CGRs belong to G19, a group that was described as putative pelagiphage based on similarity to a prophage detected in the genome of the alpha proteobacterium HIMB114, a member of the proposed order "Pelagibacterales" obtained from the coastal tropical North Pacific (Grote et al., 2012). Along similar lines, we recently found a nearly identical genome fragment in the Mediterranean (ca. 40 kb) with *>*97% identity to the cyanophage isolate S-CAM1 from the Pacific Ocean (South California Coast) which indicates that some phage clones are very stable and durable (Mizuno et al., 2013). These phage clones appear to be conserved and present in significant amounts in two independent samples separated by more than 4 years. However, even in this case, the phage population was still composed of different clones as illustrated by the sequence differences detected among individual genomes and by the uneven recruitment. It is possible that these are generalistic phages that recognize multiple host clonal lineages (Flores et al., 2011).

## **GENOMIC ISLANDS IN PHAGES AND MVIs**

Like in prokaryotic genomes, phage genomes are composed of more conserved regions that could be called "core" and regions that vary ("flexible") among otherwise closely related genomes (Angly et al., 2009). In a previous work (Mizuno et al., 2013), comparing the same sets of CGRs we detected that these genomes display a somewhat uneven similarity in which nearly identical genomic regions are juxtaposed with regions with no sequence similarity whatsoever, that is, with essentially completely different genes. Most phage genomes, even from the same location and sample, appear to be extremely plastic with large flexible genomic islands that vary in sequence from one clone to another within a framework of relatively similar genomes (Mizuno et al., 2013). This is similar to what happens in prokaryotic genomes where the flexible genome (variable from one strain to another) tends to concentrate in genomic islands. In prokaryotic cells, the flexible genomic islands tend to be also under-recruiting in metagenomes. In other words, they are also MGIs. This seems to be the case with phages and genomic islands (GIs), as shown in **Figure 2A** where we have compared the pelagiphage genome HTVC010P retrieved from Bermuda Hydrostation S (Sargasso Sea) (Zhao et al., 2013) with a cluster (C28) of three related MedDCM pelagiphage genomes and their recruitments from the metavirome. The regions that appear more specific to a particular phage lineage, reminiscent of flexible genomic islands of prokaryotes,

highlighted in gray. Vertical bars indicate RPKG values (left Y-axis, reads recruited per Kb of sequence per Gb of metagenome). Median coverage of each CGR is shown as a red line (right Y-axis). Host prediction for each CGR, wherever available, is indicated by a colored circle on top of each bar (orange: pelagiphage, green: cyanophage,

S42-C7 **(B–E)** Recruitment plots showing coverage of selected CGRs in the metavirome. These CGRs are also labeled and highlighted in the same color in the top panel. The horizontal red dotted line indicates the median base coverage and MVIs are shown as red rectangles in each plot.

also recruited at lower levels than the more conserved tracts of these genomes. Although this might be considered an obvious expectation, it reflects a persistence of the situation that is remarkable. It shows that the conserved regions among these closely related lineages transcend large geographic distances (Bermuda to the Mediterranean) and significant time gaps (4 years).

Another interesting example can be seen in **Figure 2B**, where both phage genomes in cluster C10 share around 60% of the genome at a high sequence identity (more than 90%). However, there are two regions identified in **Figure 2B** as GI-1 and GI-2 with little or no sequence similarity. GI-1 in S30-C55 recruits similarly to the rest of the genome, but the corresponding region in S46-C61 does not. In the same manner, the S46-C61 version of GI-2 recruits (albeit unevenly) while the corresponding region in S30-C55 does not recruit at all. This illustrates that some phage genomic regions are extremely labile. In this case GI-2 that codes (among other proteins) for the tail fiber and the carbohydrate binding protein (CBP) of the S30-C55 phage genome appears to have disappeared from the population; either as a result of the lineage that contained this island becoming extinct during the 4 years elapsed between the retrieval of the genomes and the metavirome, or it is present at extremely low abundance and is not detected at this level of sequencing. To a lesser degree, the same happens with GI-1 of S46-C61.

In conclusion, the results indicate clearly how the most variable regions when comparing closely related genomes are also the ones under-recruiting, similarly to what happens with cellular genomes. The patterns found for the phage GIs are very similar to those found to affect the O-chain of the LPS in Gramnegative marine bacteria (Rodriguez-Valera et al., 2009; Hooton et al., 2011) that is supposed to be a major phage receptor (Avrani et al., 2011; Hooton et al., 2011). However, a major difference is that while the gene clusters in the cellular genomes are extremely divergent, even at the level of gene content and size of the clusters, in the phages, the size and number of genes in the MVIs were quite well preserved. That is, while the cellular genomes MGIs contain different gene clusters, MVIs seemingly contain the same genes but with very different sequence.

## **GENES FOUND IN MVIs**

We found a total of 390 predicted genes in CGRs MVIs. The vast majority (252) could not be assigned any function (hypothetical proteins) (**Table 1**, Data S1). The most frequent identifiable genes in MVIs were annotated as containing domains involved in molecular recognition. Specifically, these were the putative carbohydrate binding domain containing protein (14), the tail fiber protein (12), C1q globular head like domain containing protein (10), and lectin like domain protein (9). The tail fiber protein is known to be involved in host recognition (Yu and Mizushima, 1982; Bartual et al., 2010; Garcia-Doval and Van Raaij, 2012) and many studies have reported the high variability of this protein between nearly identical phages (Angly et al., 2009; Labrie et al., 2013).

In phages there are two stages in identifying the cells to be infected (Rakhuba et al., 2010). The initial step is reversible and is carried out by the tail fiber. The second leads to irreversible binding to the host cell and is carried out by the receptor binding protein (RBP) that is highly variable depending on the phage structure and classification. The genes coding for both steps are

**Table 1 | Annotation and frequency of the most frequent genes found in MVIs.**


typically located close to each other in phage genomes, what facilitates the identification of this recognition module (Miller et al., 2003; Leiman et al., 2004). From the most frequent underrecruiting genes, we identified three as the most likely candidates to be the RBP due their well-known capability to interact with other proteins or other substrates. These are proteins containing C1q, carbohydrate binding or a lectin like domains. In addition, fibrinogen-like coiled coil and filamentous hemagglutinin-like proteins were also identified as such, considering the structural similarity with the viral receptor-recognition proteins of Influenza virus (Skehel and Wiley, 2000). In summary, out of 138 identifiable genes in MVIs, 59 (**Table 1**) are clearly related to the host recognition module and therefore could reflect the specialization in infecting different clones of host, and they are enriched in MVIs as the C-D model would predict.

Another conspicuous component of under-recruiting islands (**Table 1**) were the small and large subunits of the terminase (these proteins are involved in packaging of the phage genome into the capsid) and other genes whose involvement in host-phage interaction appears more remote (if any). For example, three internal virion proteins, a protein involved in host cell wall penetration appeared also in MVIs. The reasons for the under-recruitment of these genes are more obscure. The repeated identification of these genes within MVIs in different phage genomes makes it unlikely that their under-recruitment is artifactual and probably reflects our incomplete understanding of phage biology.

## **METAGENOMIC RECRUITMENT WITHIN THE HOST RECOGNITION MODULE GENES**

A way to glimpse at the biological meaning of the differences in recruitment throughout phage genomes is to dissect the recruitment of different regions, with well-established roles, within a single gene. Regions that under-recruit within a single gene are either very variable among the different lineages in the population (C-D) or within the same lineage at different points in time (R-Q). Since the pelagiphages are the highest recruiting phages in the metavirome (**Figure 1**), we have chosen them for more detailed recruitment analyses. Both cluster C28 and C24 contain CGRs that are predicted to prey on *Ca.* Pelagibacter due their relatedness to the pelagiphage HTVC010P (Mizuno et al., 2013). We have selected those genes in these genomes (clusters C28 and C24) that could be identified as involved in both host recognition steps of the phage life cycle (see above). Specifically, these are the tail fiber (the first step) and C1q domain containing genes (RBPs, the second step). We examined the recruitment patterns of the different functional protein domains identifiable within these genes.

In order to identify the most variable region within the tail fiber protein, we have used the tail fiber of the phage T7 of *Escherichia coli* as a reference (Garcia-Doval and Van Raaij, 2012). This protein (553 aa) shown in the top panel of **Figure 3** has a very characteristic domain structure, including a well conserved N-terminal domain (PF03906), followed by several repeats. The C-terminal domain (gp37\_C) is known to be involved in host interaction in the T7 bacteriophage. The crystal structure of the C-terminal domain has been described (**Figure 3A**) and a number of possible binding receptors suggested, leaving little doubt that

this region is involved in host recognition by T7 (Garcia-Doval and Van Raaij, 2012). The tail fiber protein sequences from the selected phage genomes of cluster C28 and C24 preserved the same domain structure (T7-like) at the N terminus (PF03096). This was followed by a repeat-pattern specific to each tail fiber protein but without any sequence similarity. A similar repetitive stretch is also seen in the T7 tail fiber protein, suggesting a similar secondary structure (**Figure 3**). No similarity was discernible in the C-terminal regions. One of the tail fiber proteins has a C-terminal YadA domain (**Figure 3C**) that is a typical adhesion like domain (Casutt-Meyer et al., 2010; Edwards et al., 2010). These observations taken together suggested that the C-terminal region of the tail fiber of these pelagiphages are much more variable than the N terminus and are likely candidates to be involved in host recognition, as shown for the Enterobacteria phage T7. Along these lines, the metagenomic coverage plot of both tail fibers shows a marked decrease for the C terminal region, particularly in the case of the representative of cluster C24, the region where the YadA domain is located.

Even more compelling is the example of the putative RBP in these phages (e.g., the C1q domain containing gene) found next or close to the tail fiber. The C1q protein is a wellknown target recognition molecule of the classical vertebrate serum complement pathway (Kishore et al., 2004; Ghai et al., 2007). **Figure 4** shows the recruitment levels of two C1q domain containing proteins from pelagiphages of both C28 and C24. These domains always under-recruited conspicuously, fitting well with the expectation of a highly variable protein that is likely involved in host recognition. In both the examples shown, the N-terminal end recruits far more than the C1q domain. We created homology-based protein structure models for both these C1q domains and mapped the recruitment from the metagenome to the protein structures. Both these C1q domains are only 38% similar to each other, and display different regions of variability (blue in the right-hand panels of **Figure 4**). It appears that these regions, that under-recruit in comparison to the rest of the sequence may be the sites of hyper-variability in these proteins. These results indicate several levels of variability even among concurrent phage lineages, ranging from large MVIs specific to each clonal lineage, and likely also existence of numerous sublineages even amongst each, as can be observed by the variation in the structural scaffold of C1q. Phage genomes are naturally constrained by a fixed size and show much more differences in specific sequence regions within genes rather than different genes altogether. This could just reflect the obvious constraints of genome size affecting phages and viruses in general.

## **CONCLUSIONS, RED QUEEN OR CONSTANT-DIVERSITY**

The identification of genomic islands in these genomes (metaviromic and flexible) and their differential persistence across time and space leads us to speculate on the nature of the dynamics in these populations. While R-Q predicts that the molecular modules involved in host-phage interaction will vary very rapidly with evolutionary time, C-D posits that many versions of such modules are present simultaneously, at any given time, within the population. The differences between these two models are subtle and linked to time. In our datasets, the time span between genome and metavirome retrieval is only 4 years, which, evolutionarily speaking, is a short period in this relatively stable habitat. There is also evidence that some host-recognition modules (e.g., tail fiber protein) are preserved over large geographic distances (S31- C64 from MedDCM and HTVC010P from Bermuda, **Figure 2**). Moreover, within this short time span, it was possible to retrieve phage genomes that do not present any MVIs (i.e., persistent lineages). However, we have also found cases where a MVI did not recruit any reads at all (for example, GI-2 in S30-C55, **Figure 2**), suggesting that this particular lineage might have disappeared altogether (or is undetectable at current sequencing depth). Even so, taken together these results appear to support persistence in lieu of continuous variability, and hence favor more a C-D, than a R-Q dynamics in these populations.

It seems clear that the genomic diversity of phage populations is outstanding and yet, some genomes appear very well preserved over long distances and time spans. Overall, the data support the C-D model (Rodriguez-Valera et al., 2009) and a more recent reformulation (Rodriguez-Valera and Ussery, 2012) in which some of us proposed that prokaryotic species and their accompanying phages form consortia that have been evolving together over many millions of years. Akin to the starters used in dairy, such consortia work well as a package that have been selected for in nature and can be found with very similar structure at similar habitats worldwide. If this hypothesis is true, the consortium that develops in the deep (below 20 m) photic zone of stratified waters such as the MedDCM could be one of the most extensive on Earth and of enormous relevance to the global ecology of the planet.

## **AUTHOR CONTRIBUTIONS**

Francisco Rodriguez-Valera conceived the work. Carolina Megumi Mizuno and Rohit Ghai performed all analyses. The manuscript was written by Francisco Rodriguez-Valera and Rohit Ghai.

## **ACKNOWLEDGMENTS**

This work was supported by projects MAGYK (BIO2008- 02444), MICROGEN (Programa CONSOLIDER-INGENIO 2010 CDS2009-00006), CGL2009-12651-C02-01 from the Spanish Ministerio de Ciencia e Innovación, DIMEGEN (PROMETEO/2010/089) and ACOMP/2009/155 from the Generalitat Valenciana and MaCuMBA Project 311975 of the European Commission FP7. FEDER funds supported this project. Rohit Ghai was supported by a Juan de la Cierva scholarship from the Spanish Ministerio de Ciencia e Innovación.

## **SUPPLEMENTARY MATERIAL**

The Supplementary Material for this article can be found online at: http://www*.*frontiersin*.*org/journal/10*.*3389/fmicb*.*2014*.* 00027/abstract

## **REFERENCES**


populations of *Alteromonas macleodii*. *Genome Biol. Evol.* 4, 1360–1374. doi: 10.1093/gbe/evs112


**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.

*Received: 08 November 2013; accepted: 16 January 2014; published online: 03 February 2014.*

*Citation: Mizuno CM, Ghai R and Rodriguez-Valera F (2014) Evidence for metaviromic islands in marine phages. Front. Microbiol. 5:27. doi: 10.3389/fmicb.2014.00027 This article was submitted to Evolutionary and Genomic Microbiology, a section of the journal Frontiers in Microbiology.*

*Copyright © 2014 Mizuno, Ghai and Rodriguez-Valera. 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.*

## The agricultural antibiotic carbadox induces phage-mediated gene transfer in *Salmonella*

#### *Bradley L. Bearson1 \*†, Heather K. Allen 2†, Brian W. Brunelle2, In Soo Lee1,3, Sherwood R. Casjens <sup>4</sup> and Thaddeus B. Stanton2*

*<sup>1</sup> Agroecosystems Management Research Unit, National Laboratory for Agriculture and the Environment, ARS, USDA, Ames, IA, USA*

*<sup>4</sup> Division of Microbiology and Immunology, Department of Pathology, University of Utah, Salt Lake City, UT, USA*

#### *Edited by:*

*Stephen T. Abedon, The Ohio State University, USA*

#### *Reviewed by:*

*Todd R. Callaway, Agricultural Research Service, USA Paul Ebner, Purdue University, USA*

#### *\*Correspondence:*

*Bradley L. Bearson, Agroecosystems Management Research Unit, National Laboratory for Agriculture and the Environment, ARS, USDA, 2110 University Drive, NSRIC-2103, Ames, IA 50011, USA e-mail: brad.bearson@ars.usda.gov*

*†These authors have contributed equally to this work.*

Antibiotics are used for disease therapeutic or preventative effects in humans and animals, as well as for enhanced feed conversion efficiency in livestock. Antibiotics can also cause undesirable effects in microbial populations, including selection for antibiotic resistance, enhanced pathogen invasion, and stimulation of horizontal gene transfer. Carbadox is a veterinary antibiotic used in the US during the starter phase of swine production for improved feed efficiency and control of swine dysentery and bacterial swine enteritis. Carbadox has been shown *in vitro* to induce phage-encoded Shiga toxin in Shiga toxin-producing *Escherichia coli* (STEC) and a phage-like element transferring antibiotic resistance genes in *Brachyspira hyodysenteriae*, but the effect of carbadox on prophages in other bacteria is unknown. This study examined carbadox exposure on prophage induction and genetic transfer in *Salmonella enterica* serovar Typhimurium, a human foodborne pathogen that frequently colonizes swine without causing disease. *S.* Typhimurium LT2 exposed to carbadox induced prophage production, resulting in bacterial cell lysis and release of virions that were visible by electron microscopy. Carbadox induction of phage-mediated gene transfer was confirmed by monitoring the transduction of a *sodCIII*::*neo* cassette in the Fels-1 prophage from LT2 to a recipient *Salmonella* strain. Furthermore, carbadox frequently induced generalized transducing phages in multidrug-resistant phage type DT104 and DT120 isolates, resulting in the transfer of chromosomal and plasmid DNA that included antibiotic resistance genes. Our research indicates that exposure of *Salmonella* to carbadox induces prophages that can transfer virulence and antibiotic resistance genes to susceptible bacterial hosts. Carbadox-induced, phage-mediated gene transfer could serve as a contributing factor in bacterial evolution during animal production, with prophages being a reservoir for bacterial fitness genes in the environment.

#### **Keywords:** *Salmonella***, bacteriophage, antibiotic, carbadox, prophage, gene transfer, transduction**

## **INTRODUCTION**

Antibiotics are used to treat bacterial infections, prevent bacterial infections, or improve feed conversion efficiency in foodproducing animals. However, antibiotics have broad and unintended, sometimes called collateral, effects on microorganisms and the microbial communities they inhabit. A microorganism's response to antibiotic exposure can be monitored by gene expression signatures that indicate the organism's physiological response to antibiotic-related stress (Brazas and Hancock, 2005). Sub-minimal inhibitory concentrations of antibiotics in particular have been shown to have the unintended effect of modulating gene expression in various bacteria (Davies et al., 2006). Transcription in the foodborne pathogen *Salmonella enterica* serovar Typhimurium is affected by sub-inhibitory concentrations of antibiotics of agricultural importance, such as quinolones (Yim et al., 2011) and tetracyclines (Brunelle et al., 2013). In both studies virulence genes were among those upregulated by these antibiotics, suggesting that low concentrations of certain antibiotics may promote rather than inhibit *S*. Typhimurium survival in the host.

In the US, salmonellae are the leading cause of bacterial foodborne morbidity and mortality for humans (Scallan et al., 2011). The prevalence of multidrug-resistant *Salmonella* isolates has increased over the last few decades, and outbreak investigations indicate that antimicrobial resistant *Salmonella* isolates are associated with an increased rate of hospitalization (Varma et al., 2005a). Furthermore, patients infected with antimicrobial resistant *Salmonella* have an increased frequency of bloodstream infections and longer lengths of hospitalization (Varma et al., 2005b). Acquisition and carriage of antibiotic resistance genes by *Salmonella* is therefore a critical factor in the degree of human morbidity and mortality caused by this pathogen.

*Salmonella* has acquired antibiotic resistance genes from the environment. Its primary habitat is within the microbial community of the intestinal tract, and this community, or gut microbiota, is a reservoir for antibiotic resistance genes (Salyers

*<sup>2</sup> Food Safety and Enteric Pathogens Research Unit, National Animal Disease Center, ARS, USDA, Ames, IA, USA*

*<sup>3</sup> Department of Biological Sciences and Biotechnology, Hannam University, Daejeon, South Korea*

et al., 2004). Sub-inhibitory antibiotics promote resistance gene transfer among gut bacteria via transposons (Shoemaker et al., 2001; Song et al., 2009), plasmids (Feld et al., 2008), and phagelike gene transfer agents (GTAs) (Stanton et al., 2008). The agricultural antibiotic carbadox is frequently used in the US during the starter phase of swine production for performance enhancement and control of enteric diseases. Carbadox is an antibacterial agent used exclusively in animals. For growth promotion and disease prophylaxis, swine feed contains 10–25 g/ton [11– 28 mg/kg or parts-per-million (ppm)] and 50 g/ton [55 mg/kg (ppm)], respectively. Carbadox, a quinoxaline-di-*N*-oxide, is mutagenic, causing base pair substitutions and frameshift mutations in DNA (Beutin et al., 1981). A range of carbadox concentrations from 0.5 to 8μg/ml (ppm) has been shown to induce prophages in Shiga toxin-producing *Escherichia coli* (STEC) (Kohler et al., 2000) and GTAs in *Brachyspira hyodysenteriae* (Stanton et al., 2008). However, it is unknown what effect carbadox would have on prophages encoded by other bacterial species, including those native to *S*. Typhimurium strains.

*Salmonella* strains have multiple prophage genomes integrated into their chromosomes. For example, the genome of *S.* Typhimurium strain LT2 contains four functional prophages: Gifsy-1 and -2 and Fels-1 and -2 (McClelland et al., 2001; Casjens, 2011). Investigation of prophages in *S.* Typhimurium indicates that many of these prophages can be induced to produce infectious virions by various environmental signals including DNA damage, antibiotics such as mitomycin C, and hydrogen peroxide (Schicklmaier et al., 1998; Figueroa-Bossi and Bossi, 1999; Schmieger and Schicklmaier, 1999; Frye et al., 2005; Garcia-Russell et al., 2009). Furthermore, prophages often encode virulence genes that enhance the pathogenesis of the bacterial strain into which the prophage is integrated (Groman, 1955; O'Brien et al., 1984; Cheetham and Katz, 1995; Waldor and Mekalanos, 1996; Figueroa-Bossi and Bossi, 1999; Mirold et al., 1999; Figueroa-Bossi et al., 2001; Ho et al., 2002; Casjens and Hendrix, 2005).

Since *Salmonella* strains usually contain multiple functional prophages and frequently colonize the swine intestinal tract, the goal of the current study was to evaluate prophage induction and genetic transfer in *S*. Typhimurium following carbadox exposure. Our research demonstrates that carbadox induced prophage production, thereby generating infectious virions capable of transferring virulence and antibiotic resistance genes via a prophage genome or generalized transduction.

## **METHODS**

## **BACTERIAL STRAINS, MEDIA, AND CHEMICALS**

Bacterial strains (**Table 1**) were grown in LB (Lennox Laboratory Supplies, Dublin, Ireland) or E minimal medium containing 0.4% glucose (Vogel and Bonner, 1956). A 5 mg/ml carbadox (Sigma-Aldrich, St. Louis, MO, USA) stock solution was made in 0.1 N NaOH and used at a final concentration of 2.5 μg/ml unless noted otherwise. Other antibiotics were used at the following concentrations: ampicillin (100μg/ml), kanamycin (50μg/ml), tetracycline (20μg/ml), chloramphenicol (30μg/ml), and carbenicillin (50μg/ml).

## *S***. TYPHIMURIUM GENE AND PROPHAGE KNOCKOUTS BY RECOMBINEERING**

Oligonucleotide primers for PCR amplification and construction of gene and prophage knockouts are listed in **Table 2**. *S.* Typhimurium gene and prophage knockouts were constructed by recombineering (recombination-mediated genetic engineering) as previously described (Bearson and Bearson, 2008; Bearson et al., 2008). Briefly, the 5 end of a gene knockout primer (bold, **Table 2**) has homology to 32–44 bp of the target gene whereas the 3 end contains universal sequences (underlined) to amplify an antibiotic resistance gene and truncate potential translation of the target gene. A gene knockout primer set was used to PCR amplify either the *neo* or the *cat* gene. Gel electrophoresis was performed on the amplification product of a knockout fragment and the respective DNA fragment was gel extracted using a Freeze'n Squeeze column (Bio-Rad, Hercules, CA). Each knockout fragment was transformed (Sambrook and Russell, 2006) into an arabinose-induced *S.* Typhimurium strain containing the pKD46 plasmid (Datsenko and Wanner, 2000). Transformants containing the knockout were selected on LB agar medium containing kanamycin. If necessary, the gene knockout with the *neo* marker was moved to another strain background by transduction using a P22 phage with a high transduction frequency. Flp mediated deletion of the *neo* or *cat* gene was performed by transferring the pCP20 plasmid into the knockout strain by either transduction or transformation followed by a procedure to screen for loss of resistance to kanamycin (Cherepanov and Wackernagel, 1995).

## **DETERMINATION OF PHAGE TITERS FOLLOWING CARBADOX INDUCTION OF BACTERIAL STRAINS LYSOGENIZED WITH THE MODEL PHAGES P22, λ, HK97, AND SF6**

Bacterial strains lysogenized with P22 (UB-1790), λ (UB-1703), HK97 (UB-1704), and Sf6 (UB-1496) were grown in LB broth at <sup>37</sup>◦C with shaking. At a density of 1 <sup>×</sup> <sup>10</sup><sup>8</sup> bacterial cells/ml, carbadox was added to a final concentration of 0.5 μg/ml (ppm). At the indicated times after carbadox addition, aliquots of the cultures were removed, shaken with several drops of chloroform to complete lysis, and titered on the permissive host. Strains used to titer the phage lysates were DB7004 (P22), 594 (λ and HK97), and UB-1458 (Sf6).

## **CARBADOX INDUCTION OF WILD-TYPE** *S***. TYPHIMURIUM**

*S.* Typhimurium strains were grown in LB 0.4% glucose at 37◦C with shaking at 180 rpm. Carbadox was added to cultures at OD600 = 0.2 at a final concentration of 2.5μg/ml (ppm). Incubation of cultures was continued to monitor for bacterial cellular lysis.

## **PHAGE TRANSDUCTION USING CARBADOX-INDUCED** *S***. TYPHIMURIUM LYSATES**

Supernatants from non-induced and carbadox-induced bacterial cultures were harvested by adding 100μl of chloroform, vortexing gently, and allowing the culture to incubate for ∼15 min with shaking. The cultures were centrifuged at 1000× g for 20 min and the supernatant was transferred to a fresh tube containing 200μl of chloroform for storage. Bacterial lysates were plated to LB to ensure that viable bacterial cells were not present. An

#### **Table 1 | Bacterial strain list.**


*\*Known antibiotic resistance phenotypes. Ap, Ampicillin; Cam, Chloramphenicol; Tet Tetracycline; Str, Streptomycin; Su, Sulfamethoxazole; Sp, Spectinomycin; Kn, Kanamycin; Nal, Nalidixic acid.*

#### **Table 2 | Primers used in this study.**


*Bold text indicates homology to the target gene/phage.*

*Underlined text indicates sequence for neo or cat amplification.*

overnight culture of the transduction recipient was grown in LB or LB glucose at 37◦C with shaking. The transduction was performed with equal volumes of both the *S.* Typhimurium recipient strain and either the non-induced or the carbadox-induced culture supernatant. The transduction was incubated at 37◦C for ∼1–3 h before plating on the appropriate selective medium. To determine the transduction frequency for transfer of the histidine operon, the transduction was plated onto E glucose minimal medium. Transduction frequency for antibiotic resistance gene transfer (*sodCIII*:: neo in Fels-1, plasmid-encoded kanamycin, and chromosomally encoded tetracycline) was determined by plating the transduction to LB containing either kanamycin or tetracycline.

### **PHAGE PURIFICATION FOR ELECTRON MICROSCOPY**

Overnight cultures were diluted 1:100 in 400 ml E glucose minimal medium and incubated at 37◦C with shaking. At OD600 = 0.5, carbadox was added to a final concentration of 2.5μg/ml, and incubation continued until lysis was achieved. Phages were purified and visualized by electron microscopy as described previously (Humphrey et al., 1997) with the following modifications. Purified phages were negatively stained by mixing phage samples with an equal volume of phosphotungstic acid (2%, pH 7.0). Samples were deposited on Formvar-coated 200-mesh carbonreinforced copper grids (Electron Microscopy Sciences, Hatfield, PA) and viewed with a FEI Tecnai G2 BioTWIN electron microscope (80 kV; Hillsboro, OR).

## **RESULTS AND DISCUSSION**

### **CARBADOX INDUCES** *S***. TYPHIMURIUM AND OTHER** *Enterobacteriaceae* **PROPHAGES TO CAUSE PHAGE REPLICATION AND CELL LYSIS**

To test whether well-characterized prophages in *S. enterica* and other *Enterobacteriaceae* species are induced by carbadox, we monitored the infectious phages produced by carbadox-treated cultures of bacterial strains that were lysogenized by the "model system" phages P22 (*S. enterica*), λ and HK97 (*E. coli*), and Sf6 (*Shigella flexneri*). In each case, at least partial clearing of the culture indicated a substantial fraction of the cells in the culture had lysed by about 2 h. All four prophages gave an approximately three-log increase in free phage in the culture (**Figure 1**), with

**FIGURE 1 | Carbadox induction of** *Enterobacteriaeae* **prophages.** Bacterial strains lysogenized with P22 (UB-1790), λ (UB-1703), HK97 (UB-1704), and Sf6 (UB-1496) were grown in LB broth at 37◦C with shaking. At a density of 1 <sup>×</sup> 108 bacterial cells/ml, carbadox was added to a final concentration of 0.5μg/ml. Phage lysates were obtained by shaking with chloroform at the indicated times after carbadox addition and titered on a permissive host. Open symbols indicate cultures with no added carbadox and closed symbols indicate cultures with carbadox added. The different induced lysogens are indicated in the figure insert.

about 10–200 progeny phage produced per bacterium that were initially present at the time of carbadox addition. Several concentrations of carbadox were tested, and 0.5 μg/ml (ppm) is shown as the minimum that gave good induction of P22. Carbadox induces the P22 *Salmonella* prophage as well as *E. coli* and *S. flexneri* prophages under these conditions, and this is not surprising given its apparent mechanistic similarities to the action of mitomycin C. Carbadox is a DNA damaging agent, and DNA damage induced by mitomycin C is known to induce the bacterial SOS pathway, which induces prophages. It is likely capable of induction of prophages from many if not all *Enterobacteriaceae* bacterial species, as well as more distantly related bacterial phyla.

To examine prophage induction by carbadox in wild-type *S.* Typhimurium isolates containing their natural prophages, we initially examined *S.* Typhimurium LT2, a strain that is widely used in the study of *Salmonella* genetics in the laboratory (Lilleengen, 1948). Cultures of wild-type LT2 in early log phase growth were exposed to carbadox. The bacterial density of the culture abruptly decreased at ∼2 h following exposure to 2.5μg/ml (ppm) of carbadox (**Figure 2**). Mitomycin C exposure is known to result in the induction of prophage Fels-1 from strain LT2 (Garcia-Russell et al., 2009). The fact that carbadox induces wild-type phage λ (above), which is only known to be induced by the SOS function of activated RecA protein (Little, 1993), strongly supports this mechanism for carbadox-mediated induction. Thus, the decrease in LT2 cell density is almost certainly due to bacterial cell lysis resulting from prophage induction.

LT2 is known to carry four prophages. Amplification of the integrase gene from DNA extracted from phage heads indicated that at least one of these, the Fels-1 prophage, was induced following carbadox exposure (Stanton and Humphrey, unpublished data). The BBS 1008 strain is an LT2 derivative from which all four prophages have been deleted. Following exposure of BBS 1008 to 2.5μg/ml of carbadox, the bacterial culture density did not decrease (**Figure 2**), indicating that these prophages are responsible for the bacterial cell lysis phenotype induced by carbadox exposure.

Examination of the purified phages by electron microscopy demonstrated the presence of mostly empty phage heads in the carbadox-induced culture of LT2 (data not shown). Treatment of *S.* Typhimurium LT2 and other strains with mitomycin C is known to induce Fels-1 plaque-forming phages with poor efficiency (Figueroa-Bossi and Bossi, 1999), so it is perhaps not surprising that whole phage particles were not seen. Nonetheless, phage heads were either greatly reduced or were not observed in the absence of carbadox for the LT2 strain and in the presence of carbadox for BBS 1008. These results confirm prophage induction in a natural wild-type *S.* Typhimurium isolate following carbadox exposure.

## **CARBADOX EXPOSURE OF** *S.* **TYPHIMURIUM LT2 PROMOTES PHAGE TRANSDUCTION INTO A SUSCEPTIBLE BACTERIAL HOST**

To monitor phage transduction frequency, a *neo* gene was inserted by recombineering into the putative virulence factor *sodCIII* on the Fels-1 prophage of strain LT2 to create strain *S.* Typhimurium BBS 565. The carbadox-induced phage lysate from BBS 565 (LT2 *sodCIII*::*neo*) efficiently transduced the kanamycin-sensitive strain BBS 243, as demonstrated by 8 <sup>×</sup> <sup>10</sup><sup>3</sup> kanamycin-resistant transductants per ml of lysate. In the absence of carbadox induction, BBS 243 remained susceptible to kanamycin following transduction with the control supernatant from BBS 565. This indicates that carbadox-induced prophages can carry genetic material from a donor into a recipient bacterial strain.

## **CARBADOX-INDUCTION OF MULTIDRUG-RESISTANT** *S***. TYPHIMURIUM DT104 RESULTS IN GENERALIZED TRANSDUCTION OF CHROMOSOMAL AND PLASMID DNA**

Generalized transduction involves the packaging of random host DNA (genomic or plasmid) into a bacteriophage particle and the transfer of that DNA to a recipient strain. Bacteriophage P22 is a generalized transducing phage that is commonly used for genetic experiments with *Salmonella* (Zinder and Lederberg, 1952; Kropinski et al., 2007). Although *S.* Typhimurium LT2 contains multiple prophages, it does not contain a P22-like prophage that performs generalized transduction. Some multidrug-resistant *S.* Typhimurium isolates like phage type DT104 harbor a P22 like prophage. This prophage has been described as PDT17, ST104, and prophage 1 (Schmieger and Schicklmaier, 1999; Tanaka et al., 2004; Cooke et al., 2008). To determine whether carbadox exposure could induce generalized transduction, *S.* Typhimurium DT104 NCTC13348 was exposed to carbadox. The culture lysed, and the phage-containing supernatant was harvested. The BBS 243 strain, a histidine auxotroph lacking the genes *hisDCBHA*, could be successfully transduced to *his*+ with this lysate. Generalized transduction was observed by growth on minimal medium lacking histidine, demonstrating the transfer of the *his* operon from the *his*+ donor (DT104) to the *his*− recipient (BBS 243) (**Table 3**). The results show that carbadox exposure of multidrug-resistant *S.* Typhimurium DT104 can stimulate generalized transduction and therefore chromosomal gene transfer.

The *Salmonella* genomic island-1 (SGI-1) is ∼43 kb and typically contains chromosomally-encoded resistance genes for multiple antibiotics including ampicillin, chloramphenicol, and tetracycline within an integron (Boyd et al., 2001; Mulvey et al., 2006). We attempted to transduce SGI-1 but were unsuccessful. Due to the size of SGI-1, transduction of this entire island by P22 into another *Salmonella* strain that does not already contain a segment of SGI-1 is inefficient [P22 packages about 43.4 Kbp of DNA (Kropinski et al., 2007)]. To overcome this experimental limitation, we utilized a DT104 derivative (BBS 1012) as our recipient strain. The BBS 1012 strain is a histidine auxotroph that has an internal deletion within SGI-1, resulting in sensitivity to ampicillin, chloramphenicol, and tetracycline due to the loss of multiple antibiotic resistance genes. In addition, BBS 1012 contains a natural plasmid that confers kanamycin resistance. Transduction of BBS 1012 with the carbadox-induced phage lysate from wild-type *S.* Typhimurium



DT104-530 (kanamycin sensitive) resulted in the transfer of tetracycline resistance following selection on tetracycline-containing LB agar medium. Transduction into the BBS 1012 strain was confirmed by growth on medium containing kanamycin and the absence of growth on minimal medium without histidine. Furthermore, the initial selection of BBS 1012 on LB medium containing tetracycline resulted in 100% co-transduction of resistance to both chloramphenicol and carbenicillin. The *floR* and *pse-1* genes encode resistance to chloramphenicol and carbenicillin, respectively, and these genes are located adjacent to *tetG* on SGI-1. These results indicate that exposing multidrug-resistant *S.* Typhimurium DT104 to carbadox can promote the transfer of numerous genes co-located within SGI-1 that encode resistance to multiple classes of antibiotics. Thus generalized transduction could participate in bacterial strain evolution by providing an assortment of antibiotic resistance genes for recombination within this important genomic region.

Although the SGI-1 integron in *S.* Typhimurium DT104 encodes multiple antibiotic resistance genes, some strains contain additional antibiotic resistance genes encoded on plasmids. Transduction of BBS 243 (kanamycin sensitive) with the carbadox-induced phage lysate from *S.* Typhimurium DT104 (745) resulted in bacterial growth on LB medium containing kanamycin, demonstrating the transfer of the plasmid encoding kanamycin resistance. Thus, carbadox exposure promoted generalized transduction of this natural plasmid as well as chromosomally-encoded antibiotic resistance genes in multidrug-resistant *S.* Typhimurium DT104.

## **CARBADOX-INDUCED GENE TRANSFER IS A GENERAL PHENOMENON THAT OCCURS IN MULTIDRUG-RESISTANT** *S***. TYPHIMURIUM STRAINS DT120 AND DT104**

Since the prevalence of multidrug-resistant *S.* Typhimurium strains has increased over the last few decades, we wanted to determine whether carbadox-induced gene transfer is unique to *S.* Typhimurium DT104 or is a general property of multidrugresistant *S.* Typhimurium strains. Phage lysates were harvested from several *S.* Typhimurium phage types following carbadox exposure and used to transduce the recipient strain BBS 243 (*hisDCBHA*) with selection on E glucose minimal medium; several different phage types were investigated as these, by definition (i.e., DT), should have varying prophage content. Carbadoxinduced phage lysates from several *S.* Typhimurium DT104 and DT120 isolates resulted in growth on minimal medium, indicating the transfer of the *his* operon to BBS 243 (**Table 3**). The results suggest that generalized transduction following carbadox induction is a common phenomenon for multidrug-resistant *S.* Typhimurium DT104 and DT120.

We PCR amplified the P22 gene *1* (encoding gp1/portal protein) from several of the DT104 and DT120 isolates that we have shown are capable of generalized transduction, suggesting that these multidrug-resistant isolates contain a P22 like prophage. To confirm that a phage capable of generalized transduction is required for carbadox-induced gene transfer, the P22-like prophage (prophage 1) was deleted from DT104 using recombineering. Gene transfer into BBS 243 was eliminated following transduction with a carbadox-induced phage lysate from the DT104 P22-like prophage knockout strain, indicating that the P22-like prophage is responsible for the generalized transduction from *S.* Typhimurium DT104. In support of this, the *S.* Typhimurium strains LT2, UK1, SL1344, and χ4232 (frequently investigated strains in the literature) are incapable of generalized transduction and do not contain a P22-like prophage. In contrast, genome scanning has identified integrated P22-like prophages in the genome sequences of isolates of *S. enterica* serovars Arizonae, Cholerae, Dublin, Hadar, Heidelberg, Houtenae, Johannesburg, Mississippi, Montevideo, Newport, Paratyphi, Rubislaw, Schwarzengrund, Tennessee, Typhimurium, Uganda, Wandsworth, and Welteverden; this suggests that P22 like prophages are common among *Salmonella* serovars and the potential for strain evolution due to generalized transduction is perhaps underappreciated. The capability of generalized transduction among *Salmonella* serovars is reinforced with the knowledge that P22 phage lysates are known to be stable for many years in the laboratory. Likewise, an ecological significance of phages in the environment is that DNA encapsulated within a phage head is protected from nucleases and therefore can survive outside of a bacterial host until encountering a recipient.

Swine environments, including swine manure, have been shown to contain abundant phage populations (McLaughlin et al., 2006; Wang et al., 2010). Bacteriophage populations present in manure could be derived principally from prophage induction of bacteria present in manure or a combination of induction from within the swine gastrointestinal tract and in manure. Prophage induction can be stimulated by various environmental signals and stresses including ultraviolet light, hydrogen peroxide, mitomycin C, and carbadox. Analysis of fecal phage metagenomes from medicated swine administered in-feed antibiotics [carbadox or ASP250 (chlortetracycline, sulfamethazine, and penicillin)] compared to non-medicated swine suggested that prophages were induced with antibiotic treatment (Allen et al., 2011). Similar work with mouse fecal phage metagenomes has shown that antibiotic treatment caused an increase in the abundance of phage-encoded antibiotic resistance genes (Modi et al., 2013). This suggests that antibiotic-induced phage-mediated transduction may contribute to antibiotic resistance gene transfer during animal production. Relatively little information is available concerning the extent of carbadox-induced prophage from bacteria, as *Salmonella* is only the third bacterial genus for which this response has been described (Kohler et al., 2000; Stanton et al., 2008). Additional information is needed to understand the capacity for carbadox to induce prophages during swine production since there is a potential for dissemination into the environment following manure application onto agricultural soils.

## **CONCLUSIONS**

Prophages are a potential environmental reservoir for bacterial fitness genes and may drive the emergence of new epidemic clones (Brussow et al., 2004). Prophages integrated in the genomes of *Salmonella* strains can encode genes associated with virulence or antimicrobial resistance (Figueroa-Bossi and Bossi, 1999; Figueroa-Bossi et al., 2001; Moreno Switt et al., 2013). Therefore, the pathogenic potential of a particular *Salmonella* strain depends in part upon the prophage repertoire integrated into the bacterial genome, and acquisition of prophages could conceivably result in enhanced bacterial virulence or survival during host colonization. In this report, we demonstrate that exposure of several different *S.* Typhimurium isolates to the agricultural antibiotic carbadox resulted in the production of transducing particles capable of transferring the individual phage genome as well as chromosomal and plasmid DNA by generalized transduction.

## **AUTHOR CONTRIBUTIONS**

Conceived and designed experiments: Bradley L. Bearson, Heather K. Allen, Sherwood R. Casjens, Thaddeus B. Stanton, and Brian W. Brunelle. Performed the experiments: Bradley L. Bearson, Heather K. Allen, Sherwood R. Casjens, In S. Lee, Brian W. Brunelle, and Thaddeus B. Stanton. Wrote and edited the manuscript: Bradley L. Bearson, Heather K. Allen, Sherwood R. Casjens, Brian W. Brunelle, and Thaddeus B. Stanton.

### **ACKNOWLEDGMENTS**

We are indebted to Stephanie Jones, Kellie Winter, Samuel Humphrey, and Briony Atkinson for outstanding technical assistance, to Judi Stasko for excellence with electron microscopy, to Kassandra Wilson and Eddie Gilcrease for measurement of phage titers, and to Kelly Hughes, Roger Hendrix, and Renato Morona for the gift of bacterial strains. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the US Department of Agriculture. USDA is an equal opportunity provider and employer.

#### **REFERENCES**


of *Salmonella enterica* serovar Typhimurium. *Microb. Pathog*. 44, 271–278. doi: 10.1016/j.micpath.2007.10.001


determinants in *Salmonella*. *Mol. Microbiol*. 39, 260–271. doi: 10.1046/j.1365- 2958.2001.02234.x


**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.

*Received: 24 September 2013; paper pending published: 02 November 2013; accepted: 23 January 2014; published online: 11 February 2014.*

*Citation: Bearson BL, Allen HK, Brunelle BW, Lee IS, Casjens SR and Stanton TB (2014) The agricultural antibiotic carbadox induces phage-mediated gene transfer in Salmonella. Front. Microbiol. 5:52. doi: 10.3389/fmicb.2014.00052*

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

*Copyright © 2014 Bearson, Allen, Brunelle, Lee, Casjens and Stanton. 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.*

## Exploring the risks of phage application in the environment

## *Sean Meaden\* and Britt Koskella*

College of Life and Environmental Sciences, University of Exeter, Penryn, UK

## *Edited by:*

Heather K. Allen, National Animal Disease Center, USA

#### *Reviewed by:*

Andres Moya, University of Valencia, Spain Natalia Ivanova, Lawrence Berkeley

### *\*Correspondence:*

National Laboratory, USA

Sean Meaden, College of Life and Environmental Sciences, University of Exeter, Penryn Campus, Treliever Road, Penryn, TR10 9EZ Penryn, UK e-mail: sm341@exeter.ac.uk

Interest in using bacteriophages to control the growth and spread of bacterial pathogens is being revived in the wake of widespread antibiotic resistance. However, little is known about the ecological effects that high concentrations of phages in the environment might have on natural microbial communities.We review the current evidence suggesting phagemediated environmental perturbation, with a focus on agricultural examples, and describe the potential implications for human health and agriculture. Specifically, we examine the known and potential consequences of phage application in certain agricultural practices, discuss the risks of evolved bacterial resistance to phages, and question whether the future of phage therapy will emulate that of antibiotic treatment in terms of widespread resistance. Finally, we propose some basic precautions that could preclude such phenomena and highlight existing methods for tracking bacterial resistance to phage therapeutic agents.

**Keywords: phage therapy, antibiotic resistance, coevolution, phage resistance, microbial communities**

"fmicb-04-00358" — 2013/11/28 — 20:30 — page 1 — #1

## **INTRODUCTION**

The selection for and subsequent evolution of antibiotic resistance in bacterial populations, both in the wider environment and during clinical treatment, presents a serious challenge to human health (Kåhrström, 2013). Although still controversial, it is increasingly clear that agricultural use of antibiotics has played a role in the continued selection for resistance genes and that the movement of these genes into pathogens of clinical relevance is possible (see Smith et al., 2009; van Cleef et al., 2010; Heuer et al., 2011; Zhu et al., 2013). Nosocomial environments also act as significant sources of antibiotic resistance and the transfer of antibiotic resistance genes into the agricultural environment has recently been demonstrated, suggesting the genetic exchange between environments works both ways (Price et al., 2012). As a result, alternative antimicrobial strategies are being sought. One such strategy is to utilize naturally occurring viral predators of bacteria: bacteriophages (phages). Lytic phages are capable of killing bacteria by invading and propagating within the host cell and then lysing open the cell to "burst" out, thus killing the bacteria. This is in contrast to temperate phages, which integrate into the genome of their hosts and can be transmitted vertically, serving as a refuge for phages in harsh environments (Svircev et al., 2011). The latter, although important to the ecology and evolution of bacterial populations, are not commonly considered for use as biocontrol and therefore will not be covered further in this review (but see Hyman and Abedon, 2010 for review of the effects of lysogeny on bacterial resistance). Despite being discovered as potential therapeutic agents over 80 years ago (d'Herelle, 1929), and their continual use in Russia and Georgia ever since (Kutter et al., 2010), few clinical trials of so-called "phage therapy" have been conducted in Western medicine (Wright et al., 2009; Sarker et al., 2012). Accordingly, no clinical phage therapy products are currently available in the West, and regulatory burdens may have dampened pharmaceutical interest, as years of research and clinical trials can cost millions of euros, presenting a formidable hurdle (Pirnay et al.,2011; Brüssow, 2012).

A more viable route to market has been provided by use of phages in agriculture and aquaculture (Jones et al., 2012; Martínez-Díaz and Hipólito-Morales, 2013). Recent review suggests that although phage therapy will not be the panacea that broad-spectrum antibiotics once were, phages could nonetheless play an important role in treating infections and maintaining food yields (Allen et al., 2013). Indeed, phage usage in agriculture has shown promise for treating numerous plant pathogens (Frampton et al., 2012), and some formulations of phages have been sold for large-scale environmental application (e.g., AgriPhageTM, Omnilytics). Likewise, the Food and Drug Administration in the USA has approved a product for the treatment of food products prior to market (ListshieldTM, Intralytix) and classed the use of phages in this specific context as "generally recognized as safe" (Food and Drug Administration (FDA), 2013). Another fruitful avenue for phage therapy may be aquaculture, an industry that has increased globally by over 10-fold in the last 30 years (Food and Agriculture Organization [FAO], 2012). Microbial diseases represent a severe threat to aquaculture productivity; accordingly, phages capable of lysing pathogens such as *Flavobacterium psychrophilum*, the causative agent of bacterial cold-water disease, have been isolated and tested as therapeutic agents (Kim et al., 2010). Combined, the scope for phage therapy to fill food production and clinical niches left vacant by redundant antibiotics seems vast. However, whilst phage therapy as a biopesticide could prove a useful tool, it also presents a risk for repeating the mistakes made with overuse of antibiotics, and the subsequently high levels of evolved antimicrobial resistance observed both in the environment and in hospitals (Zhu et al., 2013). For example, it is unclear whether introducing high concentrations of diverse phage types into the environment will select for broad resistance, making future treatments less likely to succeed. Furthermore, the role of phage-mediated selection in shaping bacterial growth rates and virulence to their hosts remains poorly understood.

Despite a number of reviews highlighting the need for increased understanding of environmental perturbations from anthropogenic antibiotic input (Martínez, 2008; Ding and He, 2010; Allen et al., 2013), the impact of antibiotic use on microbial communities has rarely been taken into account when designing treatment or application. This is particularly surprising given the known natural importance of these chemicals in shaping competition among bacterial strains (D'Costa et al., 2011). As phages are also known to select for resistant bacteria (Buckling and Rainey, 2002) and to mediate competition among strains (Bohannan and Lenski, 2000; Koskella et al., 2012), the same risks should apply to this alternative treatment. Unfortunately, little is currently known about the effects of applying high titers of phages to natural microbial communities. Most importantly, it is possible that with uncontrolled application of phages in the environment the future efficacy of phage therapy in a clinical setting could be reduced – a mistake we cannot afford given the need for new antimicrobial therapies as a result of drug-resistant pathogens (Levy and Marshall, 2004). Phage therapy in agriculture could serve as a testing ground for clinical use (Stone, 2002; Levin and Bull, 2004). However, there could also be a conflict of interest if cross-resistance to phage treatments is possible and if these resistant bacteria can spread from agricultural to clinical settings, as has been observed for antibiotic resistance (van Cleef et al., 2010; Perry and Wright, 2013). If phage therapy treatments fail, or improper use of phages in the environment goes unchecked, the use of widespread phage biocontrol in agriculture could jeopardize the future of phage therapy in hospitals. Fortunately, our understanding of phage-mediated selection is growing at a rapid pace and a new era of genomic investigation should allow monitoring of microbial communities following phage therapy. In this review, we will discuss the status of the field of phage therapy and consider the implications of phage host range and bacterial resistance. We suggest that with a few precautions phage therapy may be effective for treating bacterial infections in agriculture, aquaculture, healthcare, food production and food safety.

## **THE RISKS OF ANTIMICROBIAL USE IN AGRICULTURE**

The argument against using antibiotics as standard agricultural practice, both to improve growth rates and prevent disease, is not new (Witte, 1998) and has been extensively reviewed previously (Singer et al., 2003). However, unequivocally demonstrating increased resistance as a consequence of agricultural usage has proved elusive (Perry and Wright, 2013). A wave of new data supporting both direct and indirect routes of antibiotic resistance genes between agricultural and human populations suggests a bidirectional zoonotic exchange (Price et al., 2012). For example, recent studies have found diverse and abundant resistance genes in manure prior to disposal in the environment (Zhu et al., 2013) and a high prevalence of resistance to multiple antibiotics in enterobacteria isolated from tomato farms (Micallef et al., 2013) and in bacteria from manure-amended soils (Popowska et al., 2012). Furthermore, methicillin-resistant *Staphylococcus aureus* (MRSA) rates in workers on swine farms have been shown to be higher than for the average population in both North America and Europe (Voss et al., 2005; Khanna et al., 2008; Smith et al., 2009; van Cleef et al., 2010). Finally, calves treated with antibiotics are also more likely to carry MRSA and there is a direct association between intensity of animal contact and human MRSA carriage (Graveland et al., 2010). A similar trend is seen in aquaculture where bacteria nearer to farms were found to have higher levels of antibiotic resistance than nearby coastal regions in Italy (Labella et al., 2013). The increasing number of studies supporting the hypothesis that environmental use of antibiotics has contributed to selection for antibiotic resistance suggests that non-prudent use of antibiotics in healthcare and agriculture may reduce the effectiveness of antibiotic strategies as an essential treatment for disease.

As an alternative to antibiotic use, the application of phages in agriculture is being trialed as a biopesticide to control plant pathogens of tomato (Jones et al., 2012), citrus (Balogh et al., 2008), and onion (Lang et al., 2007) among others (reviewed in Svircev et al., 2011). For example, *Erwinia amylovora* (the causative agent of fire blight) infections are affecting a number of crop species in orchards across North America and Europe (see Malnoy et al., 2012 for review). Although antibiotics have traditionally been employed to control this disease, the emergence of streptomycin resistant strains (McManus et al., 2002) and a desire to reduce antibiotic use in the environment has led to the use of phages as an alternative. Phage biocontrol clearly has the potential to control fire blight infections, as lytic phages have been isolated that are highly infective to the pathogen, but definitive field trials are currently lacking. Given the evidenced risks of movement of antibiotic resistance genes between agricultural to human pathogens, we should ask whether the large-scale application of phages is likely to repeat these past mistakes. Until appropriate studies are conducted, the subsequent consequences of applying phages in agriculture for the spread of antibiotic resistance, the evolution of the pathogen, and the community of microbes within the plants and soil remain unknown.

## **DESIGN AND IMPLEMENTATION OF PHAGE THERAPY AND BIOCONTROL**

The process of preparing a phage therapy product for clinical use has been thoroughly described (Merabishvili et al., 2009; Gill and Hyman, 2010). **Figure 1** also describes this process for clinical and environmental samples. Briefly, environmental samples such as sewage or clinical samples from infected wounds are collected. The next step normally employs an "enrichment" process whereby the target bacterial species is added to the sample to increase the titer of phages infective to this strain. The sample is either filtered or chloroform is added to separate phages from bacteria, and individual phage "plaques" (i.e., the localized absence of bacterial growth in a lawn due to lysis) are chosen for further characterization. Transmission electron microscopy may be employed to assign family level phylogeny and genetic sequencing for finer scale taxonomic assignment, and screening of virulence factors is typically conducted. Other properties such as stability across a range of environmental conditions may be tested for optimal storage and production. Importantly, phage host range is typically tested to ensure the selected isolates have high efficacy against the pathogen of interest. However, this screening is most often done using a reference panel of laboratory stocks, rather than a large

"fmicb-04-00358" — 2013/11/28 — 20:30 — page 2 — #2

subset of bacteria from the local environment in which the phages will be applied, leading to a biased host range description. Therefore one way to reduce the possible community-level effects of applying phages would be to perform large-scale host range analyses across a biologically meaningful panel of isolates (i.e., those bacterial strains and species with which the phages are likely to interact once applied), as has been done successfully in the field of microbial ecology (Flores et al., 2011; Koskella and Meaden, 2013).

Once individual phages have been isolated and characterized, phage cocktails are produced by combining multiple, usually phylogenetically diverse, phages into one formulation. The idea behind these combined treatments is twofold: first, the use of multiple phages should increase the breadth of efficacy of the treatment to include most circulating strains of a pathogen; and second, the evolution of bacterial resistance should be slowed relative to single phage treatment. Whilst this approach could select for broadly resistant bacterial hosts, a number of studies suggest that broad resistance will carry a larger cost and therefore will not spread as rapidly (Bohannan and Lenski, 2000; Hall et al., 2012; Koskella et al., 2012). For example, strains of the plant pathogen, *Pseudomonas syringae*, evolved in the presence of three phage types were just as likely to evolve resistance against all three as those strains evolved in the presence of a single phage type. Furthermore, bacteria treated with multiple phages were no more likely to be cross-resistant to novel phages but were found to have paid a greater cost for their resistance than bacteria treated with a single phage type (Koskella et al., 2012). A similar result was observed for *Pseudomonas aeruginosa* strains treated with one versus four phages (Hall et al., 2012). This emulates the common practice of using combined antibiotic treatment to decrease the likelihood of evolved antibiotic resistance (Traugott et al., 2011; Vardakas et al., 2013). Overall, a greater understanding of the costs of resistance to phage predation and of synergistic effects among phages in controlling bacterial pathogens will allow for a more informed development and application of treatment, and ideally the prevention of widespread resistance.

## **THE IMPLICATIONS OF EVOLVED RESISTANCE TO PHAGES**

Despite the promise of many phage therapy trials, the use of phages to control bacterial pathogen begs the question: could the evolution of phage resistance mirror the evolution and spread of antibiotic resistance? Numerous studies have shown that natural phages are well-adapted to their local bacterial populations (Vos et al., 2009; Koskella et al., 2011) and that bacteria in turn adapt to resist their local phages (Kunin et al., 2008; Koskella, 2013). However, a recent review of phage resistance as a result of prolonged phage therapy (Ormälä and Jalasvuori, 2013) concludes that, as it is possible to isolate phages infective to bacteria from different geographical locations and evolutionary histories (e.g., Flores et al., 2011), long-term resistance need not be a concern as a diverse set of phages capable of infecting newly resistant strains will always be available. Local phage diversity is often high (Breitbart and Rohwer, 2005), so infective phages should be easy to isolate from just a few environmental samples. However, this parallels the problems of antibiotic discovery – the process from discovery to a useable product is arduous and expensive, so despite a ready source of infective phages few companies are investing in treatments (Brüssow, 2012). If bacterial resistance to phage infection emerges rapidly and production is slow, redundancy of treatments seems likely. As pointed out by Pirnay et al. (2011) a reactive phage therapy program that is capable of rapidly isolating, screening, and applying infective phages will be better placed to respond to phage resistance than the slow and expensive process of approval and licensing for each phage type.

Currently the maximum breadth of bacterial resistance to phage (i.e., the number of phage types a single bacterium is capable of resisting) remains largely unknown, as novel genera of phages are continually being discovered (Holmfeldt et al., 2013). For example, the ubiquitous marine bacterial clade SAR116 was thought to be so abundant as a result of escaping phage predation, but a recent finding shows that it is indeed infected by phages, and that these phages are likely to be the most abundant species on the planet (Kang et al., 2013). Our knowledge of phage ecology and evolution is still in its infancy; the exact mechanisms of infection, and in turn resistance, are often unknown and could be simultaneously diverse among strains and yet largely conserved across genera (Koskella and Meaden, 2013). There are a number of published cases of phages that are capable of infecting bacterial hosts across genera (**Table 1**), suggesting the potential for shared resistance mechanisms. Even if unlikely, evolved resistance to the few phage therapy products available to clinicians would severely impair treatment potential. This problem may be exacerbated by the more stringent control of phage products for clinical use, and thus the slow pipeline from isolation to delivery, relative to the approval of cocktails for use in agriculture. As such, rapidly responding regulation, like the measures in place for seasonal flu vaccines (Verbeken et al., 2012), could be a more effective way of countering phage therapy product redundancy.

"fmicb-04-00358" — 2013/11/28 — 20:30 — page 3 — #3


"fmicb-04-00358" — 2013/11/28 — 20:30 — page 4 — #4

Finally, the combined use of phages and antibiotics has shown great promise due to the negative pleiotropic effects of phage resistance and antibiotic resistance. Experimental evolution has demonstrated that phages applied to populations of *Pseudomonas fluorescens* that had evolved antibiotic resistance reduced population densities to a greater degree than when applied to sensitive strains (Escobar-Páramo et al., 2012). In addition, combined treatment was shown to drastically hinder the evolution of bacterial resistance over time compared to antibiotic treatment alone (Zhang and Buckling,2012). In poultry, the combination of phages and enrofloxacin resulted in lower mortality in infected birds than either treatment individually (Huff et al., 2004). Therefore, one potential step forward in controlling the spread of both antibiotic resistance and phage resistance in the environment and/or under clinical settings could be the carefully planned combination treatments of the two.

## **PHAGE-MEDIATED ATTENUATION OF BACTERIAL VIRULENCE**

A common refrain of phage therapy and biocontrol is that even if resistance does emerge, such resistance is likely to be costly, and as such would attenuate bacterial virulence in a eukaryotic host (Inal, 2003; Hagens et al., 2004). Phage resistance does seem to be correlated with a reduction in metabolic fitness both in the lab (Bohannan and Lenski, 2000; Koskella et al., 2012) and the soil (Gómez and Buckling, 2011), however the effect this may have on virulence in a eukaryotic host is largely unknown. Of the few examples, phage-resistant strains of *Yersinia pestis* have been shown to have attenuated virulence in a mouse model system, resulting in a significant increase in time to death, and in some cases complete attenuation (Filippov et al., 2011). The same phenomenon has been observed in aquaculture with a direct correlation observed between phage resistance and complete attenuation of *F. columnare* in a zebrafish system, such that one phage-sensitive phenotype resulted in 100% mortality compared to 0% in the phage-resistant phenotypes (Laanto et al., 2012).

In plant pathogens, there is also evidence to suggest that phagemediated selection might alter the infectivity and/or virulence of bacterial pathogens. Phages infecting the bacterium, *Erwinia amylovora*, were found to be most efficient at infecting strains that produced either high or low levels of exopolysaccharides (depending on the phage family examined), suggesting strong and context-dependent selection on a trait that is also known to play a key role in virulence on the plant host (Roach et al., 2013). Furthermore, Hosseinidoust et al. (2013) found that in tissue culture, phage-resistant variants of *Pseudomonas aeruginosa* actually secrete higher levels of virulence factors and caused more damage to cultured mammalian cells. With no general pattern yet emerging, phage-mediated attenuation of virulence seems hard to predict and certainly not guaranteed. It is possible that a bacterial pathogen could evolve resistance to a phage therapy product and maintain or even attain high virulence levels. Furthermore, the role of compensatory mutations to phage resistance is unknown; in the case of antibiotics such mutations can rapidly ameliorate the costs paid for resistance (Levin et al., 2000; Brandis and Hughes, 2013) and the same may be true of resistance to phages. If this turns out to be the case, phage therapy treatments that rely on the loss of costly resistance (or interactions among costly mutations to multiple phages) will be at constant risk of bacterial escape.

## **INCREASING HORIZONTAL GENE FLOW THROUGH PHAGE APPLICATION**

In addition to the direct effects of phage application on the densities and relative frequencies of bacterial pathogens, we must also be aware of the potential dangers of phage-mediated horizontal gene transfer (HGT) among pathogenic and non-pathogenic bacterial species. This is particularly relevant as HGT is an important driver of antibiotic resistance evolution (Courvalin, 1994). Given that phages facilitate horizontal gene flow through the process of transduction and that beta-lactam resistance genes have been isolated from environmental phage genomes (Colomer-Lluch et al., 2011), there is clear need to be cautious in our application of phages in the environment. Moreover it has recently been shown that antibiotic treatment itself expands the number of genes that confer drug resistance in phage metagenomes (Modi et al., 2013). Environmental perturbation with antibiotics also expanded the ecological network of phages, suggesting that drug treatment selects for greater phage-mediated transfer of resistance genes (Modi et al., 2013). Finally, antibiotic treatment of swine increased the induction of prophage – another key mechanism of HGT (Allen et al., 2011).

Phage-mediated transfer of virulence factors is also a key concern, and has been well characterized in the *Vibrio cholerae* system whereby CTX phage (among others) shapes the severity of cholera pandemics through the transfer of toxin producing virulence factors and environmental fitness benefits (Faruque and Mekalanos, 2012). Additionally, *V. cholerae* has been shown to become naturally competent on a chitin surface (similar to its environmental niche on crustacean exoskeletons), facilitating uptake of exogenous DNA from an array of sources (Meibom et al., 2005), which could be important if phages are lysing nearby pathogenic bacterial strains. Environmental perturbation through unnaturally high titers of phages could lead to high levels of transduction and horizontal gene flow with unintended outcomes, exacerbating antibiotic resistance and moving virulence factors and toxins among genomes. However, a better understanding of phage host range and host range expansion could help mitigate the spread of genes among bacterial hosts.

## **IMPACT OF PHAGE APPLICATION ON NATURAL MICROBIAL COMMUNITIES**

The importance of the microbial flora to the fitness of human hosts has become clear in recent years, with microbes playing a proposed role in obesity (Greenblum et al., 2012), oral health (Belda-Ferre et al., 2012), Crohn's disease (Manichanh et al., 2006), AIDS (Saxena et al., 2012), and even mental health (Foster and McVey Neufeld, 2013). Although less well studied, a similar role of microbiota in shaping fitness is likely to be true for agricultural plants (Peiffer and Ley, 2013) and livestock (McFall-Ngai et al., 2013). Whilst antibiotics can disrupt a large proportion of the microbial community (Dethlefsen and Relman, 2011), phage therapy may facilitate the targeted elimination of a pathogenic strain without disruption of the normal microbiota of a patient (Carlton, 1999) or a crop plant. Indeed, this argument has been a focal one to support phage therapy over antibiotic use (Loc-Carrillo and Abedon, 2011). The term "dysbiosis" has been introduced to describe a microbial flora that has become "unhealthy," typically in human disease (Tamboli, 2004). The same phenomenon is likely to hold true for natural microbial communities both within and outside of the host environments, and thus the addition of foreign phages in the form of a biopesticide could destabilize such communities, causing "dysbiosis" and potentially having subsequent effects on disease and nutrient cycling. For example, in soil where ratios of phages to bacteria are expected to be near 1:1 (Reyes et al., 2010) an influx of applied phages could well disrupt stable ecological communities important in nutrient cycling. Conversely, in the marine environment, where this ratio is closer to 1:100 (Reyes et al., 2010), an influx of applied phages is unlikely to adversely disrupt a normal microbial community. Unfortunately, by the same logic, aquaculture may suffer from lower chances of success from phage biocontrol, as the concentrations of phages needed to influence bacterial density is likely to be a limiting factor.

## **UNPREDICTABILITY OF INFECTION KINETICS**

Unlike antibiotic usage, where an effective dosage can be determined for different species and corresponding quantities used, the infection kinetics of phage therapy are less predictable (Payne and Jansen, 2003). This knowledge gap presents a challenge for effective use balanced with environmental safety, as an individual phage persists in the environment for relatively short timescales, but a lineage can reproduce indefinitely when hosts are plentiful. This continued replication makes calculations of dosage difficult. In antibiotic treatment the minimum inhibitory concentration (MIC) is crucial to informing proper antibiotic dosages; however, finding the right balance between a high enough titer of phages to be effective and not introducing excessive levels into the environment is difficult. Potentially the application of a single phage could lead to the continued replication of an infective phage, thus perpetuating not only treatment but also persistence in the environment. If there are unwanted effects of phage biopesticide there are no tools currently available to selectively remove such viruses from the environment. Also, the timing of treatment as a function of bacterial density can be crucial for a successful outcome (Payne and Jansen, 2003). Indeed, in a trial controlling*V. parahaemolyticus* infections in shrimp timing was crucial whereas changes in dosage had no effect. Reducing the multiplicity of infection made no difference whereas when treatment was delayed it was ineffective in controlling mortality (Martínez-Díaz and Hipólito-Morales, 2013), possibly as a result of non-linear infection kinetics. Such infection kinetics of phage therapy were tested *in vitro* with *Campylobacter jejuni*, a common poultry pathogen, and found to fit a non-linear model with a density-dependent proliferation threshold (Cairns et al., 2009). These studies highlight the nuances involved in effective phage therapy use and as such the use of phages cannot be analogous to that of antibiotics. Provided policy makers accept this and approach the field balancing the risks of widespread environmental use with the pressing need for a clinical alternative to antibiotics, phage biocontrol could be an integral tool in controlling bacterial diseases.

## **FUTURE OUTLOOK**

"fmicb-04-00358" — 2013/11/28 — 20:30 — page 5 — #5

Our understanding of the biology, ecology, and evolution of microbial pathogens has improved immeasurably since the advent of widespread antibiotic use. If we can learn the lessons from our mistakes with antibiotic use, phage therapy and biocontrol could form an integral tool in the fight against bacterial infections that threaten human health and food production (Pirnay et al., 2012; Allen et al., 2013). For example, the falling costs of whole genome sequencing (Kisand and Lettieri, 2013) should make tracking the evolution and spread of resistance genes in a clinical setting easier and more accurate (Didelot et al., 2012). Furthermore, advances in metagenomics may make monitoring the effects of environmental perturbations on microbial communities feasible and allow researchers to track changes over long-time periods. An attractive avenue of research for pharmaceutical companies may be the use of phage lysins – enzymes that are capable of bursting bacterial cells open "from without." Such an approach avoids the problems of infection kinetics mentioned previously and can have a broad-spectrum encompassing multiple strains of antibiotic resistant pathogens, such as MRSA and vancomycin intermediate *Staphylococcus aureus* (Gilmer et al., 2013). The downside is that lysins, unlike phages, lack the ability to counter evolve to pathogens.

Microbial biofilms present a continued risk to healthcare as they may harbor bacteria in a less metabolically active state that survive antibiotic treatment (Oliver, 2010). Phages targeting biofilms in synchrony with antibiotics may form a novel strategy, although the inherent risks of HGT still remain. It may also be possible to circumvent this cycle of treatment, selection for resistance and re-infection through the use of "social disruption" treatments that reduce bacterial virulence without selecting for resistance (Boyle et al., 2013). Phages could play a role in the reduction of biofilm and public good production – one example is an engineered phage that expresses a biofilm-degrading enzyme (Lu and Collins, 2007). Whilst this is likely to reduce the fitness of bacterial populations, selection should be weaker than that of an antibiotic. Given the difficulties faced by clinicians in treating antibiotic-resistant infections and the urgency of finding alternative therapies, the prudent use of phages should be a priority. Nevertheless we have the tools to track resistance and simple measures such as providing a diverse set of phages for treatment could help. A seasonal-vaccine-like scheme could create a treatment program that is responsive to the evolution of resistance. In short, a very different model to that of our use of antibiotics.

## **REFERENCES**


"fmicb-04-00358" — 2013/11/28 — 20:30 — page 6 — #6

Gómez, P., and Buckling, A. (2011). Bacteria-phage antagonistic coevolution in soil. *Science* 332, 106–109. doi: 10.1126/science.1198767


"fmicb-04-00358" — 2013/11/28 — 20:30 — page 7 — #7


**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.

*Received: 06 September 2013; accepted: 12 November 2013; published online: 29 November 2013.*

*Citation: Meaden S and Koskella B (2013) Exploring the risks of phage application in the environment. Front. Microbiol. 4:358. doi: 10.3389/fmicb.2013.00358*

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

*Copyright © 2013 Meaden and Koskella. 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.*

"fmicb-04-00358" — 2013/11/28 — 20:30 — page 8 — #8

## Stormwater runoff drives viral community composition changes in inland freshwaters

## *Kurt E. Williamson\*, Jamie V. Harris , Jasmin C. Green , Faraz Rahman and Randolph M. Chambers*

*Department of Biology, The College of William and Mary, Williamsburg, VA, USA*

#### *Edited by:*

*Stephen Tobias Abedon, The Ohio State University, USA*

#### *Reviewed by:*

*Lee Ann McCue, Pacific Northwest National Laboratory, USA Marla Tuffin, University of the Western Cape, South Africa*

#### *\*Correspondence:*

*Kurt E. Williamson, Department of Biology, The College of William and Mary, 3035 Integrated Science Center, Williamsburg, VA 23185, USA*

*e-mail: kewilliamson@wm.edu*

Storm events impact freshwater microbial communities by transporting terrestrial viruses and other microbes to freshwater systems, and by potentially resuspending microbes from bottom sediments. The magnitude of these impacts on freshwater ecosystems is unknown and largely unexplored. Field studies carried out at two discrete sites in coastal Virginia (USA) were used to characterize the viral load carried by runoff and to test the hypothesis that terrestrial viruses introduced through stormwater runoff change the composition of freshwater microbial communities. Field data gathered from an agricultural watershed indicated that primary runoff can contain viral densities approximating those of receiving waters. Furthermore, viruses attached to suspended colloids made up a large fraction of the total load, particularly in early stages of the storm. At a second field site (stormwater retention pond), RAPD-PCR profiling showed that the viral community of the pond changed dramatically over the course of two intense storms while relatively little change was observed over similar time scales in the absence of disturbance. Comparisons of planktonic and particle-associated viral communities revealed two completely distinct communities, suggesting that particle-associated viruses represent a potentially large and overlooked portion of aquatic viral abundance and diversity. Our findings show that stormwater runoff can quickly change the composition of freshwater microbial communities. Based on these findings, increased storms in the coastal mid-Atlantic region predicted by most climate change models will likely have important impacts on the structure and function of local freshwater microbial communities.

**Keywords: viruses, RAPD-PCR, community ecology, particle-associated, freshwater, terrestrial runoff, sediment resuspension**

## **INTRODUCTION**

Viruses are the most abundant organisms in aquatic systems, typically outnumbering bacteria 10- to 100-fold, and the majority of environmental viruses are believed to infect prokaryotes (Wommack and Colwell, 2000; Weinbauer, 2004). Viral infections can impact microbial community function by influencing abundance and species composition of host communities (Thingstad, 2000; Bonilla-Findji et al., 2009; Sandaa et al., 2009). The release of dissolved organic compounds from viral lysis of cells can also influence global nutrient cycles (Azam et al., 1983; Bratbak et al., 1992; Fuhrman, 1999; Bonilla-Findji et al., 2008). In pelagic marine waters, viruses are thought to infect primarily co-occurring bacterial hosts (Proctor and Fuhrman, 1990). Novel virus types may be introduced into freshwater systems via groundwater discharge (Yates et al., 1985; Abbaszadegan et al., 1993), sewage outfall (Krikelis et al., 1986), or land-based stormwater runoff (Rajal et al., 2007).

Storms and land-based stormwater runoff can significantly impact coastal and inland waters. Runoff conveys an astonishing array of materials from terrestrial sources to aquatic sinks. A partial list of stormwater runoff components includes solids such as sand, silt, clay, particulate organic matter, gravel, and trash; chemicals, including nitrogen and phosphorus compounds, soluble organic matter, pesticides, hydrocarbons; metals such as Zn, Pb, Hg, Fe; and microbes, including protozoa, bacteria, and, as mentioned above, viruses (U.S. EPA, 2007, 2008). Research regarding the impacts of microbes in stormwater runoff on receiving waters has historically focused on human pathogens (Geldreich, 1996; Haile et al., 1999; Ferguson et al., 2003; Arnone and Walling, 2007; Davies et al., 2008; Viau et al., 2011; Teng et al., 2012; Shapiro et al., 2013). Typically, fecal coliform plate counts (Knight et al., 2000; Jeng et al., 2005; Schoonover and Lockaby, 2006) or PCR detection of specific viral nucleic acid sequences (Rajal et al., 2007) are used to assess the microbiological quality of runoff-impacted waters, and identify potential threats to human health.

To the best of our knowledge, there has been only one investigation into the impacts of runoff on natural communities of viruses in freshwaters. Hewson et al. determined the potential contribution of terrestrially derived viruses to the viral communities of two freshwater lakes in New York state (Hewson et al., 2012). The research team constructed two metagenomic libraries for each lake, one representing surface water and one representing catchment soils. After sequence analysis of each library, specific sequences were selected to serve as indicators of viral origin, with certain sequences associated with terrestrial viruses, and certain sequences associated with aquatic viruses. By designing qPCR primers for these two groups of viruses, the team tracked the abundance of each virus group in the two lakes over a 5 week monitoring period. Aquatic virus sequences amplified consistently across all samples during the observation period, but the terrestrial virus sequences were only detected after rainfall-runoff events. These results provide strong evidence that stormwater runoff regularly introduces potentially novel virus types to aquatic habitats. However, the impacts of stormwater runoff on aquatic viral community dynamics and composition are still poorly understood.

The main goal of this paper was to determine how freshwater viral communities respond to the specific disturbance represented by storms and the influx of stormwater runoff. We first estimated the actual microbial abundance carried in stormwater runoff by sampling directly in an erosion channel during a storm. We then conducted field studies in a stormwater retention pond over the courses of two major storms (Hurricane Sandy and Tropical Storm Andrea) to quantify changes in the viral community, including the particle-associated fraction, due to stormwater runoff. Viral abundance was measured using epifluorescence microscopy and changes in viral community composition were determined using RAPD-PCR. The community and environmental data were then analyzed to identify specific environmental factors that best explained observed changes.

## **MATERIALS AND METHODS STUDY SITES AND SAMPLE COLLECTION** *Site 1: Charles City County (CCC) plantation*

Owing to the diffuse nature of non-point source runoff, sampling can be problematic. To assess the microbial load carried by surface runoff, we chose to collect and analyze water from an ephemeral channeldraininga21-haagriculturalwatershedatCCCPlantation, Charles City, VA (**Figure 1A**). Channel water was collected at discrete time intervals during Tropical Storm Lee (September 6–8, 2011) in polypropylene bottles using an ISCO 3700 automated sampler, and stored at −20◦C following retrieval. These samples were kindly provided by James M. Kaste and Gregory S. Hancock, Department of Geology, College of William and Mary.

## *Site 2: Grim dell*

The "Grim Dell" (**Figure 1B**, inset) is a 2316 m2 wet stormwater retention pond (max. depth 2.37 m) on the William and Mary campus. The pond drains a 22.77 ha watershed and discharges to Crim Dell creek and Lake Matoaka. Surface water samples were collected near the concrete weir (37◦16 16.8 N, 76◦42 59.9 W) in polypropylene bottles using an ISCO 3700 automated sampler (Teledyne Isco, Lincoln, NE) during "Superstorm" Sandy (October 26–November 1, 2012) and stored at −20◦C until analysis (approximately 1 month). Surface water samples were collected by hand in sterile polycarbonate bottles during Tropical Storm Andrea (June 3–14, 2013) and during a comparatively dry period (July 16–19, 2013). All samples were immediately frozen in liquid nitrogen and stored at −80◦C until analysis (approximately 1 month).

#### **VIRAL AND BACTERIAL ABUNDANCE**

Water samples (2 ml) were dispensed into cryovials, flashfrozen in liquid nitrogen, and stored at −80 C until use (Wen

**FIGURE 1 | Sample sites. (A)** CCC Plantation, white circle indicates sampling location, white lines indicate flow pathways; **(B)** Grim Dell stormwater retention pond (inset), white circle indicates sampling site, box indicates concrete weir, white star indicates outflow.

et al., 2004). Viral and bacterial abundances were determined as described in Hardbower et al. (2012). Briefly, dilutions of thawed, whole water samples were passed through 0.02 μm Anodisc filters (Whatman, Maidstone, England), stained with 1× SYBR Gold (Life Technologies, Grand Island, NY), and bacterial cells and virus-like particles were enumerated using epifluorescence microscopy. Particle-associated microbial abundances were determined by centrifuging 200 ml of whole water at 8000 ×g for 20 min at 4◦C to collect total suspended solids (TSS). The supernatant was decanted and the pellet was resuspended in 50 ml sterile potassium citrate buffer (Williamson et al., 2013). The slurry was sonicated in an ice bath 3 × 1 min at 60 W using a Branson 250 sonication probe outfitted with a ¼" micro-tip, with 30 s of manual shaking in between each round of sonication (Danovaro et al., 2001; Williamson et al., 2007). Aliquots were then diluted 100- to 250-fold and microbial abundances were determined as described above. Subsamples (50 ml) of the sonicated slurry were centrifuged at 5000 ×g for 20 min to pellet suspended solids and the supernatant was filtered through 0.22μm Sterivex syringe filters (Millipore, Billerica, MA) to remove bacteria. Viral concentrates (VCs) were prepared as described below.

## **WATER CHEMISTRY AND PHYSICAL PROPERTIES**

A YSI-63 hand-held multimeter was used to measure temperature and conductivity and a YSI-55 hand-held probe was used to measure dissolved oxygen in the field (YSI Inc., Yellow Springs, OH). Nutrient concentrations were determined through colorimetric assays using GF/F (glass fiber)—filtered water (Parsons et al., 1984). Rainfall data were obtained from the National Atmospheric and Oceanic Administration (http://www*.*noaa*.* gov/wx*.*html).

## **ANALYSIS OF VIRAL COMMUNITY COMPOSITION BY RAPD-PCR**

VCs were prepared using ultracentrifugation as in Hardbower et al. (2012). Aliquots (50 ml) of water samples were filtered through 0.22μm Sterivex syringe filters (Millipore, Billerica, MA) into polyallomer ultracentrifuge tubes (Beckman-Coulter, Pasadena, CA). Tubes were spun at 22,000 rpm for 2 h at 4◦C in a Beckman SW41 Ti rotor to pellet virus particles. Supernatants were carefully decanted to avoid disturbing viral pellets, and pellets were resuspended in 20–50μl TMG buffer (10 mM Tris-Cl, 10 mM MgSO4, 1% glycerol). Sodium azide (0.1% final conc.) was added to inhibit growth of any potential bacterial contamination. VCs were confirmed to be free of microbial contaminants by epifluorescence microscopy as previously described (Helton and Wommack, 2009; Winter and Weinbauer, 2010; Hardbower et al., 2012) VCs were stored at 4◦C until use, generally within 1 month of sample collection. RAPD-PCR reactions were set up using primer CRA-22 (5 -CCGCAGCCAA-3 ) and thermocycler conditions were programmed as described in Winget and Wommack (2008). Products from RAPD-PCR were separated by gel electrophoresis on 13 × 16 cm 1.8% MetaPhor agarose gels (Lonza, Aplharetta, GA) in 0.5<sup>×</sup> TBE buffer, run at 4 V cm−1. Gels were stained with 1× SYBR Safe (Invitrogen, Eugene, OR) in 0.5× TBE for 1 h prior to visualization of bands using a Kodak Gel Logic 100 imaging system. Banding patterns were analyzed using ImageQuant TL software (GE Life Sciences, Piscataway, NJ) and converted to a binary matrix representing the presence-absence of viral operational taxonomic units (OTU) (Hardbower et al., 2012).

## **STATISTICAL ANALYSES**

Environmental data were checked for normality using D'Agostino and Pearson omnibus normality test; viral and bacterial abundance values failed the test and were log10-transformed to obtain normal distributions. Binary matrix data representing the presence-absence of viral taxa (RAPD banding patterns) were converted into dissimilarity matrices (Dice method). Detrended correspondence analysis (DCA) of viral community matrices and environmental data matrices was performed using R (v.2.12.1). For all data sets, the length of the first DCA axis was *<*2 and, thus, a linear relationship between species and environmental variables was assumed (Jongman et al., 1995). To determine the most suitable set of parameters that explained viral abundance, bacterial abundance, and viral richness, stepwise multiple regression analyses were performed. Multicollinearity of multiple regression models were evaluated by calculating the variance inflation factor (VIF) as given by VIF*<sup>j</sup>* <sup>=</sup> <sup>1</sup>*/*<sup>1</sup> <sup>−</sup> *<sup>R</sup>*<sup>2</sup> *<sup>j</sup>* , using package fmsb in *R*. VIF values (*<*5 in all cases) showed that the stepwise multiple regression analysis was not affected by collinearity of the parameters. The effects of the parameters on the multiple regression models having a *p*-value *<*0.05 were assumed to be significant. Relationships between viral community data and explanatory environmental variables were analyzed by redundancy analysis (RDA). Environmental variables best describing changes in viral community composition were identified by forward selection. Explanatory variables were added until further addition of variables failed to contribute to a significant improvement to the model's explanatory power.

Dice dissimilarity matrices representing viral community data were analyzed using non-metric multidimensional scaling (NMDS) in PAST (Hammer et al., 2001). NMDS was used to detect patterns that could explain the observed similarities or dissimilarities (distances) among samples. In NMDS plots, the closer two samples are plotted together, the more similar their viral community compositions, and the more distant two samples are from each other, the more dissimilar their viral community composition. The lower the stress value, the better the goodness of fit for the overall model.

Multivariate regression trees (MRTs) were used to determine the degree to which time and different environmental factors (explanatory variables) were predictive of the viral community composition (response variable). MRT is a robust predictive method even when high-order interactions exist among explanatory variables (De'ath, 2002), thus it is appropriate to our analysis. MRTs were constructed using the mvpart package in R with 100 cross-validations to select the best tree. The tree is constructed by repeated binary splitting of the data, where each split is defined by a simple rule, usually based on one to a few explanatory variables, and forms two nodes. Splits are chosen to maximize the homogeneity of the resulting two nodes. The terminal nodes (leaves) represent the groups of data formed by the tree. The depth of the tree following each split is proportional to the variance explained by the split, and the cross-validated relative error is an indicator of the tree's value in predicting changes in the response variable, where 0 = perfect prediction and 1 = no predictive value.

## **RESULTS**

## **CCC PLANTATION: VIRAL LOAD CARRIED BY STORMWATER RUNOFF**

Tropical Storm Lee delivered 214.12 mm of rainfall over 67 h and a total of 12,700 m<sup>3</sup> of stormwater passed through the erosion channel at the CCC Plantation site during the observation period (Caverly et al., 2013). Three water samples, collected directly within the ephemeral channel, were analyzed to determine the viral load carried by stormwater runoff. The concentration of planktonic viruses varied between 1 <sup>×</sup> <sup>10</sup><sup>7</sup> and 5 <sup>×</sup> 10<sup>7</sup> ml−<sup>1</sup> (**Figure 2A**). Stormwater runoff also contained many suspended particles with microbes attached (**Figure 2B**). The particle-associated viral abundance within stormwater runoff was 2.02–7*.*<sup>8</sup> <sup>×</sup> <sup>10</sup><sup>7</sup> ml−1, comprising 26.4–75.6% of the total viral abundance in the runoff (**Figure 2A**). Given the concentrations of viruses in stormwater runoff, including both the planktonic and

particle-associated fractions, and the volume of water flux during the storm, approximately 10<sup>15</sup> virus particles passed through the erosion channel during this single storm event.

#### **RESPONSE OF GRIM DELL RETENTION POND TO HURRICANE SANDY**

Surface water samples were collected at regular 6-h intervals for the duration of Hurricane "Superstorm" Sandy (October 26–November 1, 2012). During the 138-h observation period, rainfall began at 34.5 h and averaged about 2 mm h-1 over the next 36 h (**Table 1**). Peak rainfall occurred at 72 h, with almost 28 mm of rain in the preceding 6-h period. Rainfall then tapered off, concluding by 108 h into the observation period. Water pH varied from a low of 6.17 after peak rainfall to a high of 6.84 just prior to the beginning of rainfall (**Table 1**). Conductivity followed a similar trend with the lowest value (73μS cm−1) occurring just after peak rainfall and the highest values (*>*700μS cm−1) observed prior to rainfall. Ammonium concentrations varied from below the limit of detection to a maximum of 10.1μM at the 6 h time point, but showed no clear trend over time or with regard to rainfall. Phosphate concentrations were generally lower at the beginning and end of the observation period and higher during rainfall, peaking at 4.5μM. Nitrate and nitrite concentrations were generally lower during the beginning of the observation period, ranging from 0.5 to 2.8 μM for the first 66 h, then increasing to a maximum concentration of 12.9μM after peak rainfall. The high levels of TSS in the first three time points were due to improper positioning of the ISCO sampler intake too close to bottom sediments. Taking this into consideration, suspended solids varied from a low of 6 mg l <sup>−</sup><sup>1</sup> just prior to the beginning of rainfall, with the maximum concentration of 102 mg l−<sup>1</sup> coinciding with peak rainfall, and steadily decreasing concentrations to the end of the observation period.

In spite of the large amount of precipitation, microbial abundance in the pond only varied by a factor of about 3 over the observation period, with viral abundance ranging from 1.4 to <sup>3</sup>*.*<sup>6</sup> <sup>×</sup> <sup>10</sup><sup>6</sup> ml−<sup>1</sup> and bacterial abundance ranging from 1.5 to <sup>3</sup>*.*<sup>2</sup> <sup>×</sup> <sup>10</sup><sup>5</sup> ml−<sup>1</sup> (**Figure 3A**). Viral abundance was consistently higher than bacterial abundance by a factor of about 10. Stepwise multiple regression analysis indicated that variation in viral abundance was best explained by the combined effects of conductivity, pH, and bacterial abundance (*R*<sup>2</sup> <sup>=</sup> <sup>0</sup>*.*85, *<sup>p</sup> <sup>&</sup>lt;* <sup>0</sup>*.*001), but no suitable model was found to explain changes in bacterial abundance (**Table 2**).

Based on the observation from the CCC Plantation samples that runoff contains high concentrations of particleassociated viruses, a subset of the Grim Dell time series was analyzed to determine the abundances of particle-associated microbes (**Figure 3A**). Particle-associated viral abundances varied from 1*.*<sup>32</sup> <sup>×</sup> <sup>10</sup><sup>6</sup> to 1*.*<sup>35</sup> <sup>×</sup> <sup>10</sup><sup>7</sup> ml−<sup>1</sup> and accounted for 50.9–87.4% of the total viral abundance within a given sample, while particle-associated bacterial abundances varied from 9*.*<sup>6</sup> <sup>×</sup> 104 to 4*.*<sup>6</sup> <sup>×</sup> 105 ml−1, comprising 24.5–64.8% of the total bacterial abundance. Not surprisingly, peaks in the concentration of TSS coincided with peaks in the abundance of particle-associated viruses and bacteria (**Table 1**, **Figure 3A**). As observed in the CCC Plantation samples, the abundances of particle-associated microbes during the storm were generally, but not always, higher than the abundances of planktonic microbes. Spearman correlations were performed on a subset of the data for which particleassociated viral and bacterial abundances were available. Strong and significant relationships were identified between particleassociated viruses and particle-associated bacteria (*r* = 0*.*857, *p* = 0*.*024), between particle-associated viruses and particulate phosphorus (*r* = 0*.*929, *p* = 0*.*007) and between particleassociated bacteria and particulate phosphorus (*r* = 0*.*786, *p* = 0*.*045).

VCs were prepared from a subset of all samples and viral community composition was compared based on RAPD-PCR banding patterns. Band richness varied from 7 to 11 distinct band types (**Table 1**). Stepwise multiple regression analysis indicated that viral richness was best explained by combined effects of conductivity, pH, nitrate + nitrite, and particulate phosphorus, although this model was not particularly robust (*R*<sup>2</sup> <sup>=</sup> 0*.*60, *p* = 0*.*027; **Table 2**). NMDS analysStepwise multiple regression is of changes in RAPD banding patterns demonstrated that disturbance from Hurricane Sandy resulted in changes in viral community composition (**Figure 3B**). The first three time points show relatively little change in viral community composition over the 18 h prior to rainfall. At 34.5 h, rainfall begins and at 36 h the community composition begins to shift more dramatically. Large shifts are associated with the initial rainfall (36–42 h), and following more intense periods of rain (72–84 h, **Figure 3B**). Following the conclusion of the storm, the viral



*TSS, total suspended solids; asterisk indicates questionable reading due to sample probe initially lodged in sediment; Partic P, particulate phosphorus; BDL, below detection limit; ND, not determined.*

community continued to change over time up to the 138 h time point.

MRT analysis revealed a clear relationship between time, rainfall, and viral community composition. RDA indicated that time was the single best explanatory factor in describing the viral community data (**Table 3**). Similarly, MRT indicated that time was the single best factor predictive of viral community composition (**Figure 4**). By cross-comparing the MRT with **Figure 3**, a clear relationship exists between the time intervals represented in the regression tree and storm stage (rainfall). The largest discontinuity in the viral community data occurred between 66 and 72 h, coinciding with peak rainfall, while the leaves of the tree represent pre-storm, rising intensity, falling intensity, and post-storm conditions (**Figures 3**, **4**).

A subset of samples (those for which particle-associated abundances were determined, **Figure 3A**) was analyzed to compare the community composition of the particle-associated viruses. Comparisons of RAPD banding patterns across particleassociated VCs indicated that the composition of particleassociated viral communities in the Grim Dell changed over time (**Figure 5**). Furthermore, comparison of RAPD banding patterns between the particle-associated viruses and the planktonic viruses revealed that the two viral communities were completely distinct (**Figure 5**).

## **RESPONSE OF GRIM DELL RETENTION POND TO TROPICAL STORM ANDREA**

Surface water samples were collected twice per day from June 3–14, 2013. Although this sample interval captured changes in the pond due to Tropical Storm (TS) Andrea (June 6–11, 2013), several other rainfall events occurred during the observation period (**Table 4**). Rainfall varied from 0 mm over 48 h or more to *>*60 mm within 6-h during the peak of the tropical storm. Because of the unusually high frequency of rainfall events during June, samples were also collected July 16–19, 2013 during a dry period to compare changes in the pond in the absence of rainfallrunoff (**Table 4**). Over these two observation periods, pond pH varied from 5.08 to 7.19, with the highest reading at the beginning of the observation period and the lowest value following a week of dry weather. Surface water temperature was lower in the mornings and could be ≥ 5◦C higher by the afternoon timepoint (approximately 6 h later) on the same day, likely due to the shallow depth of the pond. This daily variation was larger than the variation in the average temperature over the two observation periods. Conductivity was generally lower during rainfall events; however, decreases due to precipitation were not as dramatic as those observed during Hurricane Sandy in October, 2012 (**Tables 1**, **4**). Dissolved oxygen increased during rainfall, most likely due to increased mixing during storms. TSS were generally

error bars for aquatic viral and bacterial abundances represent SD (*n* = 2). **(B)** Non-metric multidimensional scaling (NMDS) plot of viral community composition as determined by RAPD-PCR banding patterns based on Dice dissimilarity. Data points are coded by number of hours elapsed in the observation period and correspond to samples listed in **Table 1**.

higher during periods of rainfall and decreased with time following precipitation events, with notable exceptions (**Table 4**). For example, on June 11, 2013, 40 mm of rainfall had occurred since the time point 24 h prior, but TSS actually decreased over this period.

As observed during Hurricane Sandy, viral and bacterial abundances did not change significantly due to precipitation and influx of stormwater runoff during TS Andrea (**Figure 6A**). Similar microbial abundances were observed under both wet (**Figure 6A**) and dry (**Figure 6B**) conditions, suggesting that rainfall-runoff does not significantly influence aquatic microbial abundance. Mann-Whitney *U*-tests were used to compare measurements between wet (June 3–14, 2013) and comparatively dry weather (July 16–19, 2013, **Table 4**). Viral and bacterial abundances and dissolved oxygen concentration were not significantly different in the pond during wet vs. dry weather. However, pH was significantly lower (*p* = 0*.*003) and surface water temperature and conductivity were significantly higher (*p* = 0*.*018, 0.010, respectively) **Table 2 | Stepwise multiple regression analysis of biological response variables in relation to environmental variables for Grim Dell retention pond during Hurricane Sandy (October, 2012).**


*Parameters defining the best performing models are shown.*

*Cond, conductivity; BA, bacterial abundance; Partic P, particulate phosphorus; Rain, rainfall; SE, standard error; ns, not significant.*

during the dry weather observation period. Stepwise multiple regression analysis indicated that variation in viral abundance was best explained by the combined effects of conductivity and bacterial abundance, but the explanatory power of this model was much weaker than that for viral abundance during Hurricane Sandy (*R*<sup>2</sup> <sup>=</sup> <sup>0</sup>*.*57, *<sup>p</sup> <sup>&</sup>lt;* <sup>0</sup>*.*002; **Table 5**). Variation in bacterial abundance was best explained by combined effects of viral abundance, nitrate + nitrite, phosphate, and dissolved oxygen, with a more robust model than that for viral abundance (*R*<sup>2</sup> <sup>=</sup> <sup>0</sup>*.*70, *p <* 0*.*001; **Table 5**).

Viral community composition was compared based on RAPD-PCR banding patterns as described above. Band richness varied from 6 to 15 distinct band types (**Table 4**), but stepwise multiple regression analysis was unable to extract a useful model to explain variability in band richness in terms of environmental factors (**Table 5**). NMDS analysis of changes in RAPD banding patterns suggested that disturbance from storms induced changes in viral community composition (**Figure 6C**). Large shifts were associated with periods of rainfall (e.g., 6-3\_0830 to 6-3\_1330), and followed intense periods of rain (e.g., 6-8\_1300 to 6-9\_1500 to 6-10\_0930; **Figure 6C**). Community composition dynamics for samples collected during dry conditions varied considerably from those collected during wet conditions. In the absence of precipitation, samples collected 4–24 h apart exhibited relatively little change in the viral community composition compared to any given pair of samples collected during wet conditions. A notable exception is the pair of samples 6-6\_0940 and 6-6\_1500 (**Figure 6C**), collected approximately 4 h apart during wet weather conditions. In this case, however, no rainfall was recorded for at least 48 h prior to sample collection and the small shift may represent a temporary stabilization in community composition.

RDA indicated that time was the single best explanatory factor in describing the viral community data (**Table 3**). As with the previous data set from Hurricane Sandy, MRT analysis indicated


**Table 3 | Results of redundancy analysis of viral community data in relation to environmental variables for Grim Dell retention pond during Hurricane Sandy (October, 2012) and Tropical Storm Andrea (June, 2013).**

and numbers in parentheses indicate the number of cases (samples) in each leaf. CV Error, cross-validated relative error, where 0 = perfect prediction and 1 = no predictive value.

that time was the single best factor predictive of viral community composition during TS Andrea (**Figure 7**). Cross-comparison with **Figure 6A** reveals a clear relationship between the time intervals represented in the MRT and storm stage (rainfall). The largest discontinuity in the viral community data occurred at 6-8\_1300, coinciding with peak rainfall during the observation period. With the exception of the discontinuity at 6-5\_1800, all other leaves in the regression tree coincide with either changes in rainfall intensity or lack of rainfall. The discontinuity observed at 6-5\_1800 may indicate further rearrangement of the viral community structure following the disturbance represented the pervious rainfall event.

## **DISCUSSION**

This work represents the first known study to investigate the total viral load of stormwater runoff and measure the impacts of terrestrial runoff on freshwater aquatic viral communities. While this is an important first step in understanding how aquatic viral communities respond to infiltration of runoff, several important caveats must be taken into consideration. Microbial loads in runoff will differ based on catchment land use and land cover (Sharma et al., 2013), rainfall duration and intensity (Ran et al., 2013), antecedent catchment conditions and nutrient concentrations (McCarthy et al., 2013), as well as watershed area and the nature of the receiving water. In the present study, all field sampling was performed during major storm events and quantitative analysis of storm impacts on aquatic viral communities is limited to a single small catchment (stormwater retention pond). Thus, while it is highly likely that most rainfall-runoff events will transfer viruses and other microbes to aquatic habitats, the specific rates and impacts of such transfers will not be the same for all storm events or catchments.

Several previous studies in both marine and freshwater habitats have shown RAPD-PCR to be a robust method for documenting temporal changes in viral community composition (Winget and Wommack, 2008; Helton and Wommack, 2009; Winter and Weinbauer, 2010; Hardbower et al., 2012), hence its use in the present work. As in previous studies (Winget and Wommack, 2008; Helton and Wommack, 2009; Winter and Weinbauer, 2010; Hardbower et al., 2012), care was taken to ensure VCs were free of microbial contamination prior to RAPD-PCR amplification so that banding patterns originated from viral templates, and replicate PCRs yielded reproducible banding patterns of ≥90% similarity (**Figure S1**). While RAPD-PCR cannot provide information on the total viral community and only captures a subset of viral richness and dynamics (Winter and Weinbauer, 2010), RAPD-PCR is more sensitive to community changes and enables viral community profiling with smaller water samples than the alternative profiling approach, pulsed-field gel electrophoresis (Hardbower et al., 2012).

## **PILOT STUDY: CCC PLANTATION**

One drawback to our investigation of the microbial load carried by stormwater runoff at CCC Plantation is that the data set only represents a limited sampling of runoff during a single storm event. Furthermore, we did not attempt to pair microbial abundance estimates with environmental data. In spite of these limitations, this pilot study did provide the first known estimates of total viral abundances carried by stormwater runoff. Our data indicated that stormwater runoff can carry substantial numbers of viruses, with concentrations equivalent to many other aquatic habitats (Wommack and Colwell, 2000; Weinbauer, 2004; Clasen et al., 2008). Therefore, runoff-based transport of terrestrial viruses to aquatic ecosystems may influence temporal changes in the viral community composition of inland waters.

## **IMPORTANCE OF PARTICLE-ASSOCIATED VIRUSES**

An important aspect of these pilot investigations of microbial abundances in stormwater runoff is estimation of the particleassociated fraction. Virus transport is known to be enhanced by binding to soil particles or aggregates, which themselves can be transported by mass flow (Chattopadhyay and Puls, 2000; Jin, 2000; Jin and Flury, 2002). However, relatively few studies have

examined particle-associated microbes in natural aquatic systems, particularly freshwaters. In organic marine aggregates, viral abundances have been reported anywhere from 3 <sup>×</sup> <sup>10</sup><sup>6</sup> to 8*.*<sup>7</sup> <sup>×</sup> 1010 ml−1, comprising between 0 and 40% of the total viral abundance in the water column (Weinbauer et al., 2009). In a freshwater riverine system, the viral abundances associated with suspended solid material was reported to be 2 <sup>×</sup> <sup>10</sup>5–5*.*<sup>4</sup> <sup>×</sup> <sup>10</sup><sup>9</sup> ml−<sup>1</sup> (Luef et al., 2007, 2009), comprising 0.4–35% of the total viral abundance. In the present study, the abundances of particle-associated microbes were well within the range of previously reported values. However, in samples collected from stormwater runoff and from the Grim Dell during storm events, particle-associated viruses made up a higher percentage of the total viral abundance (up to 87.4%) than has been previously reported. If particle-enhanced transport of microbes is favored over the transport of suspended microbes during rainfall-runoff events, this could explain the higher percentages of particle-associated viruses observed in our samples.

To the best of our knowledge, this study is the first to include comparison of viral community composition between the planktonic and particle-associated fractions in any aquatic ecosystem. It is important to acknowledge the limitations in our analysis here: since only single samples (*n* = 1) of particle-associated viruses were prepared from each sample, our analysis cannot account for potential within-sample variability. Thus, actual differences between samples may be smaller than they appear. Bearing in mind this limitation, cluster analysis of viral RAPD-PCR banding patterns suggested that within a given water sample, the taxonomic composition of the particle-associated viral assemblage was almost completely different from that of the planktonic viral assemblage (**Figure 5**). The taxa represented by particleassociated viruses were not a subset or expanded set of planktonic viruses, but a virtually non-overlapping set. This observation would argue against equilibrium partitioning of viruses between the aquatic and sorbed phases. Given that the samples in the present study were collected during a major disturbance (storm), it is possible that more even distributions between planktonic and particle-adsorbed viral taxa can exist under calmer conditions or in other aquatic systems.

Particle-associated viruses may be important because particle surfaces are known to be hotspots for phage-host interactions (Mari et al., 2007; Samo et al., 2008; Weinbauer et al., 2009) that can play out in any number of ways: particle-bound viruses may attack planktonic hosts, particle-bound host cells may be attacked by planktonic viruses, or particle-bound microbes infected by particle-bound viruses may release new virus types into the water column. Furthermore, depending upon the size distribution and sinking rate of introduced particles, some particle-associated viruses may serve as a reservoir for future viral infections. Thus, the ultimate impacts of particle-associated viruses entering aquatic systems through stormwater runoff have yet to be determined.

#### **RESPONSE OF PLANKTONIC VIRUSES TO RAINFALL-RUNOFF EVENTS**

Total viral and bacterial abundances in the Grim Dell did not vary more than 3-fold during storm events (**Figures 3A**, **6A,B**).


**Table 4 | Environmental and community data for Grim Dell retention pond during Tropical Storm Andrea (June, 2013) and dry weather (July, 2013).**

*All samples were collected in 2013. Samples are coded by m-d\_hhmm; Temp, temperature; dO*2*, dissolved oxygen; TSS, total suspended solids; Partic P, particulate phosphorus; V Rich, viral richness; BDL, below detection limit; ND, not determined.*

Although the CCC Plantation samples represent data from a different watershed, the viral abundances observed in stormwater runoff there were similar to the average viral abundance of the Grim Dell pond. If microbial concentrations in stormwater runoff are roughly equivalent to the receiving waters, then influx of microbes carried by stormwater would not significantly change aquatic microbial abundances per unit volume. While conductivity and bacterial abundance best explained variation in viral abundance during both observation periods, no shared factors were identified that best explained variation in bacterial abundance across the two observation periods (**Tables 2**, **5**). Differences in the extracted explanatory variables across storm events may be due to differences in the storms, themselves. For example, the data collection period for Hurricane Sandy bracketed a single, well-defined precipitation event, while the data collection period for TS Andrea encompassed several discrete precipitation events. These differences in frequency and intensity of rainfall would most likely result in biological and chemical changes in the Grim Dell of differing magnitudes and time-scales. While we suspect rainfall is the ultimate cause of variation in the environmental factors (e.g., conductivity, pH, nutrient concentrations), rainfall was not extracted as a direct explanatory variable in the model, most likely due to the large number of zero values in both time series.

Some previous studies of freshwater systems have identified strong linkages between viral and heterotrophic bacterial dynamics (Personnic et al., 2009; Cheng et al., 2010; Hardbower et al., 2012), while others have found stronger correlations between viral dynamics and photosynthetic autotrophs (Madan et al., 2005; Clasen et al., 2008; Tijdens et al., 2008). Unfortunately, measurements of chl-*a* or other indices of photosynthetic biomass were not incorporated as part of our experimental design. However, we were able to test for correlations between viral and bacterial abundances across two storm events and a short dry weather period. While no clear relationship was observed between viral and bacterial abundances during Hurricane Sandy, during

Tropical Storm Andrea and the subsequent period of dry weather, viral and bacterial abundances were significantly correlated to each other (Pearson *r* = 0*.*597, *p* = 0*.*001). These results suggest that rainfall-runoff can alter the linkages between aquatic viral and bacterial dynamics, but the impacts of individual storm events on microbial dynamics are clearly variable. Such variability is reasonable, given the differences in intensity and duration of the storms monitored.

Previous studies of temporal change in freshwater viral communities have typically described annual cycles of compositional change by comparing monthly samples (Filippini and Middelboe, 2007; Lymer et al., 2008; Tijdens et al., 2008). However, one study examined diurnal changes in freshwater viral community composition over 48 h, as well as at monthly intervals over 1 year (Lymer et al., 2008). In that study, the magnitude of diurnal changes could be as large as those of monthly changes within the same lake, suggesting that monthly sampling can underestimate the temporal variability in aquatic viral community dynamics. Our approaches in the present study represent more frequent and sustained sampling efforts, with samples gathered every 4–18 h for up to 10 consecutive days. The viral community composition in the Grim Dell was relatively stable over time periods of up to 48 h in the absence of rainfall-runoff (prior to Hurricane Sandy, **Figures 3B**, **4**; dry weather conditions, **Figures 6C**, **7**).

In at least some marine systems, microbial community dynamics exhibit reoccurring seasonal patterns in composition (Fuhrman et al., 2006; Gilbert et al., 2012; Parsons et al., 2012). Based on metagenomic analysis of marine prokaryotic communities, the reoccurrence of specific community structures in these studies has been attributed to changes in rank abundances of species that persist year-round. In other words, the taxonomic richness does not change over time, but changes in the rank abundance of individuals (distribution across taxa) leads to apparent changes in community composition (at least, as detected by coarser-grained approaches such as molecular profiling techniques), as well as changes in community function. The microbial communities of inland freshwater lakes may not share this feature, since some lakes exhibit regular annual recurrences of specific microbial taxa (Yannarell et al., 2003; Lymer et al., 2008) while others do not (Lindström, 2000; Boucher et al., 2006; Hardbower et al., 2012). The reasons behind this lake-to-lake variability are currently unclear, and indeed, the factors leading to temporal community change within any given lake remain, for the most part, poorly understood.

Our hypothesis is that changes in viral community composition observed in the Grim Dell were caused by the influx of novel virus taxa carried in stormwater runoff. This hypothesis is strongly supported by MRT analysis of viral community data. For both of the field observation periods in this study, the largest discontinuities in viral community structure coincided with peak precipitation, and the smaller discontinuities largely coincided with changes in precipitation intensity, i.e., increase, decrease, or cessation of precipitation (**Figures 4**, **7**). This hypothesis is also supported by reports from other research

**Table 5 | Stepwise multiple regression analysis of biological response variables in relation to environmental variables for Grim Dell retention pond during TS Andrea (wet + dry weather).**


*Parameters defining the best performing models are shown.*

*BA, bacterial abundance; Cond, conductivity; VA, viral abundance; dO*2*, dissolved oxygen; SE, standard error; ns, not significant.*

teams focusing on freshwater bacterial communities, rather than viruses (Lindström, 2000; Lindstrom and Bergstrom, 2005). Most recently, analysis of the small subunit ribosomal gene sequences found in watershed soils, headwater streams, and a downstream lake (Toolik Lake, AK, USA) suggested that a substantial portion of bacterial taxa found in surface freshwaters originated from terrestrial environments (Crump et al., 2012). Thus, surface runoff from watershed soils may seed downstream water bodies with new microbial taxa, including viruses.

While this is a reasonable hypothesis, storm events entail multiple simultaneous processes that may affect aquatic viral community composition, and it is difficult to isolate the impacts of runoff-borne terrestrial microbial inputs against this background. For example, runoff may have introduced chemical inducing agents into the Grim Dell, inducing bacterial lysogens to release potentially novel viral genotypes into the water column. In the present study, dramatic changes in viral community composition were observed within time intervals as short as 4 h. Well-characterized lyogenic bacteria such as *Escherichia coli* (λ) can be induced to release prophage in time intervals as small as 2 h, at least under optimal culture conditions (Little et al., 1999). However, it is unclear whether prophage induction would proceed as quickly under environmental conditions, or indeed, whether prophage release can explain all the observed variability in our data. The observance of novel viral genotypes in the water column could also have arisen through release from infected (non-lysogenic) terrestrial bacteria that were carried along in stormwater runoff. Such release may still be interpreted as the transfer of terrestrial viruses to aquatic ecosystems; rather than being transferred as extracellular particles, viruses are simply transported in a vehicle (the infected cell) and released soon after arrival.

Finally, changes in aquatic viral community composition may also be driven by resuspension of viral particles from sediments. A recent study of two lakes in upper New York state found that

virus genotypes typically associated with the water column as well as genotypes typically associated with watershed soils could both be found in lake sediments (Hewson et al., 2012). Since sediments could serve as a source of novel virus genotypes and the shallow depth of the Grim Dell (generally *<*2 m) would enable sediment resuspension during storms, we cannot rule out the possibility that changes in planktonic virus community composition were due to resuspension of viruses from sediments. In this scenario, changes in aquatic viral community composition would not be due to influx of terrestrial taxa as we hypothesized, but changes are still driven by storm impacts. While we have yet to determine the relative contributions of specific mechanisms (e.g., terrestrial runoff vs. sediment resuspension), disturbance due to storm events appears to strongly influence change in viral community composition in inland freshwaters.

Given the relationship between storm activity and aquatic viral community change, it is worthwhile considering these findings within the broader context of climate change. Most models predict increases in storm frequency and intensity in the mid-Atlantic region of the United States, where the present studies were performed (Ning et al., 2012). If the predicted changes in storm behavior are realized, this will translate to increased soil erosion, runoff, and transfer of microbes from soils to waterways, and may fundamentally change the microbial composition of inland freshwaters with unknown impacts on ecosystem function. In this regard, additional studies are needed to better characterize the impacts of storms on freshwater microbial community composition and subsequent impacts on function. Evaluating the variable impacts of stormwater runoff across different watershed types, storm durations, and storm intensities, as well as determining the ultimate impacts of these viral transport events on aquatic bacterial communities, represent important challenges for future research.

## **ACKNOWLEDGMENTS**

This work was supported by a Jeffress Memoral Trust grant (J-988) to Kurt E. Williamson and a Chappell Fellowship from the Roy R. Charles Center to Jamie V. Harris. We wish to thank Tim Russell for assistance with nutrient analyses, and James M. Kaste and Gregory S. Hancock for kindly providing runoff samples from CCC Plantation. We also wish to thank M. Drew LaMar and Matthias Leu for their helpful suggestions with statistical analyses, and Matthew Saxton and our reviewers for helpful comments in improving the manuscript.

## **SUPPLEMENTARY MATERIAL**

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

**Figure S1 | Comparison of fingerprints from randomly amplified**

**polymorphic DNA (RAPD) polymerase chain reaction (PCR).** Left half shows gel images generated by RAPD-PCR of **(A)** replicate samples and **(B)** DNase-treated samples. Viral concentrates were treated with 1 U μl −1 RQ1 DNase (Promega) according to manufacturer's instructions; in control reactions, DNase was replaced with sterile water. Sample names refer to time points in **Table 1**. R, replicate; POS, positive control; NEG, template-free control; L, molecular weight markers, indicated in base

pairs of DNA; DNase, samples treated with DNase; Ctrl, DNase-free controls. Right half shows cluster dendrograms (Dice method) of **(C)** replicate samples and **(D)** DNase-treated samples.

## **REFERENCES**


Yates, M. V., Gerba, C. P., and Kelley, L. M. (1985). Virus persistence in groundwater. *Appl. Environ. Microbiol.* 49, 778–781.

**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.

*Received: 31 October 2013; accepted: 27 February 2014; published online: 14 March 2014.*

*Citation: Williamson KE, Harris JV, Green JC, Rahman F and Chambers RM (2014) Stormwater runoff drives viral community composition changes in inland freshwaters. Front. Microbiol. 5:105. doi: 10.3389/fmicb.2014.00105*

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

*Copyright © 2014 Williamson, Harris, Green, Rahman and Chambers. 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.*

## Metagenomic and whole-genome analysis reveals new lineages of gokushoviruses and biogeographic separation in the sea

## *Jessica M. Labonté1† and Curtis A. Suttle1,2,3,4\**

*<sup>1</sup> Department of Microbiology and Immunology, University of British Columbia, Vancouver, BC, Canada*

*<sup>2</sup> Department of Earth, Ocean and Atmospheric Sciences, University of British Columbia, Vancouver, BC, Canada*

*<sup>3</sup> Department of Botany, University of British Columbia, Vancouver, BC, Canada*

*<sup>4</sup> Canadian Institute for Advanced Research, University of British Columbia, Vancouver, BC, Canada*

#### *Edited by:*

*Heather K. Allen, National Animal Disease Center, USA*

#### *Reviewed by:*

*Andrew D. Millard, University of Warwick, UK Henry M. Krisch, Centre National de la Recherche Scientifique, France*

#### *\*Correspondence:*

*Curtis A. Suttle, Department of Earth, Ocean and Atmospheric Sciences, University of British Columbia, #2178-2207 Main Mall, Vancouver, BC V6T 1Z4, Canada e-mail: suttle@science.ubc.ca*

#### *†Present address:*

*Jessica M. Labonté, Bigelow Laboratory for Ocean Sciences, East Boothbay, ME, USA*

Much remains to be learned about single-stranded (ss) DNA viruses in natural systems, and the evolutionary relationships among them. One of the eight recognized families of ssDNA viruses is the *Microviridae*, a group of viruses infecting bacteria. In this study we used metagenomic analysis, genome assembly, and amplicon sequencing of purified ssDNA to show that bacteriophages belonging to the subfamily *Gokushovirinae* within the *Microviridae* are genetically diverse and widespread members of marine microbial communities. Metagenomic analysis of coastal samples from the Gulf of Mexico (GOM) and British Columbia, Canada, revealed numerous sequences belonging to gokushoviruses and allowed the assembly of five putative genomes with an organization similar to chlamydiamicroviruses. Fragment recruitment to these genomes from different metagenomic data sets is consistent with gokushovirus genotypes being restricted to specific oceanic regions. Conservation among the assembled genomes allowed the design of degenerate primers that target an 800 bp fragment from the gene encoding the major capsid protein. Sequences could be amplified from coastal temperate and subtropical waters, but not from samples collected from the Arctic Ocean, or freshwater lakes. Phylogenetic analysis revealed that most sequences were distantly related to those from cultured representatives. Moreover, the sequences fell into at least seven distinct evolutionary groups, most of which were represented by one of the assembled metagenomes. Our results greatly expand the known sequence space for gokushoviruses, and reveal biogeographic separation and new evolutionary lineages of gokushoviruses in the oceans.

**Keywords: biogeography, ssDNA viruses,** *Microviridae***,** *Gokushovirinae***, virus diversity, ocean viruses**

## **INTRODUCTION**

Viruses are the most abundant (Suttle, 2005) and diverse (Breitbart et al., 2002; Angly et al., 2006) biological entities in the oceans. By causing lysis of specific subsets of microbial communities, they influence community composition by controlling species evenness and maintaining species richness (Hennes et al., 1995; Thingstad, 2000; Wommack and Colwell, 2000; Middelboe et al., 2001; Weinbauer, 2004; Winter et al., 2010); thereby, influencing nutrient and energy cycling (Fuhrman, 1999; Wilhelm and Suttle, 1999; Suttle, 2007). Moreover, viruses harbor an enormous pool of genetic diversity that can be exchanged among other viruses (Pedulla et al., 2003; Short and Suttle, 2005) and bacteria (Fuhrman and Schwalbach, 2003; Kenzaka et al., 2010). Despite the abundance of bacteriophages in marine systems (often *>*107 ml−1) and their important role in marine systems, relatively little is known about the distribution and composition of most groups of marine viruses.

Metagenomic approaches have provided an in-depth look at the molecular diversity of ssDNA viruses in a range of environments including marine systems (Breitbart et al., 2002; Angly et al., 2006; Bench et al., 2007), the human gut (Zhang et al., 2006; Breitbart et al., 2008; Minot et al., 2011), modern stromatolites (Desnues et al., 2008), and freshwaters (Kim et al., 2008; López-Bueno et al., 2009; Roux et al., 2012a). Recently, 608 ssDNA viral genomes were assembled from marine metagenomic data revealing far greater evolutionary diversity in ssDNA viruses than previously known (Labonté and Suttle, 2013).

Gokushoviruses are ssDNA bacteriophages belonging to the family *Microviridae* and are represented among sequences found in metagenomic data. For example, gokushovirus genomes were assembled from a wide range of environments by mining of metagenomic data, with 42 assembled from a variety of ecosystems (Roux et al., 2012b), and two others from data collected from the North Atlantic Ocean (Tucker et al., 2011), indicating the widespread occurrence of gokushoviruses. These viruses have a ∼30 nm icosahedral capsid encompassing a positive ssDNA molecule of 4.4 to 4.8 kb that encodes five major proteins. Based on the phylogeny of the major capsid protein (VP1) of isolates, the *Microviridae* are divided into two subfamiles (Brentlinger et al., 2002). Members of the *Microvirinae* (e.g., phiX174 and G4) infect enterobacteria including *Escherichia coli*(Godson et al., 1978), while members of the *Gokushovirinae* infect parasitic bacteria. The latter includes Chp1 (Storey et al., 1989), Chp2 (Liu et al., 2000; Everson et al., 2002) and Chp3 (Garner et al., 2004) that infect *Chlamydia* spp., while phiMH2K (Brentlinger et al., 2002) and SpV4 (Chipman et al., 1998) infect *Bdellovibrio* sp. and *Spiroplasma* sp., respectively.

There are no reported gokushovirus isolates, and their hosts remain unknown. Based on bacterial genomic sequences bacteria in the *Bacteroidetes* appear to be hosts for a proposed subfamily, *Alpavirinae* (Krupovic and Forterre, 2011), of previously unknown microviruses. As well, eight ssDNA phages have been isolated that are morphologically similar but genetically different to microviruses (Holmfeldt et al., 2013, 2012).

Our study examined the genetic diversity and relatedness of *Gokushovirinae*-like viruses from temperate and subtropical coastal environments. From three ssDNA-enriched metagenomic datasets we assembled and phylogenetically compared five new gokushovirus genomes. Recruitment of metagenomic reads to these genomes showed spatial differences in the most abundant gokushovirus genotypes. The genetic richness of gokushoviruses was also assessed through amplification of a ∼800 bp fragment of the conserved gene encoding the major capsid protein, VP1. These results reveal biogeographic separation and new evolutionary lineages of marine gokushoviruses, and likely reflect the underlying distributions of their hosts.

## **MATERIALS AND METHODS**

### **COLLECTION AND PREPARATION OF SAMPLES**

Samples (∼20 to ∼200 L) were collected using GO-FLO or Niskin bottles either mounted on a CTD rosette or directly on a hydrographic wire [Saanich Inlet (SI)], or by bucket from the surface (lake samples). For each sample, the viruses were concentrated ∼10–100-fold (∼200 mL final volume) using ultrafiltration (Suttle et al., 1991). Briefly, particulate matter was removed by pressure filtering (*<*17 kPa) the samples through 142-mm diameter glass-fiber filters (MFS GC50, nominal pore size 1.2 µm) and polyvinylidene difluoride filters (Millipore GVWP, pore size 0.22µm) connected in series. The viral size fraction in the filtrate was then concentrated by ultrafiltration though a 30 kDa molecular weight cut-off cartridge (Amicon S1Y30, Millipore), and stored at 4◦C in the dark until processed.

In order to integrate variation within a region, numerous virus concentrates (VCs) collected from different locations and at different times within a geographic region were combined into a single mix. Except for the SI and freshwater samples, these mixes corresponded to those used by Angly et al. (2006) in which the first ssDNA viral sequences were reported from marine metagenomic data. VCs from the Strait of Georgia (SOG) and surrounding inlets and bays were pooled into the following four mixes by combining 2 mL of each VC: BC1-1999 (23 samples), BC3-2000 (26 samples), BC4-2004 (16 samples), and BC2-Low salinity (19 samples). Similarly, samples from the Gulf of Mexico (GOM) were combined into four mixes from the eastern GOM (8 samples), northern GOM (6 samples), western GOM (6 samples), and the Texas Coast (13 samples), while samples from the Arctic Ocean were made into mixes from the Beaufort Sea (20 samples), Chukchi Sea (14 samples) and High Arctic (22 samples). To look at the diversity of freshwater gokushoviruses, two mixes were made from Chilliwack (6 samples) and Cultis (8 samples) Lakes. An extensive description of all the samples that were combined in each mix is presented in the Supplementary Material of Labonté and Suttle (2013). SI is unusual as it undergoes seasonal anoxia (Zaikova et al., 2010). For the metagenomic study, we combined surface samples from April 2007, and January, March, May, July, August, and November 2008. PCR amplifications were performed on the following nine samples from SI: 10, 120, and 150 m samples from April 2007, and surface samples from January, March, May, July, August, and November 2008.

## **ssDNA PREPARATION**

As described in Labonté and Suttle (2013), ssDNA was prepared from 10 mL of 0.22-µm filtered (PDVF; Millipore) pooled mixes from British Columbia (SOG), the GOM, the Lakes (LA), and the Arctic (ARC), or from 10 mL of each individual VCs from SI. Briefly, ssDNA was extracted using a silica column and amplified using multiple-displacement whole-genome amplification (WGA) to convert ssDNA into dsDNA. Pure amplified dsDNA was resuspended in ultrapure H2O for pyrosequencing or Tris-HCl for PCR amplification.

#### **GENOME ANALYSIS, BINNING, AND ASSEMBLY**

Metagenomic libraries were constructed from WGA ssDNA from SI, SOG, and GOM ssDNA. The purified WGA DNA was resuspended in 100µL of RNAse- and DNAse-free water (Invitrogen) and concentrated using a Millipore YM-30 Microcon centrifugal filter to a final volume of ∼50µL; 3–5 ug of DNA from each sample was sent for pyrosequencing (Roche 454 FLX instrumentation with Titanium chemistry) at Genome Québec, McGill University (SOG) and the Broad Institute at the Massachusetts Institute of Technology (GOM and SI).

The sequences were quality and linker trimmed, and assembled into contiguous sequences (contigs) using the Newbler Assembler (Roche). The individual reads and assembled sequences were compared to a database of all available genomes in GenBank (as of February 2010) from viruses belonging to the *Microviridae* using the tBLASTx algorithm with an *e*-value cut-off of 10−5. Reads with significant similarity to gokushoviruses were aligned onto the assembled contigs using the add454Reads.perl script and were reassembled into new contigs using the phredPhrap.perl script of the Consed package (Gordon, 2003). Additional contig analyses (BLAST, circularization of the genomes, annotations, alignments, and phylogeny) were performed within the Geneious Pro package v5.6 (Biomatters).

#### **PRIMER DESIGN AND PCR AMPLIFICATION**

Two forward (MicroVP1-F1, 5 -CGN GCN TAY AAY TTR ATH-3 ; MicroVP1-F2, 5 -AGN GCN TAY AAY TTR CTN-3 ) and two reverse (MicroVP1-R1, 5 -TTY GGN TAY CAR GAR AGN-3 ; MicroVP1-R2, 5 -NCT YTC YTG RTA NCC RAA-3 ) primers with respective degeneracies of 256, 215, 256, and 256 were designed from alignments of the inferred amino acid sequences of the major capsid protein (VP1) of the chlamydiaphages Chp1 (accession number NC\_001741.1), Chp2 (NC\_002194.1), Chp3 (NC\_008355.1), Chp4 (NC\_007461.1), phiCPAR39 (NC\_002180.1), and phiCPG1 (NC\_001998.1), Spiroplasmaphage Sp4 (NC\_003438.1), bdellovibriophage phiMH2K (NC\_002643.1), and the Sargasso Sea Chp1-like assembled genome (Angly et al., 2006; Tucker et al., 2011). The primers amplify a ∼800 bp VP1 gene fragment from the subfamily *Gokushovirinae* in the *Microviridae*.

Prior to use in PCR reactions, the purified WGA DNA was resuspended in 100µL of TE, and 10µL was used as a template in each PCR reaction mixture consisting of Taq DNA polymerase assay buffer [20 mM Tris·HCl (ph 8.4), 50 mM KCl], 1.5 mM MgCl2, 125µM of each deoxyribonucleoside triphosphate, 1µM of each MicroVP1-F1, MicroVP1-F2 and MicroVP1- R1 and MicroVP1-R2 primer and 2.5 U of PLATINUM Taq DNA polymerase (Invitrogen). Negative controls contained all reagents except DNA template. The samples were denatured at 94◦C for 3 min, followed by 35 cycles of denaturation at 94◦C for 30 s, annealing at 50◦C for 30 s, and elongation at 72◦C for 50 s, with a final elongation step of 72◦C for 5 min.

## **CLONE LIBRARY CONSTRUCTION AND RFLP ANALYSIS**

PCR amplicons were purified with a MinElute PCR purification kit (Qiagen), ligated into pCR2.1-Topo (Invitrogen), and used to transform chemically competent *E. coli* Top10 cells. For each sample, 30 clones were checked by colony PCR to verify that they contained an insert of the correct size. Restriction fragment length polymorphism (RFLP) analysis was then performed on 20 positive clones. For each RFLP reaction, 15µL of colony PCR product was digested with *AluI* (New England BioLabs) in a reaction containing 1 U/µg of DNA and 1× NEBuffer 4 (20 nM Tris-acetate, 50 mM potassium acetate, 10 mM magnesium acetate, 1 mM DTT, pH 7.9) by incubating at 37◦C for 16 h, followed by heat inactivation at 65◦C for 20 min. RFLP products were separated on a 2% agarose gel in 0.5× TBE (9 mM Tris base, 9 mM boric acid, 2 mM EDTA, pH 8.0) running at 110 V for ∼2 h. Sequencing of representative clones confirmed that each unique restriction pattern could be considered as an operational taxonomic unit (OTU). Forward and reverse sequences (∼800 bp) were obtained for each RFLP pattern using Big-Dye Terminator Cycle Sequencing (Applied Biosystems) and ABI 373 Stretch or ABI Prism 377 sequencers (Nucleic Acid Protein Service Unit, UBC).

## **PHYLOGENETIC ANALYSIS**

For the whole genome phylogeny, non-coding sequences were removed and the five major open reading frames were ordered. The sequences were aligned using MAFFT (Katoh et al., 2002) and maximum likelihood analysis with 100 bootstrap replicates were performed using PhyML (Guindon et al., 2010).

VP1 from the previously sequenced isolates, environmental sequences, and the degenerate PCR products from this study were trimmed to the PCR-product length (∼800 bp) and aligned using MAFFT (Katoh et al., 2002). The alignment was cured with GBlocks to remove unconserved regions that aligned with multiple gaps using the less stringent setting (allowing for smaller final blocks, gap positions within the final blocks and less strict flanking positions) (Talavera and Castresana, 2007). Bayesian phylogenetic analyses were performed on the cured alignment with MrBayes (Huelsenbeck and Ronquist, 2001). MrBayes uses a Markov chain Monte Carlo (mcmc) approach to approximate prior and posterior probabilities. Under the HKY85 substitution model with an invgamma distribution, two independent analyses of 4 (1 cold and 3 heated) mcmc chains with 20,000,000 cycles were run, sampled every 1000th cycle. The consensus tree was generated in Geneious with a burnin of 25%. Trees were viewed in Fig Tree (http://tree*.*bio*.*ed*.*ac*.*uk/software/figtree/).

## **FRAGMENT RECRUITMENT**

Recruitment of the reads from metagenomic data sets onto the assembled genomes was performed using tBLASTx with an *e*-value of 10−<sup>10</sup> and allowing only one hit per read. The metagenomic reads from the marine viromes (Angly et al., 2006) and microbialites (Desnues et al., 2008) were obtained from the CAMERA database, while the metagenomic reads from Lake Pavin and Lake Bourget (Roux et al., 2012a) were obtained from the SEED database. The environmental genomes used were Lake\_Bourget\_052, Lake\_Bourget\_523, Lake\_Pavin\_279 and 68\_Microbialite\_063 from Roux et al. (2012b), and SARssphi2 from Tucker et al. (2011).

### **NUCLEOTIDE SEQUENCE ACCESSION NUMBERS**

The five complete gokushovirus genomes as well as the 43 environmental PCR product sequences were submitted to Genbank and are available under the accession numbers KC131021- KC131025 and KC130978-KC131020, respectively.

## **RESULTS AND DISCUSSION**

## **ASSEMBLY OF COMPLETE GOKUSHOVIRUS GENOMES**

Sequence analysis of ssDNA metagenomic libraries from the SOG, SI, and the GOM recovered 1733, 374, and 194 sequences, respectively, that were significantly similar to sequences from viruses belonging to the *Microviridae*, with *>*90% of them being most similar to sequences belonging to the chlamydiamicroviruses and other gokushoviruses. From these data, five complete circular genomes were assembled with at least 3-fold coverage (two from SOG, two from SI and one from GOM). The genome sizes varied from 4062 to 5386 bp, and were uniformly shorter than those from previously sequenced isolates (**Figure 1**). Assembly of these genomes represented the accumulation of 95 reads for SOG-1, 58 for SOG-2, 53 for SI-1, 48 for SI-2, and 38 for GOM.

Even though there was only ∼30–50% similarity at the nucleotide level among the assembled genomes (**Table 1**), the chlamydiaphages and bdellovibriophage phiMH2K, the gene organization was remarkably similar among them, and included the five proteins required for replication of gokushoviruses (**Figure 1**), implying a common evolutionary origin. These comprise VP1, the major capsid protein, VP2 that is hypothesized to be involved in host recognition (Chipman et al., 1998) and virus attachment (Everson et al., 2003), VP3 that is a scaffolding protein found in the procapsid only and not in mature virions (Clarke et al., 2004), ORF4 that is a replication initiator involved in ssDNA synthesis (Liu et al., 2000; Garner et al., 2004; Salim et al., 2008), and ORF5 that is involved in DNA packaging (Liu

**FIGURE 1 | Gokushoviruses share a similar genome organization.** Whole genome phylogeny (Maximum likelihood, 100 bootstrap replicates, HKY85 model) on the ORFs of gokushoviruses rooted with the microvirus phiX174 (left) and pairwise comparisons of the five environmental gokushovirus genomes assembled from this study (bold) with the isolates and other

environmental genomes. Conserved genes are represented by colored arrows, while small overlapping genes of unknown function are represented by short black arrows. The genome similarities were visualized in ACT (Carver et al., 2005) (*e*-value *<*10−5) and the gray shading indicates the level of similarity; darker shading represents higher similarity between pairs of ORFs.

**Table 1 | Similarity matrix of the coding regions of the five environmental gokushovirus genomes assembled from this study (bold) with the isolates and other environmental genomes.**


*(i.e.* - *>90%,* - *45–89%,* -*30–45%, and <30%).*

et al., 2000; Garner et al., 2004; Salim et al., 2008). The presence of all five essential genes in the assembled genomes strongly suggests they represent complete sequences from extant viruses in the environment.

Whole genome phylogeny revealed that the environmental genomes cluster more closely with the bdellovibriophage phiMH2K, rather than the chlamydiaphages (**Figure 1**), suggesting that the host for these gokushoviruses is more closely related to *Bdellovibrio* spp., which are found in marine waters, than *Chlamydia* spp. Whole genome pairwise comparisons showed that VP2 and ORF4 are the least conserved genes, with very few regions of conservation. Moreover, there is 91–97% similarity among chlamydiamicroviruses, while only 28–49% similarity among the environmental phages. A recombination event in which ORF4 and ORF5 are inverted in phiMH2K, which infects the bacterial parasite *Bdellovibrio bacteriovorus*, and in the environmental genome SAR phi2. These genomes also cluster together suggesting a common evolutionary history (**Figure 1**).

All of the environmental genomes were shorter than those from isolates. Some, such as SI-1 and SOG-1, had multiple overlapping genes of unknown function. It is postulated (Rokyta and Burch, 2006) that ssDNA microviruses, such as the coliphages phiX174 and G4, evolve differently than dsDNA viruses because of strictly lytic life cycles, small genomes, and low rates of horizontal gene transfer (Breitbart and Rohwer, 2005; Comeau and Buenaventura, 2005; Hambly and Suttle, 2005). Novel genes were predicted to originate by overprinting rather than by horizontal gene transfer (Pavesi, 2006).

## **GENETIC RELATEDNESS AMONG GENES ENCODING THE MAJOR CAPSID PROTEIN**

To look more deeply at the genetic richness of gokushoviruses, degenerate primers were designed to amplify a ∼800 bp fragment of the gene encoding the major capsid protein (VP1) that has interspaced conserved and variable regions. For the assembled genomes, the phylogeny of VP1 is congruent with the whole genome; thus, the phylogeny of VP1 can be used to infer viral phylogeny. PCR amplification was performed on samples from the SOG (4 mixes), GOM (4 mixes), SI (9 samples), and Arctic (ARC; 3 mixes) (Table S1). No products were amplified from the ARC, LA, SOG-Low Salinity, or Eastern GOM mixes. This means that gokushoviruses were absent or at low concentrations in these samples, or that they are too divergent to be amplified by the primers.

Twenty VP1 clones from each of the 20 samples were digested using *AluI* to reveal 77 different RFLP patterns. Representative clones from each restriction pattern were sequenced (data not shown). Of these, 43 sequences were at least 98% different at the nucleotide level, and thus identified as unique gokushovirus VP1 sequences. Some sequences occurred in more than one sample from the same geographic region; for example, the sequence SI-07 was sequenced from multiple dates in SI (i.e., January, April, May and Aug), but no sequences occurred in more than one location (**Figure 2**). Some sequences were found in both SOG and SI, but no sequences were found in GOM and SOG, or GOM and SI. The 85 % of the VP1 sequences that contained the primer sequences and translated into a putative protein were kept for further analysis. The nucleotide alignment revealed multiple regions of conservation, as well as regions that were confined to specific groups, agreeing with observations made during primer design.

The amplified VP1 PCR products were compared to VP1 sequences containing the primer sequences from the assembled genomes, the chlamydiamicrovirus isolates, as well as Genbank environmental sequences from modern stromatolites (Desnues et al., 2008), freshwater (Lake Needwood, MD; Kuzmickas et al., unpublished), rice paddy soil (Kim et al., 2008) and marine genomes (Venter et al., 2004). Few sequences were similar to those from isolates. Phylogenetic analysis using Bayesian (**Figure 2**) and maximum-likelihood algorithms produced similar trees with gokushoviruses sub-divided into at least seven well-supported new groups containing more than two sequences (**Figure 2**). Five of the new clades are represented by an assembled genome or sequenced isolate (**Figure 2**). Several sequences, such as SOG4- 04 and SI-18, were too divergent to be assigned to a cluster; however, as only 20 clones were analyzed for each sample, rarer phylogenetic clusters were poorly sampled.

Sequences from a given location were usually more closely related to ones from the same location; most GOM sequences clustered within ENV6 and ENV7, while ENV2 is represented exclusively by SOG sequences. Sequences found in more than one sample also usually clustered together. For example, the sequences SI-10 and SI-07, which were found in SI on multiple dates, clustered within ENV5, along with sequences from GOM and SI-2. Collectively, these data imply that viruses in the ENV5 group are widespread in nature. Other data from modern stomatolites and marine genomes clustered together as specific phylogenetic groups.

## **HOST SPECIFICITY AND GEOGRAPHIC DISTRIBUTION OF ENVIRONMENTAL GOKUSHOVIRUS GENOMES**

Isolates in the *Gokushovirinae* infect parasitic bacteria, such as *Chlamydia* spp., *Bdellovibrio* spp., and *Spiroplasma* spp., with host specificity likely being dictated by variable genomic regions. To investigate conserved and variable motifs, metagenomic reads from our ssDNA data, as well as other viral metagenomic data sets from marine (Angly et al., 2006), freshwater (Roux et al., 2012a), and microbialite (Desnues et al., 2008) environments were recruited against environmental gokushovirus genomes (**Figure 3**). Recruitment was more even when the reads were recruited against genomes assembled from metagenomic data collected from the same region (**Figures 3**, **4**; **Figures S1**, **S2**). For the SOG-1, SOG-2, and SI-2 genomes, few reads were recruited from data sets other than those from which the genomes were assembled, suggesting that these genomes are not widespread (**Figure 3**). In contrast, reads from all of the metagenomic data sets aligned on the GOM genome (**Figure 3**), indicating a wider geographic distribution of these viruses. The high level of recruitment from other data sets on the GOM genome is also congruent with the phylogenetic clustering of the VP1 gene with other VP1 sequences that were present in multiple samples (**Figure 2**).

The distribution of reads from other environmental samples that recruited to the assembled genomes was very uneven, showing regions of higher conservation within VP1 and VP3, while few reads were recruited to the VP2 region, indicating high variability in this gene (**Figure 4**). Since the VP2 sequences of the assembled genomes differ from those of isolates, and VP2 encodes for the minor capsid protein involved in host recognition (Chipman et al., 1998), the environmental sequences are likely not from viruses infecting the genera *Chlamydia*, *Bdellovibrio*, or *Spiroplasma*. Recruitment to ORF4 was limited to the source environment for the assembled genomes, and metagenomic data from British Columbia were not recruited to ORF4 of the GOM genome. Thus, both the pairwise comparison (**Figure 1**) and the recruitment of metagenomic reads (**Figure 3**) showed that VP2 and ORF4 are less conserved. Similar patterns were observed with genomes assembled from other data sets (**Figures S1** ), showing that specific gokushovirus genotypes are restricted in distribution. , **S2**

No sequences similar to gokushoviruses were amplified from the Arctic Ocean or two lakes in British Columbia. However,

gokushovirus sequences have been found in an Antarctic lake (López-Bueno et al., 2009) and other freshwater environments (Roux et al., 2012a), suggesting that freshwater gokushoviruses differ enough in sequence that they cannot be amplified using our primers.

#### **POSSIBLE ROLE OF GOKUSHOVIRUSES IN AQUATIC ENVIRONMENTS**

The distribution of gokushovirus OTUs with respect to specific marine environments differs from observations that some viral genotypes are widely distributed Chen and Suttle, 1996; Fuller et al., 1998; Hambly et al., 2001; Short and Suttle, 2002, 2005; Breitbart et al., 2004; Labonté et al., 2009. However, Tucker et al. (2011) observed differences in the depth distribution of gokushovirus sequences in the North Atlantic Ocean that likely reflected the distribution of hosts. Most VP1 sequences that were found in more than one sample were also relatively closely related, perhaps reflecting viruses that have a broader host range, or viruses that infect widely distributed hosts. In contrast, sequences specific to a single location are probably from viruses that infect bacteria that are environment specific. Although some bacterial species are very widely distributed (Rusch et al., 2007; Biers et al., 2009), others are restricted to specific habitats (Biers et al., 2009). Hence, it is not surprising that some gokushoviruses have a very restricted distribution.

**the genomes from this study to show regions of conservation among ssDNA gokushoviruses.** Each assembled genome (GOM, SI1, SI2, SOG1, SOG2) is represented by a different panel. Each horizontal line represents a metagenomic read from ssDNA data

and Strait of Georgia (aqua) on each of the assembled genomes. Reads were recruited using tBLASTx with an *e*-value of 10−10. The position of each line represents the percent similarity of the read to the genome.

**FIGURE 4 | Fragment recruitment of reads from environmental viral metagenomes to show the regions of conservation within different environments.** Each assembled genome (GOM, SI1, SI2, SOG1, SOG2) is represented by a different panel. Each horizontal line represents a read recruited from one of the following publicly available metagenomic data sets: Gulf of Mexico (dark red), Strait of Georgia (aqua), Sargasso Sea (orange), Lake Bourget (light green), Lake Pavin (dark green), and microbialites (purple). Reads were recruited against each of the assembled genomes using tBLASTx with an *e*-value of 10−10. The position of each line represents the percent similarity of the read to the genome. VP1 is represented by a red arrow.

Based on previous work (Labonté and Suttle, 2013), gokushovirus sequences were not the most abundant viruses in our samples, and comprised only 1.6, 0.4, and 0.2% of the metagenomic reads from the SOG, SI, and GOM, respectively. In contrast, in metagenomic data from the Sargasso Sea, gokushovirus sequences comprised nearly 6% of the reads (Angly et al., 2006), while in Lake Bourget they were more than 90% of the sequences (Roux et al., 2012a). These differences may be because the small genomes of gokushoviruses permit rapid replication and high burst sizes, and allow them to dominate following a lytic event, this is consistent with the hypothesis that the most abundant marine viruses are virulent opportunists that replicate rapidly, have high burst sizes and small genomes in order to exploit rapidly growing populations of rare marine bacteria such as *Roseobacter* spp. or *Vibrio* spp. (Suttle, 2007). For example, ∼500 genomes are produced each time the chlamydiaphage Chp2 infects its parasitic host (Salim et al., 2008). High burst sizes coupled with genomes usually *<*5 kb support the idea that gokushoviruses are highly virulent and are selected for rapid population growth, which are characteristics of r-strategists. In contrast, many large DNA viruses have a low burst size, large genome and decay slowly, which are characteristics of K-strategists.

Discovering the hosts of marine gokushoviruses is a high priority in order to understand the roles that these viruses play in ecosystems. Given the challenges in culturing marine microbes, culture-independent techniques will likely be needed to determine the hosts for most of these viruses. One approach that we have tried with some success is to use fluorescence *in situ* hybridization (FISH) using labeled VP1 sequences to probe natural microbial communities. Another approach that has been used to visualize phage-infected gammaproteobacterial cells is phage-FISH (Allers et al., 2013), which could be adapted to search for cells infected by gokushoviruses. Finally, single-cell genomics (SCG) allows everything in a cell, including plasmids and viruses to be sequenced (Stepanauskas, 2012). If applied to samples with abundant gokushoviruseses, it should be possible to sequence infected cells.

This manuscript presents a new set of degenerate primers that have been used to reveal at least five new evolutionary groups of gokushoviruses, and clearly show they share a common evolutionary history with viruses that infect the obligate intracellular parasitic bacteria *Chlamydia* and *Bdellovibrio*. Phylogenetic analysis of the major capsid protein, combined with fragment recruitment of environmental metagenomic sequences shows that the distribution of some evolutionary groups of gokushoviruses is very environment dependent, whereas others are more cosmopolitan. The high-burst size, rapid replication rates and likely lytic nature of these viruses suggests that they may play an important role as mortality agents in marine systems.

### **ACKNOWLEDGMENTS**

We thank members of the Suttle laboratory for collecting and processing the samples, and the Hallam laboratory for providing filtered water from Saanich Inlet that made this study possible. This research was supported by the Natural Science and Engineering Research Council of Canada (NSERC) through a postgraduate scholarship (Jessica M. Labonté) and Discovery grants (Curtis A. Suttle). Sample collection was facilitated through ship-time grants from NSERC that supported sample collections from the Strait of Georgia (Curtis A. Suttle) and Saanich Inlet (PD Tortell and SJ Hallam), the US National Science Foundation (Gulf of Mexico), and through the Canadian Arctic Shelf Exchange Study (NSERC) and the Japan/Canada Western Arctic Climate Study. Access to sequencing was funded by the Gordon and Betty Moore Foundation through GBMF1799 to the Broad Institute, and by NSERC and the Tula Foundation using facilities at the McGill University and Genome Quebec Innovation Centre.

## **SUPPLEMENTARY MATERIAL**

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

**Figure S1 | Fragment recruitment of viral ssDNA reads onto assembled environmental gokushovirus genomes.** Each assembled genome (Lake\_Bourget\_052, Lake\_Bourget\_523, Lake\_Pavin\_279, 68\_Microbialite\_063, and SARssphi2) is represented by a different panel. Each horizontal line represents a read recruited from one of the metagenomic data sets from this study: Gulf of Mexico (dark red), Saanich Inlet (Dark blue), and Strait of Georgia (aqua). Reads were recruited against each of the assembled genomes using tBLASTx with an *e*-value of 10−10. The position of each line represents the percent similarity of the read to the genome. VP1 is represented by a red arrow.

## **Figure S2 | Fragment recruitment of reads from environmental viral metagenomes to show the regions of conservation within different environments.** Each assembled genome (Lake\_Bourget\_052,

Lake\_Bourget\_523, Lake\_Pavin\_279, 68\_Microbialite\_063, and SARssphi2) is represented by a different panel. Each horizontal line represents a read recruited from one of the following publicly available metagenomic data sets: Gulf of Mexico (dark red), Strait of Georgia (aqua), Sargasso Sea (orange), Lake Bourget (light green), Lake Pavin (dark green), and microbialites (purple) metagenomic data sets on each of the assembled genomes that recruited at using tBLASTx with an *e*-value of 10−10. The height of line represent the percent similarity of the read to the genome. VP1 is represented by a red arrow.

## **REFERENCES**


difficult to detect with DNA-binding stains. *Appl. Env. Microb.* 78, 892–894. doi: 10.1128/AEM.06580-11


in nature. *Appl. Env. Microb.* 68, 1290–1296. doi: 10.1128/AEM.68.3.1290- 1296.2002


organisms: the "killing the winner" hypothesis revisited. *Microbiol. Mol. Biol. R* 74, 42–57. doi: 10.1128/MMBR.00034-09


**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.

*Received: 12 October 2013; accepted: 06 December 2013; published online: 24 December 2013.*

*Citation: Labonté JM and Suttle CA (2013) Metagenomic and whole-genome analysis reveals new lineages of gokushoviruses and biogeographic separation in the sea. Front. Microbiol. 4:404. doi: 10.3389/fmicb.2013.00404*

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

*Copyright © 2013 Labonté and Suttle. 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.*

## Unprecedented evidence for high viral abundance and lytic activity in coral reef waters of the South Pacific Ocean

## *Jérôme P. Payet1,2 \*, Ryan McMinds1, Deron E. Burkepile3 and Rebecca L. Vega Thurber1*

<sup>1</sup> Department of Microbiology, Oregon State University, Corvallis, OR, USA

<sup>2</sup> Institute for Pacific Coral Reefs, Moorea, French Polynesia

<sup>3</sup> Department of Biological Sciences, Florida International University, Miami, FL, USA

#### *Edited by:*

Stephen Tobias Abedon, The Ohio State University, USA

#### *Reviewed by:*

Marla Tuffin, University of the Western Cape, South Africa Laurent Seuront, Centre National de la Recherche Scientifique, France

#### *\*Correspondence:*

Jérôme P. Payet, Department of Microbiology, Oregon State University, 226 Nash Hall, Corvallis, OR 97331, USA e-mail: payetj@onid.orst.edu

Despite nutrient-depleted conditions, coral reef waters harbor abundant and diverse microbes; as major agents of microbial mortality, viruses are likely to influence microbial processes in these ecosystems. However, little is known about marine viruses in these rapidly changing ecosystems. Here we examined spatial and short-term temporal variability in marine viral abundance (VA) and viral lytic activity across various reef habitats surrounding Moorea Island (French Polynesia) in the South Pacific. Water samples were collected along four regional cross-reef transects and during a time-series in Opunohu Bay. Results revealed high VA (range: 5.6 <sup>×</sup> 106–3.6 <sup>×</sup> <sup>10</sup><sup>7</sup> viruses ml−1) and lytic viral production (range: 1.5 <sup>×</sup> 109–9.2 <sup>×</sup> 1010 viruses l−<sup>1</sup> <sup>d</sup>−1). Flow cytometry revealed that viral assemblages were composed of three subsets that each displayed distinct spatiotemporal relationships with nutrient concentrations and autotrophic and heterotrophic microbial abundances. The results highlight dynamic shifts in viral community structure and imply that each of these three subsets is ecologically important and likely to infect distinct microbial hosts in reef waters. Based on viral-reduction approach, we estimate that lytic viruses were responsible for the removal of ca. 24–367% of bacterial standing stock d−<sup>1</sup> and the release of ca. 1.0– 62 g of organic carbon l−1 1 μ d<sup>−</sup> in reef waters. Overall, this work demonstrates the highly dynamic distribution of viruses and their critical roles in controlling microbial mortality and nutrient cycling in coral reef water ecosystems.

**Keywords: marine viruses, viral lysis, carbon cycling, coral reefs, South Pacific, microbial mortality, viral abundance, spatial and temporal variability**

## **INTRODUCTION**

Viruses are increasingly recognized as the most abundant and dynamic biological entities in marine ecosystems (e.g., reviewed in Fuhrman, 1999; Wommack and Colwell, 2000; Weinbauer, 2004; Suttle, 2007). Viral-mediated cell lysis can cause significant mortality of heterotrophic bacteria, cyanobacteria and eukaryotic phytoplankton (Wilhelm and Suttle, 1999; Brussaard, 2004b). Models and empirical studies have estimated that 20–50% of marine microbial communities are infected by viruses each day (e.g., reviewed in Fuhrman, 1999;Wilhelm and Suttle, 1999; Suttle, 2005, 2007). The release of organic cellular content and nutrients upon viral lysis can stimulate autotrophic and heterotrophic microbial activity (Gobler et al., 1997; Middelboe and Lyck, 2002; Weinbauer et al., 2011; Shelford et al., 2012) and increase diversity (Weinbauer and Rassoulzadegan, 2004; Motegi et al., 2013), with major effects on global biogeochemical cycles and flow of energy in the oceans (Fuhrman, 1999; Wilhelm and Suttle, 1999; Suttle, 2007). Despite their critical impact in the oceans, there is still a lack of data on the spatial and temporal dynamics of viruses and their ecological influence in marine microbial communities.

Tropical coral reefs are highly productive and diverse ecosystems yet thrive under oligotrophic conditions. Accumulating evidence also suggests that reef waters harbor abundant and active microbial communities (Moriarty et al., 1985; Gast et al., 1998; Charpy et al., 2012 and references therein) that can respond rapidly to changes in environmental conditions (Haas et al., 2011, 2013; Nelson et al., 2011; Mccliment et al., 2012). In this setting of nutrient poor conditions and high microbial abundance, viruses may play a particularly important ecological role in shaping microbial communities, with potential impacts on carbon cycling and energy transfer to higher trophic levels.

The spatial and temporal patterns of viral abundance (VA) and production have been relatively well studied in various marine environments over the past decades, but only a few studies have focused on marine viruses in tropical and subtropical reef waters (Paul et al., 1993; Seymour et al., 2005; Patten et al., 2006; Bouvy et al., 2012). These few studies have suggested that viruses are as highly dynamic and abundant as reported in higher latitude marine environments (e.g., reviewed in Vega Thurber and Correa, 2011). To our knowledge, the only study that has investigated lytic viral activity in coral reef waters estimated that lytic viruses were not a significant source of mortality for bacteria in atoll reef waters (Bouvy et al., 2012), in contrast to general findings from other marine ecosystems. More work is needed to fully elucidate the potential ecological roles of viruses in coral reef waters.

Here, we investigated viral abundance, subset composition as detected through flow cytometry (FC), and production in the south Pacific island of Moorea, French Polynesia. Particularly, we evaluated whether VA, structure, and lytic activity changed across distinct reef habitats and time. Furthermore, we used multivariate analysis to assess potential ecological factors controlling distribution patterns of VA and lytic activity; specifically, we asked whether patterns were driven by changes in trophic status of the ecosystem and/or by environmental conditions. Finally, we aimed to assess how virus-mediated mortality of heterotrophic bacteria can influence dissolved organic carbon (DOC) availability in oligotrophic reef habitats. Collectively, this novel dataset allowed us to determine whether viruses are dynamic and important players in tropical planktonic reef ecosystems.

## **MATERIALS AND METHODS**

#### **STUDY AREA**

This study was conducted at Moorea Island, in French Polynesia, in the South Pacific Ocean (**Figure 1**), during the dry season in August 2013. Moorea is a high basaltic island surrounded by barrier reefs that extend between 500 and 1500 m offshore, creating semi-enclosed back reef lagoons (e.g., reviewed in Leichter et al., 2013; **Figure 1**). For consistency herein, the semi-enclosed lagoons of individual reef platforms will be referred to as "lagoon". Eleven passes connect the lagoons to the open ocean, with some continuing near-shore as narrow deep channels (10–20 m width, 10–30 m depth). The typical reef zonation includes a fringing reef (FGR) nearest to shore (10–100 m width, <1 m depth), a shallow lagoon

(100–1000 m width, 1–6 m depth) interrupted with the occasional along-shore channel, a back reef (100–200 m width, 1–3 m depth), a reef crest (10–50 m width, <1 m depth) and an oceanward fore reef on a high downslope (50–200 m width, 2–60 m depth). On the north shore, the lagoon is connected to two narrow, 3 km long straight water bays [Opunohu Bay (OB) and Cook's Bay, <90 m deep]. Both bays are influenced by small river discharges that peak during rainfall; OB is also influenced by runoff from a nearby agricultural area that includes farming of prawns (Wolanski and Delesalle, 1995; Hench et al., 2008).

#### **SEAWATER SAMPLING**

Two types of sampling were conducted to examine spatiotemporal variability in biotic (e.g., viruses, heterotrophic, and autotrophic microbes) and abiotic (e.g., nutrients, longitude, latitude) variables. For the spatial study, samples were collected separately along four cross-reef transects located in four different geographic regions (i.e., north, east, west, and south) surrounding Moorea. For each transect, seawater was collected at four different reef habitats starting offshore and moving toward inshore: reef crest (CR), back reef (BR), lagoon (LAG), and FGR (**Figure 1**). An additional site located in the fore reef (FOR) was also collected for the north transect (**Figure 1**). FOR sites from other transects could not be sampled due to logistical constraints. For the short-term temporal study, a site located in OB, was sampled every 2–3 days from 8 to 27 August 2013.

**FIGURE 1 | Map of the study area.** Spatial samples were collected along four regional reef transects (east, west, south, and north; yellow lines). Each transect crossed distinct reef habitats: fringing reef (FGR), lagoon (LAG), back reef (BR), crest (CR), and fore reef (FOR); yellow

dots. Note the FOR was sampled only on the north transect. Temporal samples were collected every 2–3 days for 21 days in a site in the Opunohu Bay (pink dot). Aerial views of the north and south transects are shown.

**FIGURE 2 | Examples of typical flow cytograms of (A) viruses, (B) heterotrophic bacteria, and (C) phytoplankton.** Three viral subsets, with low, medium, and high nucleic acid fluorescence (V1, V2 and V3, respectively), two heterotrophic bacterial subsets with high and low nucleic acid fluorescence (HNA and LNA, respectively) and four phytoplankton subsets (Prochl, Prochlorococcus; Syn, Synechococcus; APP, autotrophic picoplankton; ANP, autotrophic nanoplankton) were distinguished using flow cytometry.

Seawater samples (0.5 l) for nutrient concentration, microbial and viral abundances (see below) were collected for all sites using high-density polyethylene bottles at ca. 0.5 m below the surface. Additional seawater samples (4 l) were collected for viral production (VP) assays (see below) using 4-l low-density polyethylene collapsible Cubitainers®. All bottles and cubitainers were acid washed (∼10% HCl), rinsed with MilliQ® water and then rinsed several times with *in situ* seawater before collection. All samples were transported back to the onshore laboratory and processed in under 2 h following collection. Samples were consistently collected between 9:00 and 12:00 h to avoid diel variation.

**Table 1 | Equations used in estimating viral production (VP), viral turnover (VT), viral-mediated mortality of bacteria (VMB), amount of organic carbon released upon viral lysis (OCr) and percent of bacterial standing stock removed due to viral lysis (BSSr).**


VA, viral abundance; to, initial time point of incubation; tf, final time point of incubation; BA, bacterial abundance; BAa, ambient bacterial abundance; BS, burst size.

#### **PREPARATION OF SAMPLES FOR FLOW CYTOMETRY**

Duplicate aliquots (1.8 ml) of each sample were dispensed into 2 ml-cryotubes containing gluteraldehyde (0.5% final concentration, electron microscopy grade, Sigma-Aldrich). Samples were fixed at 4◦C for 30 min then immediately frozen at −80◦C, before being analyzed using FC within 4–6 weeks (see below). Samples were shipped frozen on dry-ice back to Oregon State University (OSU). Due to logistic constraints, the samples could not be flashfrozen in liquid nitrogen before being stored at −80◦C. This is known to account for some virus losses, however they are reported to be minimal (<10%; Brussaard, 2004a). We therefore expect that virus loss in our samples was minimal, and that our data represent conservative estimates of VA. FC analysis of viruses, heterotrophic bacteria and phytoplankton were performed on a Becton Dickinson (BD) FACSCalibur flow cytometer (15 mW argon laser exciting at 488 nm, BD, San Jose, CA, USA), as described below.

#### **ENUMERATION OF VIRUSES AND HETEROTROPHIC BACTERIA**

Viruses and heterotrophic bacteria were enumerated separately according to standard protocols outlined in Payet and Suttle (2008) and Brussaard et al. (2010), respectively. Viruses and heterotrophic bacteria were discriminated by their signals in a bivariate scatter plot of side scatter (SSC) vs. green fluorescence (FL1, 530/30 nm), using FL1 as the threshold trigger. At least three viral subgroups were discriminated based on their relative SYBR green I fluorescence (V1, V2, and V3, respectively; **Figure 2A**). The total VA (VA) presented in this study is the sum of V1, V2, and V3. FC allowed separation and enumeration of a high nucleic acid (HNA) containing bacteria and a low nucleic acid (LNA) containing bacteria on the basis of their SSC vs. FL1 signals (**Figure 2B**; Gasol et al., 1999; Lebaron et al., 2001). Total heterotrophic bacterial abundance (BA) was calculated as the sum of HNA and LNA cells.

#### **ENUMERATION OF PHYTOPLANKTON**

Phytoplankton were enumerated using FC, following standard procedures (Olson et al., 1993; Marie et al., 2001). Just before the analysis, a mixture of yellow–green fluorescent 0.92 and 3 μm beads were added to the samples (ca. 105 beads ml−<sup>1</sup> final concentration) for internal standard. The threshold trigger was set to FL3. Phytoplankton populations were differentiated based on SSC, chlorophyll fluorescence (FL3) and phycoerithrin fluorescence (FL2, 585/42 nm) signals. In this study, FC differentiated autotrophic pico- (<2 μm) and nanoplankton (2– 20 μm); hereafter referred to as autotrophic picoplankton and autotrophic nanoplankton (APP and ANP, respectively), as well as the picocyanobacteria *Synechococcus* and *Prochlorococcus* (hereafter referred to as *Syn* and *Prochl*, respectively; **Figure 2C**). Total phytoplankton abundance (PA) was calculated as the sum of APP, ANP, *Syn,* and *Prochl*.

#### **MEASUREMENTS OF AMBIENT SEAWATER ABIOTIC VARIABLES**

Samples (90 ml) were filtered through combusted GF/C filters (Whatman GF/C, 25 mm diameter, 0.45 μm pore size) for nutrient analyses. Filters were held using acid-cleaned polycarbonate filter holders. Filter holders were attached directly to the outlet of acid-cleaned 60 cc syringes. For each sample, seawater filtrates were collected into duplicate acid-cleaned 30 ml HDPE bottles and stored upright at −80◦C until analysis at OSU within 2 months. Concentrations of dissolved inorganic nitrate plus nitrite (N+N), ammonium (NH4), and soluble reactive phosphorus (SRP) were measured using a hybrid airsegmented flow system consisting of a Technicon AutoAnalyzer II (SEAL Analytical Ltd., Milwaukee, WI, USA) and an Alpkem Rapid Flow Analyzer (Alpkem Series 300, Corp., Clackamas, OR, USA) following standard colorimetric protocols adapted from Gordon et al. (1993). In this study, we define dissolved inorganic nitrogen (DIN) as the sum of N+N and NH4 concentrations.

#### **MEASUREMENTS OF LYTIC VIRAL PRODUCTION AND ACTIVITY**

Lytic VP assays were carried out using the 4 l seawater samples collected (see above) from the FRG, BR, and CR within

habitats: the fringing reef (FRG), lagoon (LAG), back reef (BR), crest (CR), and fore reef (FOR). Note the FOR was sampled only on the north transect. Black dots indicate the relative position of the samples collected. Contour plots indicate the mean values of duplicate samples.

each transect, as well as in the OB for all the sampling dates. We used the viral-reduction approach of Winget et al. (2005) adapted from Wilhelm et al. (2002). Briefly, 900 ml of seawater was filtered through 20 μm mesh-size Nitex® screen to remove large particles. Filtered sample was then reduced to ca. 100 ml using a 0.22 μm pore-size polysulfone (PES) membrane tangential flow filter (TFF, GE Healthcare, Life Sciences). This process reduces particles <0.22 μm in diameter (i.e., most viruses infecting prokaryotes) while retaining particles ranging in size between 0.22 and 20 μm (i.e., pro- and eukaryotic microbes). The resulting retentate was subsequently washed with 900 ml of ultrafiltered (UF) seawater (<100 kDa cutoff, PES membrane TFF, GE Healthcare, Life Sciences) made from the same original seawater to further reduce viral-size particles. When ca. 100 ml of retentate

Circles represent mean values of duplicate samples collected for each

remained the sample was brought back to its original volume to produce a virus-reduced sample (i.e., 900 ml). All the TFF cartridges and tubing were cleaned with NaOH 0.1N and thoroughly rinsed with MilliQ® water and UF seawater before use. On average, the viral-reduction approach removed 65 ± 5% (range: 22–92%) of *in situ* VA and 57 ± 9% (range: 16–81%) of *in situ* BA, respectively. The same flow rates and processing times were used in all experiments.

transect.

The resulting virus-reduced sample was dispensed into triplicate sterile 50 ml conical tubes (BD Flacon) before incubation at *in situ* temperature (26 ± 1◦C) for 12–18 h in a temperaturecontrolled room in the dark. Samples (1 ml) for determination of VA and BA were collected every 3–4 h. For each individual incubation, VP was estimated from the slope of a least-square linear regression fitted to VA increases over time after correcting for *in situ* BA losses during the filtration, as described inWilhelm et al. (2002).

Viral turnover rates (VT, d−1), viral-induced mortality of bacteria (VMB, bacteria l−<sup>1</sup> d−1), percentage of bacterial standing stock removed (%BSSr, d−1) and extracellular dissolved organic carbon released (OCr, μgCl−<sup>1</sup> d−1) were calculated as in Wilhelm et al. (2002) and Payet and Suttle (2013; **Table 1**). We used a burst size (BS) of 30 viruses per lytic event, which was close to the average BS estimates of 28 reported in South Pacific tropical waters (Bouvy et al., 2012) and of 24 reported for marine environments (Parada et al., 2006). A cellular carbon quota of 20 fg C per marine bacterium was used to convert BA into organic carbon units (Lee and Fuhrman, 1987).

#### **DATA ANALYSIS AND STATISTICS**

Differences among mean biotic/abiotic variables during timeseries and across reef-transects were tested by Kruskal–Wallis (KW) analysis of variance on ranks, as the data did not meet the assumptions of normal distribution and homoscedasticity needed for analysis of variance tests (Zar, 1999). When KW tests were significant, the Dunn's *post hoc* test was performed to evaluate within group differences.

The distance-based linear model (DistLM) analysis (Legendre and Anderson, 1999; Anderson, 2004) was carried out to examine which biotic or abiotic variables were potential predictors of spatiotemporal variations in the viral variables (i.e., VA, V1, V2, V3, and VP). For this analysis, Bray–Curtis dissimilarity matrices of log-transformed datafor a selected viral variable were fitted against the abiotic (see below) and biotic (i.e., log-transformed BA, HNA, LNA, PA,*Prochl*, *Syn*, ANP, and APP) variables. A forward selection procedure based on Akaike's Information Criterion with a secondorder bias correction for small sample size (AICc) measures of fit was used to determine which explanatory variables could best predict selected viral variables (see Burnham and Anderson, 2004).

Highly correlated explanatory variables (*r* > 0.9) were omitted for the DistLM procedure. *P*-values were obtained using 999 random permutations of the data. For the spatial dataset, the abiotic variables included nutrient concentrations (log-transformed DIN and SRP) and coordinates (latitude and longitude). For the temporal dataset, the abiotic variables included log-transformed nutrients and time (number of days after first sampling). All the abiotic variables were normalized prior to DistLM procedure. Statistical analyses were performed using RStudio Version 0.97.551 (http://www.rstudio.org/; Racine, 2012) and PRIMER 6 with the PERMANOVA+ add-on (PRIMER-E, Plymouth, UK; Clarke and Gorley, 2006; Anderson et al., 2008). Means ± standard deviations (SD) are reported in the text for specific data sets.

## **RESULTS AND DISCUSSION**

This study examined spatial and short-term temporal variability of VA and lytic activity in relation to changes in microbial population abundances and environmental conditions, for the first time, in coral reefs surrounding Moorea Island, in the South Pacific Ocean. The results show high spatial heterogeneity and relatively low temporal changes in VA and lytic activity, concomitant with shifts in microbial host population dynamics. Overall, our data suggest that viral-induced lysis can exert

strong controlling influences on heterotrophic BA, with implications for nutrient and carbon fluxes in these oligotrophic ecosystems.

## **ENVIRONMENTAL CONDITIONS**

Relatively low mean nutrient concentrations measured at all sites confirmed the oligotrophic nature of this reef ecosystem, with SRP and DIN averaging 0.35 ± 0.08 μM and 0.49 ± 0.24 μM, respectively (data not shown). Although no significant differences in nutrient concentrations were detected among reef habitats, DIN concentrations were higher in the FRG (0.48 ± 0.14 μM) relative to the LAG (0.40 ± 0.13 μM) and BR (0.38 ± 0.09 μM), CR (0.34 ± 0.09 μM), and FOR (0.41 ± 0.10 μM; KW, *p* > 0.05). In the OB, DIN concentrations (0.93 ± 0.15 μM) were ∼1.9-fold higher, but not significantly different, than in the FRG (KW, *p* > 0.05). During the time-series, DIN concentrations remained relatively stable in the OB with low date-to-date variations (range: 0.7- to 1.2-fold). No significant differences were detected between sampling dates (KW, *p* > 0.05).

## **HETEROTROPHIC BACTERIA DISPLAY SPATIAL VARIABILITY AND SHORT-TERM TEMPORAL STABILITY**

Bacterial abundance averaged 2.7×10<sup>5</sup> <sup>±</sup>1.5×10<sup>5</sup> cells ml−<sup>1</sup> and did not vary significantly among the four transects (KW, *p* > 0.05;

(V3). See **Figure 3** for legend.

**Figures 3A–D**). However, consistent spatial trends emerged within transects, according to reef habitat (**Figure 4A**). On average, BA decreased 1.5-fold from the FOR toward the CR and a ∼2-fold from the CR toward the BR and LAG. Similar decreasing trends in BA from the FOR toward the BR were reported during a long-term study in Moorea (Nelson et al., 2011). The authors hypothesized that low wave-driven circulation and long water turnover time in the BR could increase encounter rates between bacteria and heterotrophic benthic organisms, resulting in low abundances. BA increased ∼3-fold in the FRG relative to the BR and LAG (KW with Dunn's test, *p* < 0.05; **Figure 4A**). Given that DIN concentrations were ∼1.2-fold higher in the FRG relative to the LAG and BR, it is likely that some micro-gradients in nutrient availability may have occurred in the FRG. For instance, small inputs from terrestrial runoff may have increased nutrient availability, which in turn stimulated heterotrophic BA. For example, Weinbauer et al. (2010a) reported increases of heterotrophic BA in response to small increases in nutrient availability due to terrestrial runoffs in another reef ecosystem. It is also likely that organic matter releases from benthic organisms may have caused micro-gradients in nutrient and carbon availability that stimulated ambient heterotrophic bacteria in the FRG. This is consistent with previous findings (Kline et al., 2006; Haas and Wild, 2010; Haas et al., 2011) that coral and macroalgal exudates of organic matter can enhance heterotrophic microbial abundance and activity, with implications for community structure (Haas et al., 2011, 2013; Nelson et al., 2013).

Consistent with BA, LNA and HNA containing cells followed decreasing trends along transects (**Figures 3E–L** and **4B,C**) and averaged 1.9 <sup>×</sup> 105 <sup>±</sup> 1.4 <sup>×</sup> <sup>10</sup><sup>5</sup> and 9.2 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>±</sup> 4.2 <sup>×</sup> <sup>10</sup><sup>4</sup> cells ml−1, respectively. HNA cells were more abundant than LNA cells (*t*-test on ranks, *p* < 0.01) and contributed 65 ± 9% of total BA along all transects. The proportion of HNA cells was higher in the FRG relative to the LAG and BR (KW with Dunn's test, *p* < 0.05), contributing up to 81% of total BA. In contrast, the proportion of LNA cells increased oceanward (KW, *p* >0.05) and contributed up to 44% of total BA in the FOR. Seymour et al. (2005)reported similar increased proportions of HNA cells in proximity to corals. HNA cells are reported to be large contributors to heterotrophic microbial activities, particularly in nutrient-replete conditions (Gasol et al., 1999; Lebaron et al., 2001; Servais et al., 2003). Therefore these results may indicate increased heterotrophic microbial activity in the FRG. However, recent surveys also have shown LNA cells, which are members of the abundant alphaproteobacterial clade SAR11, can be highly active, particularly in nutrient-depleted conditions (Zubkov et al., 2001; Jochem et al., 2004; Longnecker et al., 2005; Mary et al., 2008; Hill et al., 2010; Gomez-Pereira et al., 2013).

In the OB, BA averaged 5.6 <sup>×</sup> 105 <sup>±</sup> 0.9 <sup>×</sup> <sup>10</sup><sup>5</sup> cells ml−<sup>1</sup> and was significantly higher than sites along the transects (KW, *p* < 0.05; **Figure 5A**). HNA cells outnumbered LNA cells in the OB, representing 67 ± 8% of the total BA. This higher heterotrophic microbial abundance, and presumably activity, parallel the general pattern observed in the adjacent Cook's Bay (Nelson et al., 2011), that has relatively similar hydrological settings. In the OB, there were notable increases in suspended organic matter in the water column, as evidenced by reduced visibility (<5 m). Near-bottom currents that continuously re-suspend silty bottom of the bay, concomitant with small terrigenous inputs from a river near the tip of the bay, may explain such increases in suspended particles, as was reported previously in the OB and Cook's Bay (Wolanski and Delesalle, 1995; Hench et al., 2008; Nelson et al., 2011). Thus, potential nutrient supply in the form of suspended organic matter may have been important in sustaining high heterotrophic microbial abundance in the OB. The time-series in the OB revealed small temporal oscillations in BA with low date-to-date changes (range: 0.7- to 1.5-fold) and relatively stable proportions of HNA cells (KW, *p* > 0.05; **Figure 5A**), suggesting relative homogeneity in the abundance structure of the heterotrophic microbial communities.

## **AUTOTROPHIC MICROBES DISPLAY SPATIAL VARIABILITY AND SHORT-TERM TEMPORAL STABILITY**

For the transects PA averaged 3.6 <sup>×</sup> 104 <sup>±</sup> 3.0 <sup>×</sup> 104 cells ml−<sup>1</sup> and displayed consistent spatial trends (**Figures 6A–D**). Although PA gradually decreased ∼5-fold from the FOR toward the FRG (**Figure 7A**), no significant differences were detected among reef habitats (KW, *p* > 0.05). Similar to BA, there were dynamic spatial shifts in the phytoplankton community, with *Prochl* cells (3.4 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>±</sup> 2.9 <sup>×</sup> 104 cells ml−1) dominating total PA relative to *Syn* (2.6 <sup>×</sup> 103 <sup>±</sup> 2.6 <sup>×</sup> 103 cells ml−1), APP (1.1 <sup>×</sup> 103 <sup>±</sup> 0.9 <sup>×</sup> 103 cells ml−1) and nanoplankton (ANP; 142 <sup>±</sup> 108 cells ml−1) cells. All phytoplankton subsets except ANP cells followed decreasing trends from the FOR toward the CR and BR (**Figures 6E–T**), likely due to increased encounter rates with benthic filter feeders, as mentioned above. In the FRG, the dominant *Prochl* subset continuously decreased while the abundance of *Syn*, ANP, and APP cells markedly increased (**Figures 7B–E**). Organic matter and nutrient supply from benthic exudates may have stimulated microalgae with higher nutrient requirements in the FRG. However it is noteworthy that *Prochl*, which are known to cope better than *Syn* and eukaryotic phytoplankton in nutrient-depleted seawater (e.g., reviewed in Scanlan et al., 2009), were still prevailing in the FRG. This suggests that the nutrient supply was not sufficient to shift the phytoplankton community from *Prochl*-dominated toward *Syn*-dominated communities, as has been reported in other coastal tropical reef waters (e.g., reviewed in Charpy et al., 2012).

Similar to BA, higher overall mean PA was measured in the OB (2.2 <sup>×</sup> <sup>10</sup><sup>5</sup> <sup>±</sup> 1.1 <sup>×</sup> <sup>10</sup><sup>5</sup> cells ml−1; KW, *<sup>p</sup>* <sup>&</sup>lt; 0.05). In contrast to the transects, *Syn* was the most abundant phytoplankton, followed by *Prochl* (58 ± 4% and 39 ± 5%, respectively; **Figures 5B,C**). This shift in autotrophic microbial community structure and abundance supports the hypothesis that suspended organic matter and associated nutrients may have stimulated autotrophic cells with high nutrient demand. During the time-series in the OB, PA displayed similar temporal oscillations to BA with relatively low date-to-date changes (range: 0.4- to 2.5-fold; **Figure 5B**; KW, *p* > 0.05). Phytoplankton subsets also remained relatively unchanged (KW, *p* > 0.05), indicating a relatively stable autotrophic microbial community.

## **VIRUSES DISPLAY SPATIAL VARIABILITY AND SHORT-TERM TEMPORAL STABILITY DRIVEN BY BIOTIC AND ABIOTIC FACTORS**

Along cross-reef transects, VA averaged 1.2 <sup>×</sup> 107 <sup>±</sup> 0.5 <sup>×</sup> 107 viruses ml−<sup>1</sup> and was within ranges previously reported for coral reef waters (Paul et al., 1993; Seymour et al., 2005; Patten et al., 2006; Bouvy et al., 2012). Viral abundance followed similar spatial trends within the four transects (**Figures 8A–D**), with lowest values in the BR and LAG (**Figure 9A**). On average, VA was ∼2.2-

to 2.4-fold higher in the FRG and CR, relative to the BR and LAG (**Figure 9A**). DistLM for VA vs. biotic (i.e., BA and PA) and abiotic (i.e., DIN, SRP, longitude and latitude) variables indicated that these variables contributed to 24 and 38% of spatial variability in VA, respectively (**Table 2**). Among biotic variables, BA was the best predictor of VA, while SRP and DIN were the best abiotic predictors (**Table 2**). Interestingly, these results suggest that microbial host abundance only partially explained spatial VA distribution along transects, and that other unmeasured ecological processes may have influenced the distribution of VA. For example, small-scale changes in hydrological conditions may have influenced host distribution and metabolic activities, with implications for host-virus dynamics and subsequent spatial patterns of VAs. Alternatively, increases in VT time relative to host microbial

**Table 2 | Results of separate distance-based linear model (DistLM), with forward procedure, fitting viral abundance (VA), low nucleic acid viral V1 subset abundance, medium nucleic acid viral V2 subset abundance, high nucleic acid viral V3 subset abundance and viral production (VP) against biotic and abiotic variables.**


The analyses were conducted separately for biotic and abiotic variables as well as for samples collected along four cross-reef transects and during the time-series in the Oponuhu Bay. Pseudo-F and p-values were obtained by permutations (n = 999). P-values at a significance level of 0.05 are in *bold*. Proportional and cumulative percentages of variance are also reported. BA, bacterial abundance; LNA, low nucleic acid bacteria; HNA, high nucleic acid bacteria; SRP, soluble reactive phosphorus; DIN, dissolved inorganic nitrogen; PA, phytoplankton abundance; APP, autotrophic picoplankton; Syn, Synechococcus; Prochl, Prochlorococcus.

abundance may have dampened the relationship between viruses and host microbes.

Higher VA was measured in the OB relative to other reef sites, with an overall mean value of 1.9 <sup>×</sup> <sup>10</sup><sup>7</sup> <sup>±</sup> 0.9 <sup>×</sup> 107 viruses ml−1. During the time-series, VA displayed relatively similar temporal oscillations to BA and PA, with relatively low date-to-date changes (range: 0.3- to 2.4-fold; **Figure 10A**; KW, *p* > 0.05). DistLM for temporal VA vs. biotic (i.e., BA and PA) and abiotic (i.e., SRP, DIN and time) variables indicated that these variables explained 32 and 45% of temporal variability, with PA and SRP as their main predictors (**Table 2**). As outlined in the above section, the relatively weak relationships among VA, the measured biotic and abiotic variables suggest that other unmeasured variable(s) may be partially responsible for the observed temporal variability in virus–host dynamics, with implications for temporal distributions of VA.

## **SPATIOTEMPORAL DISTRIBUTION OF VIRAL SUBSETS INDICATES DYNAMIC VIRAL ASSEMBLAGES**

Based on their fluorescence properties, FC analysis revealed at least three viral subsets (i.e.,V1,V2, andV3), withV1 andV2 contributing to most of the VA (70 and 25%, respectively). This is consistent with subsets reported in other marine ecosystems (Baudoux et al., 2007; Evans et al., 2009; Brussaard et al., 2010; Mojica et al., 2014), however it should be noted that in some environments, only V1 and V2 are readily detected (Patten et al., 2006; Seymour et al., 2006; Payet and Suttle, 2008).

Along transects, V1 and V3 displayed relatively high spatial variability (**Figures 8E–H, M–P**), with significant increases in V1 in the FRG relative to the LAG and significant increases in V3 in the CR relative to the BR (KW with Dunn's test, *p* < 0.05; **Figures 9B,D**). V2 displayed relatively low spatial variability along transects (**Figures 8I–L**), with no significant differences detected among reef sites (KW, *p* > 0.05; **Figure 9C**).

Distance-based linear model for spatial V1, V2, and V3 abundances indicated significant relationships with biotic (i.e., HNA, LNA, *Syn*, *Prochl*, APP, and ANP) and abiotic (i.e., DIN, SRP and coordinates) variables (**Table 2**). Biotic and abiotic variables explained 39 and 47% of spatial variability in V1, respectively, with both LNA and *Prochl,* and both SRP and DIN as the best predictors (**Table 2**). Previous studies have reported that LNA cells contain a larger proportion of small bacteria from the alphaproteobacterial clade SAR11; these SAR11 bacteria typically co-occur with *Prochl* in nutrient-depleted waters (Hill et al.,2010; Gomez-Pereira et al., 2013). Therefore it may be that spatial patterns in V1 abundances were more associated with changes in these autotrophic and heterotrophic bacterial subsets across the reef. Biotic and abiotic variables explained 27 and 34% of spatial variability in V2, with HNA and SRP as the best predictors (**Table 2**). This suggests that viruses in the V2 subset were associated with heterotrophic

viral-mediated mortality of bacteria (VMB) and proportion of bacterial standing stock removed due to viral lysis (BSSr), and **(C)** amount of organic carbon released upon viral lysis (OCr) and viral turnover (VT) in the Oponuhu Bay. Error bars represent standard deviations.

microbes, with presumably high metabolic activity as outlined above. Biotic and abiotic variables explained 50 and 44% of the temporal variability in V3, respectively, with *Syn*, APP, and SRP as the best predictors (**Table 2**). This indicates that V3 is comprised of viruses that are associated with autotrophic microbial host cells.

Similar to temporal patterns in microbial community structure, proportions of V1, V2, and V3 remained relatively stable during the time-series in the OB, with only small date-to-date changes

(**Figure 10A**). This suggests that viral community structure was relatively homogeneous over time. DistLM for temporal V1, V2, and V3 abundances revealed significant relationships with biotic and abiotic variables in the OB (**Table 2**). Similar to temporal variations in overall VA, associations between viral subsets and biotic/abiotic variables tended to be stronger than those observed for transects. For V1, biotic and abiotic variables explained 21 and 39% of temporal variability, respectively, with *Syn* cells and SRP as the main predictors (**Table 2**). For V2, biotic and abiotic variables explained 35 and 51% of temporal variability, with LNA cells and SRP as the main predictors (**Table 2**). For V3, biotic and abiotic variables explained 50 and 44% of temporal variability in V3 abundances, respectively, with *Syn,* APP and SRP as the best predictors (**Table 2**). Different best predictors in the OB indicate that V1, V2, and V3 subsets may be influenced by different ecological factors than those within transects. However, V3 subset had similar predictors in both the OB and transects, suggesting these viruses are associated with changes in autotrophic microbial communities.

#### **SPATIOTEMPORAL DISTRIBUTION OF LYTIC ACTIVITY SUGGESTS VIRUSES IMPACT MICROBIAL MORTALITY AND CARBON CYCLING**

Along transects, estimates of lytic VP and VT averaged 7.3 <sup>×</sup> 109 <sup>±</sup> 4.2 <sup>×</sup> 109 viruses l−<sup>1</sup> <sup>d</sup>−<sup>1</sup> and 0.6 <sup>±</sup> 0.3 d−1, respectively (**Figures 11A,B**), and were within ranges previously reported for other marine ecosystems (Wilhelm et al., 2002; Poorvin et al., 2004; Winget et al., 2005; Weinbauer et al., 2009; Payet and Suttle, 2013). In general, VP and VT followed similar spatial trends along transects (**Figures 11A,B**), with highest and lowest values in the FRG and BR, respectively. Biotic (i.e., PA and BA) and abiotic (i.e., SRP, DIN and coordinates) variables explained 28 and 19% of spatial variability in VP, respectively, with LNA cells and SRP as the main predictors (**Table 2**). This implies phage infection of smaller bacteria may have been important, and is consistent with recent evidence showing phages are associated with highly abundant and small bacteria in the SAR11 and SAR116 clades in the oceans (Kang et al., 2013; Zhao et al., 2013).

Consistent with VP, VMB, BSSr and OCr displayed similar spatial trends across reef habitats (**Figures 10C–E**). On average, VMB removed 2.5 <sup>×</sup> 108 <sup>±</sup> 1.4 <sup>×</sup> 108 bacteria l−<sup>1</sup> <sup>d</sup>−<sup>1</sup> (range: 5.1 <sup>×</sup> 107– 5.2×10<sup>8</sup> bacteria l−<sup>1</sup> <sup>d</sup>−1), which accountedfor an estimated BSSr of 72 <sup>±</sup> 31% d−<sup>1</sup> (range: 24–133% d−1) and OCr of 4.9 <sup>±</sup> 2.7μg C l <sup>−</sup><sup>1</sup> d−<sup>1</sup> (range: 1.0–10.4 μgCl−<sup>1</sup> d−1) along cross-reef transects.

Assuming mean ambient DOC concentrations of 68μM as previously measured in this reef ecosystem (e.g., Nelson et al., 2011), viral lysis may contribute to ca. 2–15% of the pool of DOC in these reef waters. This implies that viral lytic activities may be important in fueling labile organic carbon and associated nutrient supply to other non-infected microbes in these oligotrophic reef waters. This is especially true for the BR, which is known to be a particularly carbon depleted habitat (Nelson et al., 2011). Therefore, slow release of organic matter upon viral lysis maybe providing essential carbon and nutrient supply for microorganisms in the BR and sustain the low levels of heterotrophic and autotrophic microbes observed in this study.

For the time-series in the OB, VP, and VT estimates mirrored temporal patterns in VA and averaged 4.1 <sup>×</sup> 1010 <sup>±</sup> 3.1 <sup>×</sup> 1010 viruses l−<sup>1</sup> <sup>d</sup>−<sup>1</sup> and 2.8 <sup>±</sup> 2.4 d−1, respectively (**Figures 11B,C**). VP and VT displayed small temporal oscillations, with no significant differences detected between sampling dates (KW, *p* > 0.05). Temporal changes in biotic (i.e., PA and BA) and abiotic (i.e., SRP, DIN and time) variables explained 82 and 15% of variability VP, respectively, with APP, *Syn,* and SRP as the main predictors (**Table 2**). These results indicate that viral lytic production was strongly associated with changes in primary producers in the OB.

Consistent with VA and VP, higher estimates for VMB, BSSr, and OCr were detected in the OB compared to transects (**Figures 11B,C**). On average, VMB was responsible for the removal of 1.4 <sup>×</sup> 109 <sup>±</sup> 1.0 <sup>×</sup> <sup>10</sup><sup>9</sup> bacteria l−<sup>1</sup> <sup>d</sup>−1. This corresponds to an estimated BSSr of 152 <sup>±</sup> 108% d−<sup>1</sup> (range: 42–367% d−1) and OCr of 27.8 <sup>±</sup> 19.6 <sup>μ</sup>gCl−<sup>1</sup> <sup>d</sup>−<sup>1</sup> (range: 6.9–61.2 μgCl−<sup>1</sup> d−1) in the OB. VMB, BSSr, and OCr followed similar temporal trends throughout the time-series, with relatively low date-to-date changes (**Figures 11B,C**). Similar to patterns in BA, PA and VA, these results suggest that viral lytic activity was relatively stable over time in the OB.

Again, assuming a mean DOC concentration of 68 μM (e.g., Nelson et al., 2011), these cellular lysis products were responsible for between 10 and 90% of ambient DOC levels in the OB. Thus, viral infection of heterotrophic bacteria may be an important source of DOC and associated nutrients for other non-infected microbes in the OB.

### **METHODOLOGICAL CONSIDERATIONS**

These results should be interpreted in the context of several limitations. We used FC to identify viral subsets, according to an established protocol that has been used in other studies (Brussaard et al., 2010). Although significant trends among viral subsets and biotic communities were detected, their identities still remain unclear. Recently, Martínez-Martínez et al. (2014) genetically characterized three viral subsets with relatively similar FC signatures to those reported in this study. While increased proportions of viruses infecting eukaryotic phytoplankton were detected from the V1 and V3 subset (Martínez-Martínez et al., 2014), these subsets still contained significant proportions of phages, highlighting the poor resolution of FC in distinguishing particular viruses associated with certain host cells. Thus, it is possible that this low resolution may have masked potential correlations among viruses and their hosts, explaining relatively weak relationships among viruses, biotic and abiotic variables. However, FC has become a standard methodology for studying of virus–host interactions, but can be extended and complemented by other microscopic techniques and molecular approaches.

During the viral-reduction approach, filtration steps required for reducing virus–host contact rates may have altered nutrient availability and microbial processes, potentially influencing results. In addition, estimates of BS and cellular carbon quota used to infer viral-induced mortality and carbon cycling are likely to fluctuate across gradients of microbial productivity. Despite these caveats, the viral-reduction approach has been successfully applied in various aquatic environments and has been shown to be a robust and straightforward approach for estimating viral lytic activity (Weinbauer et al., 2010b and references therein). The repeatable temporal patterns measured during our 3 week time-series study suggest that this viral-reduction approach can be applied to investigate lytic viral activity in response to variability in host abundance and environmental conditions.

Further work is needed to improve detection of viruses infecting microbial hosts in natural assemblages. In particular, the development of high-throughput methods to routinely detect virally infected microbial hosts from environmental samples, in conjunction with developments of specific molecular probes to target potential host–virus systems will provide novel insights on viral dynamics and their impacts in the oceans. Recently, the application of a new culture-dependent and independent approaches has allowed direct detection of viruses infecting microbial host isolates and offered new exciting perspectives for enabling simultaneous detection of host–virus interactions (Deng et al., 2012; Allers et al., 2013). However, the applicability to study a broader range of natural samples has still to be shown.

## **SUMMARY**

This study is the first to report the abundance, distribution and ecological impact of viruses in the coral reef waters of Moorea Island. Our data revealed distinct short-term spatiotemporal changes in VA and activity and demonstrated that these changes were linked to microbial host abundances and environmental variables. This work also confirmed general findings from other studies which have suggested that small shifts in host abundance and activity may be important in driving VA and lytic activity in marine systems (Fuhrman, 1999; Weinbauer, 2004; Seymour et al., 2005; Suttle, 2007; Brussaard et al., 2008; Clasen et al., 2008; Payet and Suttle, 2008; Rowe et al., 2008; Evans et al., 2009; Winget et al., 2011; Winter et al., 2012; Payet and Suttle, 2013).

Analysis of short-term temporal patterns in VA and lytic production in OB indicated persistent VA and infection. These findings confirm recent time-series studies that have also observed steady-state temporal dynamics of lytic viral activity (Winget and Wommack, 2009; Winget et al., 2011).

Highest VA and lytic activity as well as highest microbial host abundances were reported in FGRs as well as in OB, likely due to microgradients in nutrient availability.

Viral lysis was estimated to kill a significant fraction of heterotrophic microbes (%BSSr: 24–367%) daily. These mortality estimates are substantially higher than those estimated by Bouvy et al. (2012) in another reef ecosystem in French Polynesia, but were in agreement with other studies in other marine ecosystems (Wilhelm et al., 2002; Winget et al., 2005, 2011; Evans et al., 2009; Winget and Wommack, 2009; Evans and Brussaard, 2012; Payet and Suttle, 2013). Given that Bouvy et al. (2012) used frequency of visibly infected cells to infer mortality estimates through transmission electronic microscopy, it may be that this approach underestimated lytic viral impacts, as it heavily relies on specific conversionfactors and potentially lacks of resolution due to sample preparations (e.g., see Weinbauer et al., 2002).

Notably, our data demonstrate that viral lysis substantially contributes to the overall pool of DOC (OCr: 1.0–62 μgCl−<sup>1</sup> d−1) available to other microbes in these oligotrophic coral reef waters. Our estimates of OCr due to viral lysis were ca. 1- to 90-fold higher than previous reports in oligotrophic polar waters (Evans and Brussaard, 2012; Payet and Suttle, 2013), but within the range of other studies in marine environments (Wilhelm et al., 2002; Winget et al., 2005, 2011).

In conclusion, this study demonstrates that viruses have a key role in both top down and bottom up control of microbial communities in coral reef seawater.

#### **ACKNOWLEDGMENTS**

We would like to thank the staff of USR 3278 CRIOBE CNRS-EPHE for the use of their facilities, logistical support and valuable input in Moorea. We would also like to thank the reviewers for their thoughtful suggestions on the manuscript. This research was supported by the U.S. National Science Foundation (NSF) grants OCE-0960937 and OCE-1130786 (Rebecca L. Vega Thurber and Deron E. Burkepile) and by the Institute for Pacific Coral Reefs (IRCP) – Tahiti Perles research grant (Jérôme P. Payet).

## **REFERENCES**


transport and retention processes on Moorea, French Polynesia. *Oceanography* 26, 52–63. doi: 10.5670/oceanog.2013.45


Zar, J. H. (1999). *Biostatistical Analysis*. Upper Saddle River, NJ: Prentice-Hall.

Zhao, Y., Temperton, B., Thrash, J. C., Schwalbach, M. S., Vergin, K. L., Landry, Z. C., et al. (2013). Abundant SAR11 viruses in the ocean. *Nature* 494, 357–360. doi: 10.1038/nature11921

Zubkov, M. V. B., Fuchs, B. M., Burkill, P. H., and Amann, R. (2001). Comparison of cellular and biomass specific activities of dominant bacterioplankton groups in stratified waters of the Celtic Sea. *Appl. Environ. Microbiol.* 67, 5210–5218. doi: 10.1128/AEM.67.11.5210-5218. 2001

**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.

*Received: 14 May 2014; accepted: 01 September 2014; published online: 23 September 2014.*

*Citation: Payet JP, McMinds R, Burkepile DE and Vega Thurber RL (2014) Unprecedented evidence for high viral abundance and lytic activity in coral reef waters of the South Pacific Ocean. Front. Microbiol. 5:493. doi: 10.3389/fmicb.2014.00493*

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

*Copyright © 2014 Payet, McMinds, Burkepile and Vega Thurber. This is an openaccess 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.*

## Filamentous phages of Ralstonia solanacearum: double-edged swords for pathogenic bacteria

## *TakashiYamada\**

Department of Molecular Biotechnology, Graduate School of Advanced Sciences of Matter, Hiroshima University, Higashi-Hiroshima, Japan

#### *Edited by:*

Heather K. Allen, National Animal Disease Center, USA

#### *Reviewed by:*

Loren John Hauser, The University of Tennessee, Oak Ridge National Laboratory, USA Jasna Rakonjac, Massey University, New Zealand

#### *\*Correspondence:*

Takashi Yamada, Department of Molecular Biotechnology, Graduate School of Advanced Sciences of Matter, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima 739-8530, Japan e-mail: tayamad@hiroshima-u.ac.jp

Some phages from genus Inovirus use host or bacteriophage-encoded site-specific integrases or recombinases establish a prophage state. During integration or excision, a superinfective form can be produced. The three states (free, prophage, and superinfective) of such phages exert different effects on host bacterial phenotypes. In Ralstonia solanacearum, the causative agent of bacterial wilt disease of crops, the bacterial virulence can be positively or negatively affected by filamentous phages, depending on their state. The presence or absence of a repressor gene in the phage genome may be responsible for the host phenotypic differences (virulent or avirulent) caused by phage infection. This strategy of virulence control may be widespread among filamentous phages that infect pathogenic bacteria of plants.

**Keywords: filamentous phage, integration, phytopathogen,** *Ralstonia solanacearum,* **virulence change**

"fmicb-04-00325" — 2013/10/31 — 16:31 — page 1 — #1

## **FILAMENTOUS PHAGES AND PATHOGENIC BACTERIA**

Bacteriophages belonging to the genus *Inovirus* are filamentous particles containing a circular single-stranded (ss) DNA genome. Infection with this kind of phage does not cause host cell lysis, but establishes a persistent association between the host and phage, producing and releasing phage particles from the growing and dividing host cells. In general, the genome of inoviruses, represented by *Escherichia coli* F-pillus specific phage Ff (f1, fd or M13), is organized in a modular structure, in which functionally related genes are grouped together (Horiuchi et al., 2009; Rakonjac et al., 2011). Three functional modules are always present: the replication module, the structural module, and the assembly and secretion module. The replication module contains the genes encoding rolling-circle DNA replication and singlestrand DNA (ssDNA) binding proteins *gII*, *gV*, and *gX* (Horiuchi et al., 2009). The structural module contains genes for the major (*gVIII*) and minor coat proteins (*gIII*, *gVI*, *gVII*, and *gIX*), and gene *gIII* encodes the host recognition or adsorption protein pIII (Wang et al., 2006). The assembly and secretion module contains the genes for morphogenesis and extrusion of the phage particles (*gI* and *gIV*; Marvin, 1998). Gene *gIV* encodes protein pIV, an aqueous channel (secretin) in the outer membrane, through which phage particles exit from the host cells (Marciano et al., 1999). Some phages encode their own secretins, whereas others use host products (Davis et al., 2000).

Because inoviruses coexist with their host cells, infection by these phages can mediate conversion of the host bacterial phenotypes in various ways. In pathogenic bacteria of either animals or plants, virulence is frequently affected by phage infection. For example, infection of *Xanthomonas campestris* pv. *oryzae* NP5850 by the filamentous phages Xf and Xf2 enhanced virulence, possibly because of overproduction of extracellular polysaccharides (EPS) by the phage-infected bacterial cells (Kamiunten and Wakimoto, 1982). Tseng et al. (1990) also reported that infection of *X*. *campestris* pv. *campestris* by the filamentous phage Lf increased virulence by promoting EPS production. Filamentous phages are assembled at the host cell surface and secreted into the environment. However, once then cells form colonies on the semi-solid medium (and possibly within the liquid medium), some fractions of secreted phage population are bound to stay trapped in the colony, potentially accumulating to high concentrations and forming a matrix surrounding the cells in the colony. These trapped phage particles may serve to cross-link cells to give high densities and induce biofilms. This situation was reportedfor small colony variantformation in *Pseudomonas aeruginosa* depending on phage Pf4 activity (Webb et al.,2004; Rice et al., 2009). More direct involvement of filamentous phages in host virulence is well characterized in *Vibrio cholerae.* The pathogenicity of this severe diarrheal disease-causing bacterium depends on two key virulence factors, the toxin co-regulated pilus (TCP) and cholera toxin. Cholera toxin genes are encoded on the filamentous phage CTXφ and introduced into bacterial cells by phage integration mediated by the host *dif*/XerCD recombinase system (Huber and Waldor, 2002; Davis and Waldor, 2003). In *Ralstonia solanacearum*, infection by φRSS1 induced the early expression of *phcA*, a global virulence regulator, and also enhanced twitching motility (Addy et al., 2012b).

Contrasting with these virulence-enhancing effects of φRSS1, loss of virulence was also reported in *R. solanacearum*. *R. solanacearum* completely lost virulence through infection with two other filamentous phages φRSM1 and φRSM3 (Addy et al., 2012a). Many virulence factors were significantly reduced in φRSM-infected cells. These opposing effects of different filamentous phages on *R. solanacearum* virulence makes it an ideal study model system for understanding the effect of filamentous phage on their hosts. Here I will describe the role of filamentous phage in the virulence of *R. solanacearum* and suggest a causative relationship between a phage-encoded transcriptional repressor and *R. solanacearum* pathogenicity.

## *Ralstonia solanacearum* **AND BACTERIAL WILT**

*Ralstonia solanacearum* is a Gram-negative β-proteobacterium that causes bacterial wilt disease in many important crops including tomato, potato, tobacco, and eggplant. Because of its wide geographic distribution and unusually broad host range (more than 50 plant families), it is responsible for significant crop losses worldwide (Hayward, 2000). Once the bacteria enter a susceptible host, they colonize the intercellular spaces of the root cortex and vascular parenchyma. The bacteria eventually enter the xylem and spread into the upper parts of the plant, causing wilt (Vasse et al., 2000; Kang et al., 2002; Yao and Allen, 2007). The development of bacterial wilt disease depends on bacterial pathogenicity and virulence (Carney and Denny, 1990; Denny, 2006). *R. solanacearum* virulence is additive, complex, and involves the production of multiple virulencefactors (Schell,2000; Genin and Boucher,2002). For example, exopolysaccharide I (EPSI), a large nitrogen-rich acidic exopolysaccharide (Lavie et al., 2002), is thought to be an important virulence factor. It enhances the speed and extent of stem infection spreading from the root (Saile et al., 1997) and is presumed to cause wilting by restricting water flow through xylem vessels (Garg et al., 2000). In addition to EPSI, *R. solanacearum* secretes enzymes that degrade the plant cell wall through the type II secretion system (T2SS). Pectinolytic enzymes fragment pectin into oligomers, which act as a substrate for bacterial growth (Tans-Kersten et al., 2001). The breakdown of pectin enhances virulence by facilitating bacterial movement through pectin-rich regions such as vascular bundles (Gonzalez and Allen, 2003). Cellulolytic enzymes also facilitate bacterial invasion of roots and/or penetration of xylem vessels by degrading cellulosic glucans in the cell wall (Liu et al., 2005). In addition to T2SS-mediated secreted proteins, the type IV pilus (Tfp) is believed to be another virulence factor of *R. solanacearum* (Davis and Waldor, 2003). This protein forms a surface appendage that is responsible for twitching motility and polar attachment to host cells or to plant roots, and enhances the severity of wilt disease (Liu et al., 2001; Kang et al., 2002).

Expression of the pathogenesis and virulence genes in *R. solanacearum* is controlled by a complex regulatory network (Schell, 2000; Genin and Boucher, 2002; Denny, 2006) and is drastically affected by various environmental factors. The regulation is outlined as follows: the transcriptional regulator PhcA plays a critical role in the regulatory network. Abundant PhcA activates production of multiple virulence factors such as EPSI and cell wall degrading enzymes (CWDE). PhcA is activated by a quorum sensing system mediated by the two-component regulatory system PhcS/PhcR that responds to thereshold levels of 3-OH palmitic acid methylester (3-OH PAME), an autoinducer of quorum sensing that controls virulence. Therefore, the levels of 3-OH PAME, cell density, as well as cell surface nature all affect virulence in *R. solanacearum*.

## **THREE STATES OF FILAMENTOUS PHAGE φRSS WITH DIFFERENT EFFECTS ON HOST VIRULENCE**

φRSS1 was isolated from a soil sample collected from tomato crop fields (Yamada et al., 2007). φRSS1 particles have a flexible filamentous shape 1,100 ± 100 nm in length and 10 ± 0.5 nm in width, giving a morphology resembling coliphage Ff (M13, f1 or fd; Buchen-Osmond, 2003; ICTVdB). The φRSS1 particles contain a ssDNA genome (6,662 nt; DDBJ accession no. AB259124), with a GC content of 62.6%. There are 11 open reading frames (ORFs), located on the same strand (**Figure 1A**). The φRSS1 gene arrangement is consistent with the general arrangement of Ff phages. Genomic Southern blot hybridization showed several examples of φRSS1-related sequences integrated in the genomes of various *R. solanacearum* strains (Yamada et al., 2007). A φRSS1-related phage (designated φRSS0) was induced and isolated from one such crosshybridizing strain (C319) by infection with another phage (jumbo phage φRSL1). The DNA sequence of φRSS0 was very similar to φRSS1, but contained an extra 626 nt at φRSS1 position 6,628, next to the intergenic region (IG), giving an entire genomic size of 7,288 nt (GenBank accession no. JQ408219). Within the φRSS0 extra region, an ORF (ORF13) of 468 nt, corresponding to 156 amino acid residues, in a reversed orientation compared with the other ORFs, was found (**Figure 1A**). The amino acid sequence of ORF13 showed similarity to DNA-binding phage transcriptional regulators (accession no. B5SCX5, E-value = 1e-29).

Using inverse PCR with the new phage nucleotide sequences as primers, the prophage (φRSS0)-junctions (*att*L and *att*R) in strain C319 were obtained and their nucleotide sequences determined. It was found that both *att*L and *att*R contained repeated elements, corresponding to the *dif* sequence of *R. solanacearum* GMI1000 (Carnoy and Roten, 2009). This repeated sequence, 5- - TATTT AACAT AAGAT AAAT-3- , was also found at the 3 end of ORF13 on the RSS0 genome, suggesting that it serves as *att*P.

Taken together, these results indicated that φRSS1 (with a genome size of 6,662 nt) is a truncated form of the larger phage φRSS0 (with a genome size of 7,288 nt). The 626 nt φRSS0 sequence missing from φRSS1 contains *att*P (corresponding to the *dif* sequence) and ORF13, a possible regulatory gene (Yamada, 2011; Tasaka et al., unpublished). φRSS0 is integrated at the *dif* site, similar to CTXφ of *V. cholerae*, which uses the host XerCD recombination system (Huber andWaldor, 2002). ORF13 encoded on φRSS0 may function as a phage immunity factor, because strain C319 (φRSS0 lysogen) is resistant to secondary infection by φRSS0. C319 is susceptible to φRSS1, thus φRSS1 (without ORF13) may be an escaped superinfective phage. These three states of φRSS phages and their interrelationships are shown in **Figure 1B**.

Upon infection by the φRSS1 phage, the host *R. solanacearum* cells showed several abnormal behaviors, including less turbidity and frequent aggregation in the liquid culture, less coloration of colonies on plates, and a decreased growth rate (approximately 60% of the normal rate). More interestingly, φRSS1-infected cells showed enhanced virulence on tobacco (Yamada et al., 2007) and tomato plants (Addy et al., 2012b). In the case of strain C319 (φRSS0 lysogenic), inoculated tobacco plants showed wilting symptoms of grade 2–3 at 14 days post-inoculation (p.i.), whereas tobacco plants inoculated with φRSS1-infected C319 cells

"fmicb-04-00325" — 2013/10/31 — 16:31 — page 2 — #2

"fmicb-04-00325" — 2013/10/31 — 16:31 — page 3 — #3

arrows oriented in the direction of transcription. The functional modules for replication (R), structure (S), and assembly secretion (A-S) are indicated according to the M13 model (Marvin, 1998). The region containing the attP sequence is also indicated. **(B)** Interrelationship between three states of φRSS phages. The phage genomic DNA is shown in a circular form where most genes are not shown. φRSS0 is equipped with a 626-nt element containing ORF13, within which attP (dif) is located. This element is

borders are indicated as attL and attR, respectively. This integration (reversible) is mediated by the host XerCD system. φRSS1 may be produced directly from RSS0φ. The three states of φRSS0-type phage (φRSS0, φRSS1, and φRSS0 prophage) affect host R. solanacearum cells differently after infection, especially in host virulence. Compared with wild-type virulence (+), φRSS1 enhances (++) and φRSS0 reduces (-/±) the host virulence.

wilted earlier; grade 2–3 symptoms were observed at 10 days p.i. and plants were almost dead at 14 days p.i. (Yamada et al., 2007). Effects on host virulence by infection with φRSS0 in its free form (not prophage) were also examined. To make wilting symptoms clear, tomato-tropic *R. solanacearum* strain (MAFF 106603) in tomato experimental system was used. The cells were infected with either φRSS0 (free) or φRSS1. The physiological features of φRSS0-infected *R. solanacearum* MAFF 106603 cells were almost the same as φRSS1-infected MAFF 106603 cells, except that the φRSS0-infected cells formed colonies of more mucoid appearance on CPG plates. When MAFF 106603 (wild-type) cells were inoculated into the major stem of tomato plants, all 12 plants showed wilting symptoms as early as 3 days p.i. and died 5–7 days p.i. φRSS1-infected cells of MAFF 106603 inoculated into tomato in the same way caused wilting earlier, at 2 days p.i., and all 12 plants died by 5 days p.i. In contrast, tomato plants inoculated with φRSS0-infected cells showed wilting symptoms much later: most plants (10 of 12) survived after 7 days and a few plants did not show any symptoms until 23 days p.i. Therefore, φRSS0 infection caused reduced virulence in host bacterial cells (Tasaka et al., unpublished). The virulence-enhancing effects by φRSS1 infection can be explained as follows: surface-associated φRSS1 particles (or phage proteins) may change the surface nature (hydrophobicity) of host cells to generate a high local cell density, resulting in early activation of *phcA*, the global virulence regulator, or lack of *orf13*, which is absent from the φRSS1 genome (Addy et al., 2012b). The reduced virulence observedforφRSS0-infected cells may be caused by the function(s) of ORF13 encoded by φRSS0. These results are summarized in **Table 1**.

## **EFFECT OF FILAMENTOUS PHAGE φRSM ON VIRULENCE OF** *Ralstonia solanacearum*

φRSM1 is also a soil-isolated filamentous phage 1,400 ± 300 μm long and 10 ± 0.7 nm wide (Yamada et al., 2007). The infection cycle of φRSM1 phage resembles that of φRSS1. The genome of φRSM1 is 9,004 nt long (DDBJ accession No. AB259123) with

#### **Table 1 |Three states of filametous phages and their effects on host virulence.**


a GC content of 59.9%. There are 12 putative ORFs located on the same strand and three on the opposite strand. The φRSM1 genes are shown in **Figure 2A**, in comparison with the conserved gene arrangement of M13-like phages (Kawasaki et al., 2007). Here, ORF13, ORF14, and ORF15 (reversely oriented) are inserted between ORF12, corresponding to pII as a replication protein, and ORF1, corresponding to a ssDNA-binding protein like pV, in the putative replication module. ORF13, ORF14, and ORF15 show amino acid sequence similarity to a proline-rich transmembrane protein, a resolvase/DNA invertaselike recombinase, and a putative phage repressor, respectively (Kawasaki et al., 2007;Addy et al., 2012a). There are two additional ORFs (ORF2 and ORF3) between the replication and structural modules. The functions of these ORF-encoded proteins are not known. In genomic Southern blot hybridization, two different types of φRSM1-related prophage sequences were detected in *R. solanacearum* strains. Strains of type A include MAFF211270 and produce φRSM1 itself, and strains of type B (giving different restriction patterns) are resistant to φRSM1 infection, but are susceptible to φRSM3 (see below). By determining the nucleotide sequences of junction regions of the φRSM1-prophage in the MAFF 211270 chromosomal DNA, an *att*P/*att*B core sequence was identified as 5- -TGGCGGAGAGGGT-3- , corresponding to positions 8,544–8,556 of φRSM1 DNA, located between ORF14 and ORF15. Its nucleotide sequence is identical to the 3- -end of the host *R. solanacearum* gene for serine tRNA(UCG) in the reverse orientation. A φRSM1-like prophage (type B) in strain MAFF 730139, designated φRSM3, was obtained by PCR amplification using appropriate primers containing these *att* sequences (Askora et al., 2009). Compared with the φRSM1 genome, the φRSM3 prophage sequence (8,929 nt) is 75 nt shorter. The sequences show 93% nucleotide identity and major differences are found within two regions; positions 400–600 and positions 2,500–3,000 in the φRSM1 sequence. The former region corresponds to ORF2, which is inserted between the replication module (R) and the structural module (S), and has no similarity between the two phages. The latter falls into the possible D2 domain of ORF9 (pIII), which determines the host range. All other ORFs identified along the φRSM3 are highly conserved between two phages (over 90% amino acid identity). It is interesting that the amino acid sequence of ORF14 (putative DNA invertase/recombinase) is 100% identical in the two phages. The gene arrangement of φRSM3, which is almost the same as φRSM1, is shown in **Figure 2A**.

As described above, the genomes of φRSM phages are sometimes integrated in the host genome. Askora et al. (2011) demonstrated that the integration is mediated by the phageencoded recombinase (ORF14 of φRSM1/φRSM3), which has significant homology to resolvases/DNA invertases (small serine recombinases), with *att*P/*att*B corresponding to the 3 end of the host serine tRNA(UCG) gene in the reverse orientation. This is the first case of filamentous phages demonstrated to integrate into the host genome by its endogenously encoded integrase (Askora et al., 2012). The same unit of integration (φRSM Int/*att*P) was found in a *Ralstonia pickettii* 12J phage and in *Burkholderia pseudomallei* 668 prophages (Askora et al., 2012). Together with these phages, it would not be surprising if similar Int-containing filamentous phages occur widely in nature.

Infection by φRSM1 or φRSM3 establishes a persistent association between the host and the phage. Upon infection by φRSM phages, the host cells showed some abnormal behaviors and characteristics, such as frequent aggregation, dark coloration, and relatively small colony size, as observed in φRSS infection. When cells of MAFF 106611 (φRSM3 lysogenic strain) or MAFF 106603 (not lysogenic) were inoculated into tomato plants, all 20 plants showed wilting symptoms as early as 3 days p.i., whereas none of the 20 tomato plants inoculated with free-φRSM3-infected cells (for example, MAFF 106603) showed any wilting symptoms until 4 weeks p.i. (Addy et al., 2012a). This loss of virulence effect of φRSM3 infection can be explained in three ways: (i) reduced twitching motility and reduced amounts of type IV pili (Tfp), (ii) lower levels of β-1,4-endoglucanase (Egl) activity and EPS production, and (iii) reduced expression of certain virulence/pathogenicity genes (*egl*, *pehC*, *phcA*, *phcB*, *pilT,* and *hrpB*) in the infected cells (Addy et al., 2012a). This is supported restoring virulence in φRSM3 lysogen by deletion of φRSM3 encoded *orf15*, the gene for a putative repressor-like protein, was disrupted (Addy et al., 2012a). Therefore, ORF15 of φRSM3 may repress host genes involved in pathogenicity/virulence and consequently result in loss of virulence. With different strains as hosts, φRSM1 also gave similar results. The φRSM states and interaction with the host genome can be depicted similarly to φRSS phages, as shown in **Figure 2B**. These results are summarized and compared with the three states of φRSS in **Table 1**.

## **PERSPECTIVE AND HYPOTHESIS**

"fmicb-04-00325" — 2013/10/31 — 16:31 — page 4 — #4

As seen here, for *R. solanacearum*, filamentous phages such as φRSS and φRSM are double-edged swords; sometimes they help bacteria to infect plants by enhancing bacterial virulence, and sometimes they interrupt bacterial infection of plants by repressing host genes involved in virulence. The contradictory effects of these phages may largely depend on the presence or absence of a phage-encoded regulatory protein. Two questions arise here: (i) How does the regulatory affect on the host genes; working alone, with other phage factors, or with host factors? (ii) How does such a regulatory gene become acquired by or lost from the phage genome? Concerning the first question, as shown in **Figure 1B**, attP is located within ORF13 on φRSS0 DNA, and after integration at attB on the host genome, a truncation of ORF13 (at the C-terminus) occurs. By creating a new stop codon in the reading frame, the size of ORF13 reduced from 156 to 130 aa with a 26-aa truncation at the C-terminus (Tasaka et al., manuscript in preparation). A DNA-biding motif (Helix-Turn-Helix) is located in the N-terminal moiety and the C-terminal region may have some regulatory function (such as ligand-binding). This suggests a functional difference of the ORF13 protein before and after integration. One possibility is that the full length ORF13 (ORF15 in φRSM phages) expressed from free phages may function to preferentially regulate host genes and the truncated (or modified) form expressed from the prophage may function to stabilize the prophage state and phage immunity, protecting against infections by related phages (Hypothesis 1). This hypothesis is compatible

"fmicb-04-00325" — 2013/10/31 — 16:31 — page 5 — #5

with the observation that once a φRSS and φRSM prophage state was established, the phage genomic DNA and phage particles seldom appeared in the lysogenic strains. Like φRSM, the DNA or the phage particles are not identified in the lysogen, even though the orf15 encoding the putative repressor ORF15 is not changed before and after the host integration. Because ORF14 integrase (serine recombinase) of φRSM phages likely mediates both integrative and excessive recombinations (Askora et al., 2011), some additional factors are required to mediate prophage replication or excision. The function, regulatory mechanism, and effect on virulence of φRSS orf13 or φRSM orf15 remain to be investigated by direct expression of the corresponding gene in an appropriate host strain. In our preliminary trial where the coding region of ORF15 of φRSM3 (ORF13 of φRSS0) was expressed from a plasmid under the control of lacP and introduced into appropriate host strains, no transformants with a correct construct appeared (colonies that appeared on the selection plates after transformation all contained deleted inserts). One of the explanation for this is putative toxic effect of ORF13 and ORF15 on the host when expressed under these conditions. Some additional factors encoded on the phage genome may be involved in the appropriate regulation, interacting with ORF13 or ORF15 (Hypothesis 2). Further studies with mutated constructs of ORF13 or ORF15 are required to test these hypotheses.

As for question of loss of a repressor protein, a 626-nt sequence unit containing *orf13* and *att*P detected inφRSS0 and missingfrom φRSS1 plays a crucial role in φRSS dynamics. The origin of such a sequence and the mechanism how it comes in or out of the phage are largely unknown. However, the possibility of two forms from a phage is important. Apparently, φRSS1-infected bacterial cells have an advantage in the pathogenic lifestyle. Nevertheless, the virulence is not always necessary for this soil-borne bacterium. Infection of φRSS0 provides the host cells with a sophisticated mechanism to control their virulence. Similar mechanisms may function in other pathogenic bacteria (Hypothesis 3). To test this hypothesis, various systems involving pathogenic bacteria and their filamentous phages should be examined. For example, φRSS1-like superinfective phage Cf1tv spontaneously appeared from the Cf1t lysogenic strain of *Xanthomonas campestris* pv. *citri* (Kuo et al.,1994). Unfortunately, nucleotide sequence information is not available for this phage. Similar kinds of phage involvement in host virulence regulation may be universal, because φRSS- or φRSM-related sequences are frequently found in various bacterial genomic sequences in the databases, including *R. pickettii* (CP001645), *Ralstonia syzygii* (FR854090), *Burkholderia rhizoxinica* (FR687359), *Pectobacterium wasabiae* (CP001790), and *Erwinia carotovora* (BX950851). There are also other filamentous phages that have lysogenic cycles, including *X. campestris* phages Cf1c (Kuo et al.,1991), Cf1t (Kuo et al.,1987a,b), Cf16v1 (Dai et al., 1980), and φLf (Lin et al., 2000); *Xylella fastidiosa* phage Xfφ-f1 (Simpson et al., 2000); *Yersinia pestis* phages CUSφ-2 (Gonzalez and Allen, 2003) and Ypfφ (Derbise et al., 2007); Nf of *Neisseria meningitidis* (Kawai et al., 2005), and *V. cholerae* phages VGJφ (Campos et al., 2003) and VCYφ (Xue et al., 2012). The host bacteria of these phages are plant or animal pathogens.

## **ACKNOWLEDGMENTS**

This study was supported in part by Research and Development Projects for Application in Promoting New Policy of Agriculture, Forestry, and Fisheries (No. 250037B), and by JST/BIOTEC Strategic Research Cooperative Program on Biotechnology.

## **REFERENCES**


of the exopolysaccharide I biosynthetic operon of *Ralstonia solanacearum*. *J. Bacteriol.* 182, 6659–6666. doi: 10.1128/JB.182.23.6659-6666.2000


"fmicb-04-00325" — 2013/10/31 — 16:31 — page 6 — #6


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

*Received: 15 August 2013; accepted: 14 October 2013; published online: 04 November 2013.*

*Citation: Yamada T (2013) Filamentous phages of Ralstonia solanacearum: doubleedged swords for pathogenic bacteria. Front. Microbiol. 4:325. doi: 10.3389/fmicb. 2013.00325*

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

*Copyright © 2013 Yamada. 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.*

"fmicb-04-00325" — 2013/10/31 — 16:31 — page 7 — #7

## Microarray tools to unveil viral-microbe interactions in nature

#### *Fernando Santos 1, Manuel Martínez-García1, Víctor Parro2 and Josefa Antón1 \**

*<sup>1</sup> Departamento de Fisiología, Genética y Microbiología, Universidad de Alicante, Alicante, Spain*

*<sup>2</sup> Departamento de Evolución Molecular, Centro de Astrobiología (INTA-CSIC), Madrid, Spain*

#### *Edited by:*

*Heather K. Allen, National Animal Disease Center, USA*

#### *Reviewed by:*

*Loren John Hauser, University of Tennessee Knoxville, USA Carl James Yeoman, Montana State University, USA*

#### *\*Correspondence:*

*Josefa Antón, Departamento de Fisiología, Genética y Microbiología, Universidad de Alicante, Apartado 99, 03080 Alicante, Spain e-mail: anton@ua.es*

The interactions between viruses and their microbial hosts play a central role in the control of microbial communities in nature. However, the study of such interactions within the uncultured majority is technically very challenging. Here, we review how microarray tools can be used to analyze the interactions between viruses and their microbial hosts in nature, away from laboratory pure culture-based models. We show examples of how DNA arrays have been used to study the expression of viral assemblages in natural samples, and to assign viruses to hosts within uncultured communities. Finally, we briefly discuss the possibilities of protein and glycan arrays to gain insight into the ways microbes interact with their viruses.

**Keywords: phage-bacteria interactions, microarray analysis, virus ec, microbiology, virus-, virus-host interaction**

## **INTRODUCTION**

As stated in the editorial of this issue of *Frontiers in Microbiology*, the nature of life is change (Allen and Abedon, 2013). Changes occur at every level of biological complexity and are especially fast and dramatic within the microbial world. These changes affect the complex networks of interactions between viruses and their hosts that structure microbial communities in nature. Here, we review how microarrays can be used to analyze the interactions between uncultured viruses and their microbial hosts in nature, away from laboratory pure culture-based models. In addition to providing some examples on the successful use of this technology, we offer some ideas about its potential applications.

Microarray technology was developed in the mid-1990s (Schena et al., 1995). Microarrays/chips are created by the immobilization of molecules in discrete spatial locations (**Figure 1**), normally at high density, on a solid surface. Depending on the nature of the molecule immobilized on the physical substrate, there are different types of microarrays, such as DNA, protein, glycan, small molecule, cell, or even tissues (Berard et al., 2012). The power of microarrays lies in their unprecedented capacity to simultaneously interrogate tens to hundreds of thousands of immobilized probes.

Microchips have been used extensively to study viral biology, many times addressing ecologically relevant questions. Most of these studies are circumscribed to the analysis of isolated viruses or to issues related to the role of viruses in human health and would therefore not fall within the scope of this mini review. However, many of the examples of microarray applications in the field of human virology could be (or have been) "upscaled" from individual populations to the whole viral community present in the samples under study. A similar rationale boosted the birth of microbial metagenomics, which benefited from the sequencing of the human genome.

The first uses (**Table 1**) of microarray tools in microbial viral ecology focused on the description of the diversity and dynamics of viral assemblages in environmentally relevant isolated viralhost pairs. In these examples, DNA probes encompassing the complete viral genomes were immobilized on a chip and used to interrogate their expression within infected cells or to conduct viral comparative genomic studies of isolated viruses.

Combined with metagenomics, microarrays have been used to analyze viral diversity and dynamics within natural samples (**Table 1**). For instance, Breitbart et al. (2008) constructed custom arrays with sequences retrieved from a week-old infant virome to follow the evolution of the viral assemblage during the first 2 weeks of the infant life. More recently, in a very elegant application of microarrays in viral ecology, Snyder et al. (2010) used a cellular CRISPR spacer-based microarray for detection of viruses in *Archaea*-dominated acidic hot springs samples, taken from Yellowstone National park (USA). The idea behind the design is that the spacers contained within the direct repeat units in the CRISPR loci of cellular genomes often derive from viruses that infected the cell in the past. Therefore, the spacer sequences are in a way a record of the infection history of the cell. These sequences were retrieved from environmental cellular metagenomes (by *in silico* search and PCR amplification) and immobilized as probes into a chip to interrogate the viral assemblage. This array proved to be useful both in the detection of unknown viruses and in the monitoring of viral dynamics within the community.

Nowadays, given the lower costs and relative ease of Next Generation Sequencing techniques, microarrays are no longer needed for addressing the issues mentioned above. However, this approach still offers some unparalleled possibilities, as described below.

#### **METAVIROTRANSCRIPTOMIC STUDIES**

Array-based viral metatranscriptomic studies overcome the fact that there is no straightforward, unambiguous way to extract viral mRNA from bulk mRNAs and make the sequencing of all mRNAs unnecessary. Thus, using the same terms presented in



*aThere are many examples on the study of phage-bacteria interactions related to disease and food production that are not considered here. All the examples provided use different kinds of DNA microchip, except that of Barr et al. which uses a glycoarray (see text).*

*bAlthough this study could not be strictly considered as "microbial viral ecology," we have included it here because it is a widely used model of interaction between phage and bacterial host.*

*cGene Transfer Agents.*

Allen and Wilson (2008) it no longer seems likely that microarrays will be used to ask "who is there?" in the field of viral ecology, although they can definitely help to answer the ecologically relevant questions "what is a particular component of the community doing and how is it doing it?"

Kunin et al. (2008) analyzed viral and microbial assemblages in two geographically distant sludge bioreactors (one in the USA, and one in Australia) that were enriched with a single bacterial species *Candidatus* Accumulibacter phosphatis (CAP). Comparison of viral and microbial metagenomes suggested that phages were active in the sludge ecosystem. To confirm this, the US sludge was monitored at three time points spanning 3 months using expression arrays constructed from selected predicted genes from both phage and bacterial metagenomes. Hybridization of cDNA from the three time-course points with the array indicated that a large number of viral genes were highly expressed and that some phages persisted for long periods in the sludge. These results indicated that the bacterial community in the sludge was under "persistent local predation pressure." Although microarrays had been used previously for studying the activity of viruses infecting humans, this is the first example of a metatranscriptomic approach specifically targeting viral mRNAs since, normally, in metatranscriptomic studies both viral and microbial mRNA are analyzed together and separated only *in silico*. As the authors state, this work illustrated the value of combining sequence and gene expression data from the bacterial and viral fractions to address ecologically relevant questions.

Applying a similar rationale, we analyzed the viral expression in one hypersaline environment using a metatranscriptomic approach in which clones from a metaviromic library were immobilized on a microarray and used as probes against total mRNA extracted from the hypersaline community (Santos et al., 2011). The immobilized metaviromic library had been prepared previously by cloning sheared viral DNA into plasmids. The metaviromic sequences were analyzed and classified according to their putative hosts based on GC content and dinucleotide frequency analysis, and thus each clone in the array corresponded to a genomic fragment of a virus tentatively assigned to a microbial host. The array was then hybridized with bulk mRNA extracted from the natural assemblage (previously converted into cDNA), which contained all the RNA from the microbial genomes expressed at the time of sampling as well as the viral RNA produced during infection. We found that the viral groups that had the highest hybridization signal in the array (i.e., highest expression levels compared to the rest of immobilized viruses) were those related to high GC content haloarchaea and *Salinibacter* representatives, which were minor components in the environment. In addition, we analyzed the changes in expression of the immobilized viruses when the natural samples were subjected to two stress conditions (ultraviolet radiation and dilution). More specifically, we also analyzed the expression of some viral genome fragments carrying non-synonymous single nucleotide polymorphisms (SNPs) and saw that differences in sequences could be related to changes in expression levels under various analyzed conditions. We suggested that the viral assemblages could include very closely related viruses, for which we proposed the term "ecoviriotypes," which would respond differentially to changes in environmental parameters.

## **ASSIGNING VIRUSES TO HOSTS**

Quite recently, we combined microarrays with *single cell genomics* (SCG) to assign virus-host pairs within natural uncultured microbial communities without the need for culture techniques (Martínez-García et al., in press). We investigated virus-host pairs within an *Archaea*-dominated hypersaline sample. First, we constructed a fosmid viral metagenomic library. This cloning approach was chosen because the optimum insert size for fosmids (i.e., between 30 and 45 kb) corresponded to the size of most haloviral genomes detected in the analyzed samples. Then, individual cloned viral genomes were immobilized in an array or "virochip" (which, in a way, resembles the one designed by Wang et al., 2002, that contained probes targeting different families of human pathogenic viruses). In parallel, genomic contents of uncultured individual cells present in the sample were retrieved by fluorescence activated cell sorting and their genomes amplified by multiple displacement amplification. Finally, individual cell genomes were hybridized against the virochip. Infected cells and immobilized viruses yielding positive hybridization were then sequenced and characterized. With this new approach, we described the first uncultured virus infecting the ubiquitous, uncultured Nanohaloarchaeal group in hypersaline environments, demonstrating the usefulness of combining these two high-throughput technologies (microarrays and SCG).

The advances of SCG as a mean to decipher genetic information of uncultured cells is spurring the development of *single virus genomics* (Allen et al., 2011) for unveiling genomic viral diversity in nature, one virus at a time. So, following the same rationale as our above mentioned study, massive microarrays could be constructed with thousands of immobilized single viral and prokaryote genomes, as a way to assign viruses to hosts in the majority of uncultured microbes, disentangling the complex virus-host network interactions within uncultured assemblages.

## **PROTEIN AND GLYCAN ARRAYS**

*Protein arrays* can be either analytical (to check for the presence of given proteins in the analyzed samples) or functional (to query for properties of the immobilized proteins) (Uzoma and Zhu, 2013). They can range from immobilized peptides to more complex systems as yeast two hybrid or phage display arrays. Analytical protein microarrays have been used for biomarker detection in virus infection, in the search for novel serological biomarkers in disease, and in the study of lectin-glycan interactions to characterize mammalian cell envelopes. Functional protein microarrays can be used for identifying protein-protein, protein-lipid, proteinantibody, protein-small molecule, protein-DNA (transcription factors), protein-RNA, protein-lectin interactions; for identifying substrates or enzymes corresponding to different modifications; and for profiling the immune response (reviewed in Hsu and Mahal, 2009; Ben-Ari et al., 2013; Uzoma and Zhu, 2013).

One specific type of protein microarrays uses antibodies as probes. They have been widely used in the diagnosis and monitoring of disease but also have applications that extend beyond the field of biomedicine. For instance, Parro et al. (2011) designed and built a series of instruments called SOLID (for "Signs Of Life Detector") for the automatic *in situ* detection and identification of substances that can be applied in the search for life in extreme environments and planetary exploration (Parro et al., 2011). This kind of microarrays are multiplex immunosensors with antibodies targeting polymeric biomolecules such as lipoteichoic acids, peptidoglycan, DNA, exopolysaccharides, proteins or whole cells, and can be used for biomarker and community profiling in environmental samples or for veterinary and biomedical applications (Palacín et al., 2012).

Although, to the best of our knowledge, protein microarrays have not been used to study phage-microbe interactions, they have been used to address complex aspects of the interaction between human viruses and their target cells. Therefore, why couldn't they be used to study the interactions of phages with their host bacterial communities? For instance, virion proteins immobilized on microarrays could be used to "fish" host envelope components involved in phage recognition. Furthermore, immunosensors constructed with antibodies against different phages could be used for monitoring the presence of these phages in environmental samples.

Finally, *glycan* (or carbohydrate) arrays allow the evaluation, in a high-throughput manner, of interactions between carbohydrates and proteins (including antibodies), viruses and cells (see Liang and Wu, 2009; Rillahan and Paulson, 2011). These arrays can be used to probe hundreds of receptor-ligand interactions in one experiment. They have been applied to the study of the activity of carbohydrate modifying enzymes and the roles of glycan in the detection of diseases such as AIDS and other viral and infectious diseases and in vaccine development. These chips can be built either by immobilizing known glycans or with so-called shotgun glycomics, which uses bulk glycans extracted from the cells of interest (Song et al., 2011).

In 2004, Blixt et al. (2004) published the description of an array of more than 200 carbohydrates representing major glycan structures of glycoproteins and glycolipids. This "glycochip" can be used for profiling the specificity of a diverse range of mammalian, plant, viral and bacterial glycan binding proteins (lectins), as well as of antibodies and intact viruses. Currently, the Consortium for Functional Glycomics has available an array with 610 diverse glycans found in mammalian cells that can be used for a wide range of applications within the biomedical field, such as the profiling of immune response to infection or the glycomics of several diseases of different etiology (www*.*functionallycomics*.* org). However, applications of the array extend beyond strictly disease-oriented analyses and have been used, for instance, to unveil interactions between the gut microbiota and the host environment (Garrido et al., 2011).

In a recent pioneering work within the field of phage ecology (Barr et al., 2013), the glycoarray was applied to the analysis of the interactions of bacteriophage and gut mucus. These authors demonstrated that mucus surfaces in metazoans were enriched in bacteriophages compared to the surrounding environments and that this enrichment occurred by binding between mucin glycopetides and immunoglobulin-like domains exposed on the capsid surfaces. Based on the results of their study, the authors proposed the bacteriophage adherence to mucus (BAM) model according to which phages provide a non-host derived antimicrobial defense on the mucosal surfaces of diverse metazoan hosts.

The presence of Ig-like proteins has been found in structural protein of approximately one fourth of sequenced tailed phages where they were proposed to facilitate adherence to bacterial cell surfaces during infection (Fraser et al., 2007). Given the relevance of carbohydrates in the phage-host recognition, it seems feasible that this approach may also be useful within the field of phage biology. For instance, in the above mentioned example of Kunin et al. (2008) the CAP strains detected in the bioreactors in the USA and Australia, although very closely related, presented highly variable extracellular polymeric substances (EPS) gene cassettes. EPS can act as a first line of defense against phage predation and, conversely, lytic bacteriophages can encode strain-specific EPS degrading polysaccharases (Sutherland, 2001). Conversely, carbohydrate binding and lectin-like domains have frequently been found within the metaviromic islands in marine phages (see the paper by Mizuno et al., 2014). Carbohydrate microarrays could then be a very powerful tool to explore this kind of glycan-mediated phage-microbe interaction.

## **FINAL REMARKS**

**Figure 1** summarizes all the applications of the different kinds of microarrays presented here. Our point has not been to convince readers to use microarrays but to illustrate the power of the combination of several high throughput techniques to unveil relevant aspects of virus-microbe interactions in the past and, hopefully, in the future. We hope that this goal has been accomplished.

## **ACKNOWLEDGMENTS**

Our current studies with viral microarrays are supported by projects CGL2012-39627-C03-01 (to Josefa Antón) and AYA2011-24803 (to Víctor Parro) of the Spanish Ministry of Science and Innovation, which are co-financed with FEDER support from the European Union.

## **REFERENCES**


**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.

*Received: 10 February 2014; accepted: 16 June 2014; published online: July 2014. Citation: Santos F, Martínez-García M, Parro V and Antón J (2014) Microarray tools to unveil viral-microbe interactions in nature. Front. Ecol. Evol. 2:31. doi: 10.3389/ fevo.2014.00031 07*

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

*Copyright © 2014 Santos, Martínez-García, Parro and Antón. This is an openaccess 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.*