It is widely known that Mollicutes have limited biosynthetic capabilities (Fraser et al., 1995; Sasaki et al., 2002; Vasconcelos, 2005); these are primarily for energy acquisition, with synthetic pathways being considerably reduced or absent. As such, many nutrients must be obtained from external sources (e.g., the host organism). Metabolic and genomic analyses have been carried out to elucidate the principles of metabolic regulation and adaptation to environmental conditions by Mollicutes (Lazarev et al., 2011; Kube et al., 2014; Vanyushkina et al., 2014). The *Acholeplasma* genome encodes a wide variety of ATP-binding cassette (ABC) transporters, and most genetic modules are implicated in the conversion of intermediates and are linked to ABC transporter genes for the uptake of oligopeptides and amino acids (Kube et al., 2014). This demonstrates that *Acholeplasma* are capable of incorporating oligopeptides and amino acids released from dead cells and their components. However, this raises questions about the assimilation of labile carbon sources by heterotrophic bacteria in slurry, such as why only the rRNA recovered from *Acholeplasma* was 13C-labeled despite most heterotrophic bacteria having access to dead cells and their components, and why 13C-labeled rRNA was recovered after a very short (2-h) incubation period. *Acholeplasma* has the capacity to generate dNTPs from adenosine, guanosine, uracil, and thymine (Kube et al., 2014), and *A. laidlawii* was reported to incorporate nucleic acid precursors during the first 4 h of incubation in medium (McIvor and Kenny, 1978). The detection of 13C-labeled rRNA within a 2-h incubation period leads us to speculate that *Acholeplasma* incorporated 13C-labeled nucleic acid precursors released by dead *E. coli* cells. Although most heterotrophic bacteria may incorporate 13C-labeled dead cells and their components, a longer incubation time would be required to obtain 13C-labeled rRNA from these bacteria, with the exception of *Acholeplasma*.

*Acholeplasma axanthum*, the closest relative to metabolically active bacteria identified in our study, has been isolated from a variety of habitats and hosts, including the lung of swine (Stipkovits et al., 1974). The genotypic patterns of these strains differed markedly from one another, suggesting considerable heterogeneity within the same species (Razin et al., 1983). Growth in varied environments requires regulatory switches and cross-talk between metabolic pathways to allow rapid adaptation to environmental changes in nutritional flux. The *A. laidlawii* genome contains genes encoding components of signal transduction pathways and the SOS response system, which are closely related to the regulation of gene expression and mutagenic response to stress (Lazarev et al., 2011). These functions allow *A. laidlawii* to adapt to various environmental conditions. It is unclear whether the *Acholeplasma* spp. detected in this study originated from the body of pigs or from the manure slurry itself; nonetheless, the results suggest that

# REFERENCES

Aoyagi, T., Kimura, M., Yamada, N., Navarro, R. R., Itoh, H., Ogata, A., et al. (2015). Dynamic transition of chemolithotrophic sulfur-oxidizing bacteria inresponse to amendment with nitrate in deposited marine sediments. *Front. Microbiol.* 6:426. doi: 10.3389/fmicb.2015.00426

Linking the phylogeny and function of microorganisms is one of the most fundamental challenges in microbial ecology. Generating defined mixed cultures of isolated microorganisms is one promising approach to achieve this end (Kato et al., 2005, 2014). However, in natural and engineered environments, unidentified microbial components may nonetheless play critical ecophysiological roles. In this study, we identified the first-toarrive scavengers in aerated pig manure slurry by 13C-labeled bacterial biomass probing of rRNA. *Acholeplasma* spp. were identified for the first time as a major heterotroph that assimilates the dead cells and their components in pig manure slurry. Substrate competition has been found in many types of biological processes, in which physiologically similar microorganisms coexist and fill ecological niches owing to the difference in the substrate availability and affinity (Hori et al., 2006; Aoyagi et al., 2015). Due to the limited biosynthetic capabilities, free-living *Acholeplasma* spp. most likely require biomolecules released from dead bacterial cells as substrates. Actually, they were capable of incorporating the dead cells and their components more rapidly than the other microbial populations, which allowed the survival and prevalence of *Acholeplasma* spp. in pig manure slurry. These findings provide deeper understanding of carbon flow mediated by microorganisms in copiotrophic environments and unveil one of the mechanisms underlying the microbial ecosystem development.

# AUTHOR CONTRIBUTIONS

DH designed the research. DH and TH wrote the manuscript. DH, TA, and TH carried out the experiments and analyzed the data.

# ACKNOWLEDGMENTS

This work was supported by a Japan Society for the Promotion of Science KAKENHI (grant number 21580338). We thank Dr. Satoshi Tobita and Dr. Kenji Fujieda for δ13C analysis of slurry samples by isotope ratio mass spectrometry.

# SUPPLEMENTARY MATERIAL

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

Atobe, H., Watabe, J., and Ogata, M. (1983). *Acholeplasma parvum*, a new species from horses. *Int. J. Syst. Bacteriol.* 33, 344–349. doi: 10.1099/00207713- 33-2-344

Chen, J. G., Lou, D., and Yang, J. F. (2011). Isolation and identification of *Acholeplasma* sp. from the mud crab. Scylla serrata. *Evid. Based Complement. Alternat. Med.* 2011, 209406. doi: 10.1155/2011/209406

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

*Copyright © 2015 Hanajima, Aoyagi and Hori. 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 effects of elevated CO2 concentration on competitive interaction between aceticlastic and syntrophic methanogenesis in a model microbial consortium

## *Souichiro Kato1,2,3 \*, Rina Yoshida4 , Takashi Yamaguchi 4 , Tomoyuki Sato1, Isao Yumoto1 and Yoichi Kamagata1,2*

1 Bioproduction Research Institute, National Institute of Advanced Industrial Science and Technology, Sapporo, Japan

2 Division of Applied Bioscience, Graduate School of Agriculture, Hokkaido University, Sapporo, Japan

3 Research Center for Advanced Science and Technology, The University of Tokyo, Tokyo, Japan

4 Department of Civil and Environmental Engineering, Nagaoka University of Technology, Nagaoka, Japan

#### *Edited by:*

Shin Haruta, Tokyo Metropolitan University, Japan

#### *Reviewed by:*

Hans Karl Carlson, University of California, Berkeley, USA Spyridon Ntougias, Democritus University of Thrace, Greece Shawn E. McGlynn, Tokyo Metropolitan University, Japan

#### *\*Correspondence:*

Souichiro Kato, Bioproduction Research Institute, National Institute of Advanced Industrial Science and Technology, 2-17-2-1 Tsukisamu-Higashi, Toyohira-ku, Sapporo, Hokkaido 062-8517, Japan e-mail<EMAIL_ADDRESS>
Investigation of microbial interspecies interactions is essential for elucidating the function and stability of microbial ecosystems. However, community-based analyses including molecular-fingerprinting methods have limitations for precise understanding of interspecies interactions. Construction of model microbial consortia consisting of defined mixed cultures of isolated microorganisms is an excellent method for research on interspecies interactions. In this study, a model microbial consortium consisting of microorganisms that convert acetate into methane directly (Methanosaeta thermophila) and syntrophically (Thermacetogenium phaeum and Methanothermobacter thermautotrophicus) was constructed and the effects of elevated CO2 concentrations on intermicrobial competition were investigated. Analyses on the community dynamics by quantitative RT-PCR and fluorescent in situ hybridization targeting their 16S rRNAs revealed that high concentrations of CO2 have suppressive effects on the syntrophic microorganisms, but not on the aceticlastic methanogen. The pathways were further characterized by determining the Gibbs free energy changes (-G) of the metabolic reactions conducted by each microorganism under different CO2 concentrations. The -G value of the acetate oxidation reaction (T. phaeum) under high CO2 conditions became significantly higher than −20 kJ per mol of acetate, which is the borderline level for sustaining microbial growth.These results suggest that high concentrations of CO2 undermine energy acquisition ofT. phaeum, resulting in dominance of the aceticlastic methanogen. This study demonstrates that investigation on model microbial consortia is useful for untangling microbial interspecies interactions, including competition among microorganisms occupying the same trophic niche in complex microbial ecosystems.

**Keywords: model consortia, methanogenesis, acetate, thermodynamics, CO2 concentration**

#### **INTRODUCTION**

In natural and engineered environments, many species of microorganisms coexist by interacting with each other. Comprehension of interspecies interactions is essential for describing the features of complex microbial ecosystems, and competition among microorganisms occupying similar trophic niches is a conventional and significant aspect of such interspecies interaction. Coexistence of multiple microorganisms with similar trophic niches is regarded as one of the major factors to confer functional stability and resiliency on microbial ecosystems (Loreau et al., 2001; Deng, 2012). It is important to grasp how the population of each microorganism changes depending on a specific environmental disturbance. Most microbial ecological research has assessed the effects of specific environmental factors on competitive interactions among multiple microbial species by observing the transition of abundances of each microorganism responding to environmental disturbances. Although this approach has produced

many excellent outcomes, existence of non-target microorganisms and uncontrollable environmental factors in the systems often hamper precise understanding of the effects of specific environmental factors on the competitive interactions among target microorganisms.

Construction of microbial model consortia, in which interspecies interactions in ecosystems are reproduced by defined co-culture of isolated microorganisms, is appreciated as a worthwhile method to investigate microbial interactions (Haruta et al., 2009; De Roy et al., 2014; Großkopf and Soyer, 2014). For instance, the complex phenomenon of bacterial competition as being similar to rock-paper-scissors among colisin-producing, colisin-resistant, and colisin-sensitive strains was untangled by constructing model co-culture systems (Kerr et al., 2002; Nahum et al., 2011). Kato et al. (2005, 2008) constructed model microbial consortia composed of 4–5 bacterial strains, in which all members stably coexisted for long period of time, and demonstrated

that existence of both positive and negative interspecies interactions among the members make these consortia stable. The construction of model consortia is a specific and beneficial feature of microbiological research fields, which will also be effective for proof-of-concept studies for theories in the field of macro-ecology (Haruta et al., 2009, 2013).

Methanogenesis from organic compounds is a complex microbial process accomplished by catabolic interactions among different trophic levels of microorganisms (Schink, 1997; Jones et al., 2008; Kato and Watanabe, 2010). Among the sequential biodegradation processes, acetate is the most important intermediary metabolite (Schink,1997). Methanogenic acetate degradation proceeds by either aceticlastic methanogenesis or syntrophic acetate oxidation. The aceticlastic pathway is solely mediated by aceticlastic methanogens (Jetten et al., 1992). On the contrary, syntrophic acetate oxidation pathway requires cooperative interactions of two different types of microorganisms: acetate is first oxidized to H2 and CO2 by syntrophic acetate-oxidizing bacteria (SAOB), and then hydrogenotrophic methanogens convert the products to CH4 (Zinder and Koch, 1984). As the acetate oxidation reaction is endoergonic under the standard conditions and is feasible only under extremely low H2 partial pressure, acetate oxidation by SAOB requires H2 elimination by hydrogenotrophic methanogens (Karakashev et al., 2006; Hattori, 2008). These two different acetate-degrading methane-producing pathways and organisms involved can co-exist, but diverse environmental factors, such as temperature, pH, salinity, toxic compounds, and concentrations of substrates determine one pathway and organisms to dominate over the other (Nüsslein et al., 2001; Shigematsu et al., 2004; Karakashev et al., 2006; Hao et al., 2013; Kato et al., 2014).

In our previous studies, we demonstrated that the syntrophic pathway is the dominant methanogenic acetate degradation pathway in underground, thermophilic petroleum reservoirs (Mayumi et al., 2011). We further demonstrated that aceticlastic pathway becomes dominant under high CO2 concentrations, which mimicked carbon capture and storage field conditions (Mayumi et al., 2013), whereas syntrophic acetate oxidation dominated over aceticlastic reactions under low CO2 concentrations. Since CO2 is either substrate or product of aceticlastic methanogenesis, acetate oxidation, and hydrogenotrophic methanogenesis, high CO2 concentration alters the thermodynamics of each methanogenic reaction, which may cause the observed transition between syntrophic and aceticlastic methanogenic pathways. However, all the data were based on the analyses of complex microbial communities in field samples thus many other factors that affect the community shift could not be ruled out.

In the present study, the effect of CO2 concentrations on methanogenic microorganisms were assessed by using a defined inorganic medium and a defined methanogenic consortium which is comprised of three organisms, i.e., SAOB, hydrogenotrophic methanogen and aceticlastic methanogen, namely, which contains two different acetate-degrading methanogenic pathways. The experiments allowed to precisely show the CO2 concentrations to be a crucial factor affecting the dominance of respective pathways and organisms.

#### **MATERIALS AND METHODS**

#### **MICROORGANISMS AND CULTURE CONDITIONS**

*Methanosaeta thermophila* DSM6194T (Kamagata and Mikami, 1991) and *Thermacetogenium phaeum* DSM12270T (Hattori et al., 2000) were obtained from the Deutsche Sammlung von Mikroorganismen und Zellkukturen GmbH (Braunschweig, Germany). *Methanothermobacter thermautotrophicus* strain TM was isolated from a thermophilic anaerobic methanogenic reactor in Japan (Hattori et al., 2000). Routine cultivations were conducted at 55◦C with 68-ml capacity serum vials containing 20 ml of a bicarbonate-buffered inorganic medium (pH 7.0; Kato et al., 2014) under an atmosphere of N2-CO2 [80/20 (v/v)] without shaking. Pyruvate (40 mM) or 200 kPa H2-CO2 [80/20 (v/v)] was supplemented as energy and carbon sources for the pure cultures of *T. phaeum* and *Methanothermobacter thermautotrophicus*, respectively. Sodium acetate (40 mM) was utilized as an energy and carbon source for the pure culture of *Methanosaeta thermophila*, the defined co-culture of *T. phaeum* and *Methanothermobacter thermautotrophicus*, and the tri-culture of the three strains. The tri-culture was constructed by simultaneously inoculating 1 and 2 ml of the early-stationary phases of pure culture of *Methanosaeta thermophila* and the defined coculture of *T. phaeum* and *Methanothermobacter thermautotrophicus,* respectively, into the 20 ml of the medium. Although the long term stability of the tri-culture was not been tested, coexistence of the three microorganisms in the batch culture was confirmed.

#### **CULTURES WITH DIFFERENT CO2 CONCENTRATIONS**

Three culture conditions were prepared to examine the effects of CO2 concentrations on the microorganisms. For each condition, the media were supplemented with different concentrations of sodium bicarbonate and the gas phases were replaced with N2/CO2 mixed gas with different volume ratios, as described in **Table 1**. The medium was bubbled with the respective deoxygenated gas with 100 ml min−1 for 5 min and immediately capped with a butyl rubber stopper and an aluminum cap. The medium pH was adjusted to 7.0 by adding 1N NaOH solution before the cultivation, and the fluctuation of pH value throughout the cultivation was less than 0.2. For pH measurement, 100 μl of the medium was sampled with syringes and the pH value was determined using a compact pH meter B-212 (Horiba). The concentration of CO2 in the aqueous phase [*caq* (M)] was calculated according to Henry's law (*caq* = *kp*), where *k* is the Henry's low constant (0.019 for CO2 at 55◦C) and *p* is the partial pressure of CO2 in the gas phase (atm). Then the bicarbonate concentrations were calculated based on the equilibrium formula (H2CO3 = H+ + HCO− 3 ) with the

equilibrium constant of 4.47 × 107. The culture experiments were conducted in triplicate and the student's *t*-test was used for the statistical analyses.

Growth of *Methanosaeta thermophila* and *Methanothermobacter thermautotrophicus* in pure and mixed cultures was determined by measuring methane production. Growth of *T. phaeum* pure culture was determined by measuring acetate production from pyruvate. The partial pressure of CH4 was determined using a gas chromatograph GC-2014 (Shimadzu) as described previously (Kato et al., 2014). The partial pressure of H2 was determined using a trace reduction gas analyzer TRA-1000 (ACE Inc.) according to the manufacturer's instruction. The concentrations of organic acids were determined using a high performance liquid chromatography (D-2000 LaChrom Elite HPLC system, HITACHI) equipped with Aminex HPX-87H Ion Exclusion column (BIO-RAD) and L2400 UV detector (HITACHI).

#### **FLUORESCENT** *IN SITU* **HYBRIDIZATION (FISH)**

Microbial cells of the tri-cultures in the early stationary phases were collected by centrifugation, fixed with 4% paraformaldehyde in phosphate buffered saline (PBS; 137 mM NaCl, 2.7 mM KCl, 8.1 mM Na2HPO4, 1.5 mM KH2PO4, pH 7.2) and left for 6 h at 4◦C. The samples were washed three times with PBS, immobilized on glass slides, and dehydrated by successive passages through 50, 70, 80, 90, and 100% ethanol (3 min each). The following oligonucleotide probes complementary to specific regions of 16S rRNA were utilized for hybridizations: (i) Alexa488 labeled EUB338, specific for the domain *Bacteria* (Amann et al., 1990) and (ii) TexRed-labeled ARCH917, specific for the domain *Archaea* (Loy et al., 2002), (iii) Alexa594-labeled MSMX860, specific for the order *Methanosarcinales* (Raskin et al., 1994), and (iv) Alexa488-labeled MB311, specific for the order *Methanobacteriales* (Crocetti et al., 2006). Hybridizations were performed at 46◦C for 3 h with hybridization buffer (0.9 M NaCl, 0.1 M Tris-HCl, pH 7.5) containing 5 ng μl −1 of each labeled probe. The specificity of each probe was confirmed by FISH observations using pure cultures of the three microorganisms used in this study even with the hybridization buffer not containing formamide. The washing step was done at 48◦C for 30 min with washing buffer (0.2 M NaCl, 0.1 M Tris-HCl, pH 7.5). The samples hybridized with the probes were observed with a fluorescent microscope Provis AX70 (Olympus).

#### **QUANTITATIVE RT-PCR (qRT-PCR)**

Microbial cells were harvested from the mid-logarithmic phases by centrifugation at 10,000 X *g* and 4◦C. Total RNA was isolated using ISOGEN II reagent (Nippon Gene, Japan) combined with a beadbeating method, as described previously (Kato et al., 2014). Total RNA was purified using an RNeasy Mini kit (Qiagen) with DNase treatment (RNase-free DNase set, Qiagen) as described in the manufacturer's instructions. The purified RNA was spectroscopically quantified using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies). The PCR primers used for quantitative RT-PCR (qRT-PCR) were designed with Primer3 software (http://simgene.com/Primer3) and are listed in **Table 2**. Quantification of 16S rRNA copy numbers in the defined mixed culture

**Table 2 | Quantitative RT-PCR (qRT-PCR) primers designed and used in this study.**

were performed by one-step real-time RT-PCR using a Mx3000P QPCR System (Stratagene) and RNA-direct SYBR Green Realtime PCR Master Mix (Toyobo) as described previously (Kato et al., 2014). At least three biological replicates were subjected to qRT-PCR analysis, and at least two separate trials were conducted for each sample. Standard curves were generated with serially diluted PCR products (103–108 copies ml−1) amplified using the respective primer sets and were used to calculate the copy number of rRNA in the total RNA samples.

#### **RESULTS AND DISCUSSION**

#### **EFFECTS OF CO2 CONCENTRATIONS ON THE MODEL METHANOGENIC CONSORTIUM**

As the model consortium performing methanogenic acetate degradation, we utilized a defined mixed culture of an aceticlastic methanogen (*Methanosaeta thermophila*), a hydrogenotrophic methanogen (*Methanothermobacter thermautotrophicus*), and a SAOB (*T. phaeum*; **Table 3**). These microbial species were originally isolated from a thermophilic methanogenic digester (Kamagata and Mikami, 1991; Hattori et al., 2000) and are regarded as representative species for the methanogenic acetate degradation reactions that occur in various natural environments such as high-temperature petroleum reservoirs (Pham et al., 2009; Mayumi et al., 2011, 2013) and thermophilic methanogenic digesters (Sekiguchi et al., 1998; McHugh et al., 2003; Hori et al., 2011).

To adequately assess the effects of CO2 concentration itself, media with different supplementation of CO2/HCO− 3 were prepared (**Table 1**). The initial concentrations of total CO2/HCO− 3 in the cultures, designated as [CO2]initial, were 5.0, 50.7, or 113.4 mmol l−1. The model consortium composed of *Methanosaeta thermophila*, *Methanothermobacter thermautotrophicus*, and *T. phaeum* was cultivated under the three different [CO2]initial conditions to evaluate their methanogenic acetate degradation abilities (**Figure 1**). A stoichiometric production of CH4 from acetate in a 1:1 molar ratio was observed in all culture conditions tested. Both acetate consumption and CH4 production rates slightly decreased with increasing the [CO2]initial (**Figures 1A,B**). Interestingly, the partial pressure of H2, which is an important intermediate of syntrophic

#### **Table 3 |The metabolic reactions and the respective standard Gibbs free energy changes (G◦') of the microorganisms utilized in this study.**

aThe -G◦ values were calculated according to the reference (Thauer et al., 1977).

*Methanothermobacter thermautotrophicus***, and** *Methanosaeta thermophila***.** Time courses of acetate **(A)**, CH4 **(B),** and H2 **(C)** concentrations during cultivation on acetate with different [CO2]initial are shown. Data are presented as means of three independent cultures, and error bars represent SDs.

acetate degradation, significantly decreased with increasing the [CO2]initial (**Figure 1C**). This observation suggests that syntrophic methanogenic microorganisms are influenced by elevated CO2 concentrations.

To assess the influence of the elevated CO2 concentrations on each methanogenic pathway, the relative abundances of each microorganism in the exponentially growing cultures of the model consortium with the different [CO2]initial were evaluated by FISH and qRT-PCR analyses. The qRT-PCR analysis clearly demonstrated the decrease of the abundances of *Methanothermobacter thermautotrophicus* and *T. phaeum* in the higher [CO2]initial cultures (**Figure 2**). The FISH analysis also demonstrated that the relative abundances of *Methanothermobacter thermautotrophicus* and *T. phaeum* in the cultures with higher CO2 concentrations are significantly lower than those in the low CO2 cultures (**Figure 3**). These results indicate that the syntrophic methanogenic pathway is more strongly influenced by the elevation of CO2 concentrations compared to the aceticlastic pathway.

**consortium with different [***-***CO2]initial.** The 16S rRNA copy numbers of each microorganism in the mid-exponential phases were determined by the qRT-PCR analysis. The abundance of each microorganism was normalized against those of the cultures with [CO2]initial of 5.0 mmol l−1, and plotted against the respective [CO2]initial values. Data are presented as the means of three independent cultures, and error bars represent SDs.

#### **EFFECTS OF CO2 CONCENTRATIONS ON THE ACETICLASTIC AND SYNTROPHIC PATHWAYS**

To confirm the differences in the suppressive effects of elevated CO2 concentrations on the two methanogenic pathways, the pure culture of *Methanosaeta thermophila* and the defined co-culture of *Methanothermobacter thermautotrophicus* and *T. phaeum* were separately cultivated in the media with the different [CO2]initial (**Figure 4**). The growth of *Methanosaeta thermophila* was barely affected by the elevated CO2 concentration: the methanogenic rate in the [CO2]initial of 113.4 mmol l −1 cultures decreased only about 10% compared to the cultures with [CO2]initial of 5.0 mmol l−1 (**Figures 4A,C**). On the contrary, the methanogenic rate of the syntrophic co-culture in the [CO2]initial of 113.4 mmol l−1 dropped to less than half of that in the cultures with [CO2]initial of 5.0 mmol l−1 (**Figures 4B,C**). These observations confirm the assumption that the syntrophic acetate degradation pathway is more susceptible to elevated CO2 concentrations than the aceticlastic pathway.

#### **EFFECTS OF CO2 CONCENTRATIONS ON THE PURE CULTURES OF** *Methanothermobacter thermautotrophicus* **AND** *T. phaeum*

One possible explanation for the suppressive effects of CO2 on the syntrophic methanogenesis is the susceptibility of

*Methanothermobacter thermautotrophicus* and/or *T. phaeum* to some environmental alterations induced by increased CO2 or to CO2 itself. To evaluate this possibility, pure cultures of *Methanothermobacter thermautotrophicus* and *T. phaeum* were cultivated in media with different [CO2]initial (**Figure 5**). No significant differences were observed for the growth of both *Methanothermobacter thermautotrophicus* and *T. phaeum* under the different CO2 conditions tested. These results suggest that elevated CO2

independent cultures and error bars represent SDs.

concentrations negatively affect the microbial activity only when *Methanothermobacter thermautotrophicus* and *T. phaeum* are in a syntrophic relationship.

#### **EFFECTS OF CO2 CONCENTRATIONS ON THE THERMODYNAMICS OF EACH REACTION**

The other possible explanation for the suppression of syntrophic methanogenesis by elevated CO2 concentration is alterations of thermodynamic conditions of each microbial reaction. A minimum energy required for biochemical energy conversion is estimated at around –20 kJ mol−1 (Schink, 1997), while some anaerobic microorganisms have been reported to thrive under more thermodynamically restricted conditions (Jackson and McInerney, 2002; Nauhaus et al., 2002). The value was estimated from the energetics of ATP formation (around −70 kJ mol−1 under the physiological conditions; Jetten et al., 1991; Tran and Unden, 1998) and the number of protons transported to ATP formation (between 3 and 4; Maloney, 1983; Stock et al., 1999). Since syntrophic methanogenesis from acetate

is one of the least exergonic microbial metabolisms (Schink, 1997), it is no wonder that only slight perturbations on the thermodynamics induce deteriorations of the syntrophic methanogenesis.

To evaluate the influences of elevated CO2 concentrations on the thermodynamic properties, -G values of metabolic reactions conducted by each microorganism in the model consortium were determined using the data-set of metabolite concentrations shown in **Figure 1**. The -G values of the aceticlastic methanogenesis conducted by *Methanosaeta thermophila* were not significantly influenced by the elevated CO2 concentrations (**Figure 6**). The average -G values during the logarithmic growth phase (day 2– 5) with the [CO2]initial of 5.0, 50.7 and 113.4 mmol l−1 were −47.7 ± 3.5, −44.9 ± 2.6, and −44.6 ± 2.0 kJ mol−1, respectively, which are substantially lower than the -G value required for microbial energy acquisition.

The -G values of the hydrogenotrophic methanogenesis catalyzed by *Methanothermobacter thermautotrophicus* were also largely not altered with different CO2 settings and were constantly lower than −20 kJ mol−1 (**Figure 6**). The average -G values during the logarithmic growth phases with the [CO2]initial of 5.0, 50.7, and 113.4 mmol l−1 were −24.6 ± 1.0, −27.2 ± 1.0, and −26.0 ± 1.2 kJ mol−1, respectively. Since CO2 is the substrate for hydrogenotrophic methanogenesis, lower -G values under the higher CO2 conditions are expected. However, the decrease in H2 partial pressures under the higher CO2 conditions (**Figure 1C**) compensates for the positive effects of increase in CO2 concentration.

On the contrary, elevation of CO2 concentrations significantly influenced the -G values of the acetate oxidation reaction performed by *T. phaeum* (**Figure 6**). While the average -G value during the logarithmic growth phases with the [CO2]initial of 5.0 mmol l−1 (−23.1 ± 2.7 kJ mol−1) was less than the borderline -G value of −20 kJ mol−1, those with the [CO2]initial of 50.7 and 113.4 mmol l−1 (−17.8 ± 1.3 and −18.7 ± 1.4 kJ

mol−1, respectively) exceeded the borderline. As acetate oxidation reaction produces 2 mol of CO2 from 1 mol of acetate, it is rational that this reaction is strongly influenced by the elevation of CO2 concentration. The decrease in the partial pressure of H2, the other metabolic product of acetate oxidation, is expected to compensate for the negative effects of increase in CO2. However, the decrease in H2 partial pressure would be limited by the minimum threshold for H2 consumption by *Methanothermobacter thermautotrophicus*. The minimum thresholds for H2 utilization by hydrogenotrophic methanogens have been reported as around 5–10 Pa (Lovley, 1985; Thauer et al., 2008). However, considering the energy required for active growth, H2 partial pressure of around 10–15 Pa observed in the increased CO2 conditions in this study may be the minimum H2 threshold for the syntrophic interaction. Actually, if the H2 partial pressure in the cultures with [CO2]initial of 113.4 mmol l−1 at the logarithmic growth phase (day 5) becomes 10 Pa, the -G value becomes > −20 kJ mol−1 (−19.7 ± 0.3 kJ mol−1) from the actual value of −25.1 ± 1.4 kJ mol−1 (with H2 partial pressure of 16.4 ± 1.7 Pa). These results clearly demonstrated that high concentrations of CO2 thermodynamically constrain the acetate oxidizing reaction, which results in the deterioration of syntrophic methanogenesis from acetate.

#### **CONCLUSION**

This is the first paper to evaluate the influence of elevated CO2 concentration on the two different methanogenic acetate degradation pathways, namely aceticlastic and syntrophic pathways, using a model microbial consortium. As expected from the observations based on *in situ* environments with complex microbial communities, high concentrations of CO2 suppressed the syntrophic pathway rather than the aceticlastic pathway. Thermodynamic calculations revealed that the acetate oxidation reaction is more intensely constrained by elevated CO2 concentrations. This study exemplified the importance of even slight changes in the -G values of microbial metabolisms in anaerobic biota. Furthermore, this study demonstrated that the construction of model microbial consortia is useful for assessing competitive interspecies interactions even in anaerobic, methanogenic environments.

#### **AUTHOR CONTRIBUTIONS**

Souichiro Kato, Tomoyuki Sato, and Yoichi Kamagata designed the research. Souichiro Kato, Rina Yoshida, Takashi Yamaguchi, Tomoyuki Sato, and Isao Yumoto carried out the experiments and analyzed the data. Souichiro Kato and Yoichi Kamagata wrote the manuscript.

#### **ACKNOWLEDGMENTS**

This work was supported by grants from the Japan Society for the Promotion of Science (JSPS). We thank Dr. Mia Terashima and Hiromi Ikebuchi for helpful comments by critical reading and technical assistance, respectively.

#### **REFERENCES**

Amann, R. I., Krumholz, L., and Stahl, D. A. (1990). Fluorescent-oligonucleotide probing of whole cells for determinative, phylogenetic, and environmental studies in microbiology. *J. Bacteriol.* 172, 762–770.

**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: 28 August 2014; accepted: 13 October 2014; published online: 30 October 2014.*

*Citation: Kato S, Yoshida R, Yamaguchi T, Sato T, Yumoto I and Kamagata Y (2014) The effects of elevated CO2 concentration on competitive interaction between aceticlastic and syntrophic methanogenesis in a model microbial consortium. Front. Microbiol. 5:575. doi: 10.3389/fmicb.2014.00575*

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

*Copyright © 2014 Kato, Yoshida, Yamaguchi, Sato, Yumoto and Kamagata. 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.*

# Recombination Does Not Hinder Formation or Detection of Ecological Species of Synechococcus Inhabiting a Hot Spring Cyanobacterial Mat

Melanie C. Melendrez 1 \*, Eric D. Becraft 1 , Jason M. Wood1 , Millie T. Olsen1 , Donald A. Bryant 2 , John F. Heidelberg3 , Douglas B. Rusch4 , Frederick M. Cohan5 and David M. Ward1

*1 Department of Land Resources and Environmental Science, Montana State University, Bozeman, MT, USA, 2 Department of Biochemistry and Molecular Biology, Pennsylvania State University, University Park, PA, USA, 3 Department of Biological Sciences, College of Letters, Arts and Sciences, University of Southern California, Los Angeles, CA, USA, 4 Informatics Group, J. Craig Venter Institute, Rockville, MD, USA, 5 Department of Biology, Wesleyan University, Middletown, CT, USA*

#### Edited by:

*Hiroyuki Futamata, Shizuoka University, Japan*

## Reviewed by:

*Kim Marie Handley, University of Chicago, USA Michael Kevin Watters, Valparaiso University, USA*

\*Correspondence: *Melanie C. Melendrez<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Systems Microbiology, a section of the journal Frontiers in Microbiology*

Received: *31 August 2015* Accepted: *21 December 2015* Published: *14 January 2016*

#### Citation:

*Melendrez MC, Becraft ED, Wood JM, Olsen MT, Bryant DA, Heidelberg JF, Rusch DB, Cohan FM and Ward DM (2016) Recombination Does Not Hinder Formation or Detection of Ecological Species of Synechococcus Inhabiting a Hot Spring Cyanobacterial Mat. Front. Microbiol. 6:1540. doi: 10.3389/fmicb.2015.01540* Recent studies of bacterial speciation have claimed to support the biological species concept—that reduced recombination is required for bacterial populations to diverge into species. This conclusion has been reached from the discovery that ecologically distinct clades show lower rates of recombination than that which occurs among closest relatives. However, these previous studies did not attempt to determine whether the more-rapidly recombining close relatives within the clades studied may also have diversified ecologically, without benefit of sexual isolation. Here we have measured the impact of recombination on ecological diversification within and between two ecologically distinct clades (A and B') of *Synechococcus* in a hot spring microbial mat in Yellowstone National Park, using a cultivation-free, multi-locus approach. Bacterial artificial chromosome (BAC) libraries were constructed from mat samples collected at 60◦C and 65◦C. Analysis of multiple linked loci near *Synechococcus* 16S rRNA genes showed little evidence of recombination between the A and B' lineages, but a record of recombination was apparent within each lineage. Recombination and mutation rates within each lineage were of similar magnitude, but recombination had a somewhat greater impact on sequence diversity than mutation, as also seen in many other bacteria and archaea. Despite recombination within the A and B' lineages, there was evidence of ecological diversification within each lineage. The algorithm Ecotype Simulation identified sequence clusters consistent with ecologically distinct populations (ecotypes), and several hypothesized ecotypes were distinct in their habitat associations and in their adaptations to different microenvironments. We conclude that sexual isolation is more likely to follow ecological divergence than to precede it. Thus, an ecology-based model of speciation appears more appropriate than the biological species concept for bacterial and archaeal diversification.

Keywords: population genetics, speciation, Synechococcus, cyanobacteria, multi-locus sequence typing, Ecotype Simulation, ecotype, recombination

# INTRODUCTION

For much of the last century, theories of the origin of species have revolved around genetic exchange (Mayr, 1942, 1963; Dobzhansky, 1951; Coyne and Orr, 2004). In highly sexual organisms such as animals and plants, genetic exchange is frequently seen as a cohesive force holding together the individuals within a species population, thus preventing divergence into irreversibly separate lineages (Cohan, 2013). The role of genetic exchange in bacteria and archaea, where recombination is not connected to reproduction (Levin and Bergstrom, 2000), has been vigorously debated. Recombination has been viewed by some as a rampant process that can erode the distinctness of nascent species (Doolittle and Zhaxybayeva, 2009); moreover, recombination is seen as a cohesive force binding a species together, as in the case of sexual species. Speciation theories based on Mayr's Biological Species Concept see the breaking of recombination between populations as a critical step in population divergence; in this view, populations can diverge ecologically and irreversibly only after their genetic exchange is blocked (Coyne and Orr, 2004).

On the other hand, theories of ecological speciation hold that ecological divergence is the critical step in the origin of diversity, and that recombination plays only a minor role in the subsequent fate of newly split, ecologically distinct populations (Van Valen, 1976; Mallet, 2008; Schluter, 2009). According to a theory developed long ago by Haldane, rare introduction of genes from one population to another cannot reverse the adaptive divergence between populations (Haldane, 1932). This is because natural selection against niche-specifying genes from other populations can limit these foreign genes to negligible frequencies, especially if the rate of recombination between populations is very low. The equilibrium frequency of a recombinant niche-specifying gene is equal to the rate of recombination between populations (cb ) divided by the intensity of selection (s) against the gene (Mallet, 2008; Wiedenbeck and Cohan, 2011). The equilibrium frequency of recombinant genes is thus extremely low, as recombination rates are not known to exceed the mutation rate by more than a factor of 10 (Vos and Didelot, 2009), and mutation rates are known to be very low in nearly all bacteria and archaea studied (about 10−6 per gene per generation; Drake, 2009). It can therefore be inferred that recombination rates are much too low to reverse adaptive divergence of microbial populations, and that recombination is unlikely to prevent ecological divergence among the most closely related bacteria and archaea (Wiedenbeck and Cohan, 2011).

The conflict between recombination-based and ecology-based theories of speciation has motivated several recent studies on the role of recombination in the evolution of bacterial and archaeal species. In a study of coexisting strains of Sulfolobus islandicus isolated from the same hot spring, Cadillo-Quiroz and colleagues found depressed rates of recombination between two major clades formed by these organisms compared to recombination within either clade (Cadillo-Quiroz et al., 2012). Inter-clade recombination was taken to be insufficient to prevent the divergence of the clades into easily distinguishable clusters, and the authors claimed this as support for the Biological Species Concept in archaeal speciation. Ellegaard et al. came to a similar conclusion based on comparison of genomes derived from closely related Wolbachia strains (97–98% 16S rRNA identity) isolated from the same host fly (Ellegaard et al., 2013). Shapiro et al. obtained similar results for closely related clades of marine Vibrio cyclitrophicus strains that may be ecologically distinct, based on their isolation from particles of different sizes (Shapiro et al., 2012). They proposed a model in which ecological specialization, based on horizontal genetic transfer, leads to occupancy of different habitats; this causes a decrease in recombination between the ecologically distinct populations, allowing the ecological populations to further diverge (Shapiro et al., 2012; see also Retchless and Lawrence, 2010; Polz et al., 2013). In all these studies the reduced recombination between the clades under investigation was taken as evidence that recombination must be lowered to allow irreversible divergence among close relatives into stably coexisting, ecologically distinct populations. However, these studies did not address whether ecological divergence may have also occurred within each clade, where recombination rates are higher.

Here we test whether ecological divergence can occur among organisms so closely related that recombination between them is not reduced. Our model system is the set of unicellular cyanobacteria (Synechococcus) inhabiting the microbial mats of Mushroom Spring, Yellowstone National Park. This hypothesis emerged from a series of studies based on variation found in single genetic loci (Ferris and Ward, 1997; Ramsing et al., 2000; Ferris et al., 2003; Melendrez et al., 2011). First, closely related 16S rRNA genotypes were distributed differently along the effluent channel (genotypes A," A', A, B', and B along a temperature gradient from ∼72–74◦C to ∼50◦C Ferris and Ward, 1997) and with depth in the photic zone (Ramsing et al., 2000), suggesting differential adaptations to temperature and light that were later confirmed by studying isolates with genotypes representative of in situ populations (Allewalt et al., 2006). Second, analyses based on the internal transcribed spacer separating 16S and 23S rRNA genes revealed that 16S rRNA sequence variation was insufficient to detect all ecologically distinct populations (Ferris et al., 2003; Ward et al., 2006), limiting the utility of 16S rRNAbased methods for detecting closely related ecotypes (e.g., Eren et al., 2013). Third, analyses of several more rapidly evolving protein-encoding genes revealed numerous closely related and ecologically distinct populations, the number of which depended on the depth of coverage and on the degree of sequence divergence of the gene used (Becraft et al., 2011; Melendrez et al., 2011). It is clearly important to appreciate that the amount of molecular divergence between these closely related species is very small (see Discussion).

Our aim in these studies has been to discover the fundamental "ecotypes" within the A and B' clades (originally demarcated as 16S rRNA clusters), where we define an ecotype as a group of organisms that is ecologically distinct from other ecotypes and whose members are ecologically interchangeable (Kopac and Cohan, 2011). More specifically, different ecotypes are expected to coexist indefinitely owing to their ecological distinctness, while lineages within an ecotype are too similar to be able to coexist indefinitely (Kopac et al., 2014). In our recent work we hypothesized the demarcations of putative ecotypes (PEs) from sequence clusters using a theory-based algorithm (Ecotype Simulation, based on the Stable Ecotype Model of species and speciation, Kopac et al., 2014; Cohan, 2016), which simulates the evolutionary dynamics of ecotype formation and the forces that purge sequence diversity within ecotypes (Koeppel et al., 2008). Unlike other studies that have used traditional or intuitive sequence cutoffs to define sequence clusters of interest (e.g., "main cloud" roughly equivalent to 3% 16S rRNA cutoff, see Rosen et al., 2015), we have used this dynamic approach to let the natural variation, itself, reveal the clusters that most likely fit the dynamic properties of species, given the resolution allowed by the sequence data. In recent pyrosequencing analyses of a gene encoding an essential photosynthesis protein (PsaA), Ecotype Simulation demarcated 18 B'-like, 14 A-like, and 5 A' like Synechococcus PEs (Becraft et al., 2015). Importantly, these analyses demonstrated the ecological distinctness of the PEs: (i) the most abundant PEs had unique environmental distributions, (ii) the PEs showed unique gene expression patterns, (iii) cultivated strains with identical 16S rRNA sequences (see Nowack et al., 2015) from PEs with different depth distributions showed different light adaptations and acclimation responses (Nowack et al., 2015), and (iv) comparative genomic analyses of these strains are beginning to reveal the genetic mechanisms underlying these adaptations (Olsen et al., 2015). Moreover, the membership within PEs was shown to be ecologically interchangeable: (v) individuals predicted to belong to the same PE were distributed similarly along ecological gradients and (vi) the membership within a PE showed coordinated responses to environmental change (Becraft et al., 2015; Olsen et al., 2015).

Recently, Rosen et al. (2015) presented a view of the same Mushroom Spring Synechococcus populations that suggests recombination prevents the formation of such PEs. These authors pyrosequenced numerous PCR-amplified genes to study the genetic diversity within the A and B' lineages of the present model system (that is, the A and B' lineages of Synechococcus in Mushroom Spring) Although, Rosen et al. described their study as a multi-locus sequence analysis, any given organism was sampled for only one locus (i.e., multiple loci were unlinked in these analyses). Rosen et al. (2015) inferred from patterning of single-nucleotide polymorphisms (SNPs) in individual amplicons that recombination and mutation occur at similar rates within the sampled population. Viewing recombination to be "ubiquitous" and, based on comparison of neutral drift and sexual recombination models, they inferred that their "results were inconsistent with the presence of multiple ecotypes" within the A and B' lineages and conjectured that the A and B' lineages of Synechococcus were each behaving "as effectively [a] quasisexual species. . . occupying a broad environmental niche." However, as noted above we previously showed these lineages to each contain multiple, ecologically distinct populations (Becraft et al., 2015; Nowack et al., 2015; Olsen et al., 2015). Clearly, these Synechococcus populations provide a test case for understanding the roles of recombination and ecological specialization in speciation within these lineages. Here we will demonstrate in a true multilocus analysis that each lineage has diversified into ecologically distinct populations in spite of recombination between populations.

Single-locus analyses may suffer from the effects of recombination, such that the evolutionary history of the gene may deviate from that of the organisms (Feil et al., 2001, 2004). Thus, we developed a novel multi-locus sequence analysis (MLSA) approach for cultivation-free sampling and analysis of population structure. MLSA provides a buffer against the effects of recombination in any single gene, increases molecular resolution, and allows detection of recombination events that involve both whole genes and parts of genes because the genes investigated are known to be linked on the same genome.  Furthermore, the development of methods that permit access to large segments of genomes has extended its application to uncultivated organisms (Kashtan et al., 2014). Our cultivation-independent MLSA approach was based on the use of bacterial artificial chromosome (BAC) libraries (Shizuya et al., 1992; Tao et al., 2002; Liles et al., 2008) to sample multiple genes from individual genomes of Synechococcus inhabiting the microbial mat of Mushroom Spring.

BACs were derived from libraries constructed from DNA that was extracted from samples collected at either a 60◦C site or a 65◦C site. Consequently, like the other studies mentioned above, our focus populations have diverged into two major lineages (i.e., the A and B' lineages of Synechococcus, defined by 16S rRNA sequence variants A and B', see above) inferred to rarely recombine between lineages and known to be ecologically distinct (Ward et al., 1998, 2006; Klatt et al., 2011). BACs were screened to identify clones containing a Synechococcus A/B lineage-specific 16S rRNA region. Multiple loci were selected and PCR-amplified from these BAC clones and the sequence data obtained were analyzed using Ecotype Simulation of both concatenated MLSA and single-gene phylogenies, and alternatively, using eBURST analysis, which has routinely been used in MLSA studies to demarcate variants into clusters (see Methods; Feil et al., 2004).

Ecotype Simulation and eBURST use very different approaches to infer populations from sequence variation. Ecotype Simulation is a theory-based algorithm demarcating clades that are consistent with the dynamics of ecotype formation and purging of diversity within ecotypes; eBURST clusters variants on the basis of a pre-determined cutoff—that members of a cluster must be identical at all or all but one of the loci studied; however, concatenated MLSA sequences have been used in phylogenetic analyses (independent from eBURST analyses) to assist in demarcating species clusters (Maiden et al., 1998; Hanage et al., 2005). In previous analyses of this type, suspected recombinant sequences have simply been removed from Ecotype Simulation analyses, on the assumption that these recombinant sequences do not reflect the history of the organism at the locus undergoing recombination and, if not removed, might lead to overestimation of PEs (Koeppel et al., 2008). In this study, however, we address the ability of Ecotype Simulation to predict PEs when recombined sequences are included.

The extent of recombination was assessed in single-gene and concatenated MLSA alignments using programs designed to detect recombination events and phylogenetic incongruency (also assessed by visual comparisons), and by SNP pattern analyses. The ecological distinctness of PEs, predicted either by MLSA or single-locus analysis, was determined through association with the temperature (60◦ or 65◦ ) from which the BACs were obtained. We also tested for finer-scale habitat associations by tracking the distributions of genetic variants of an MLSA locus using pyrosequencing of DNA from samples collected over a range of vertical microenvironments in the 63– 65◦C mat (Becraft et al., 2015). In some cases genomic sequences were available for isolates representative of PEs demarcated in MLSA analyses (Olsen et al., 2015). This permitted comparison of distribution of the MLSA locus with psaA distributions that had previously been determined on the same samples (Becraft et al., 2011).

We will demonstrate not only that recombination has occurred more frequently within than between the A and B' lineages, but that it has neither prevented ecological diversification within each of these lineages nor our ability to detect ecotypes. We suggest that an ecology-based model of speciation is more appropriate than the biological species concept for bacterial speciation in hot spring Synechococcus, and possibly many other bacteria and archaea in diverse environments.

# RESULTS

# BAC Clone Libraries

Metagenomic libraries constructed from the 60◦ and 65◦ mat samples from the Mushroom Spring effluent channel contained 304,128 and 64,512 BAC clones with average insert sizes of 90 and 120 kb, respectively (see Supplemental Data Sheet Section I). BLAST analyses of paired-end sequences of 9216 randomly selected clones from each BAC library showed that the BAC libraries contained the same predominant taxonomic composition as small-insert metagenomes previously produced from the same DNA (Klatt et al., 2011; also see Supplemental Data Sheet Section II, Supplemental Data Tables 3, 4, and Supplemental Data Presentation Figure 2). While the BACs yielded poorer recovery of Synechococcus genes than the small-insert metagenomes, they provided broad coverage of Synechococcus strain JA-3-3Ab [accession: CP000239.1] and JA-2-3B'a(2-13) [accession: CP000240.1] reference genomes, which are representative of A-like and B'-like Synechococcus 16S rRNA lineages (Bhaya et al., 2007; **Figure 1A**, Supplemental Data Presentation Figure 3). BAC clones containing a Synechococcus strain A-like or B'-like 16S rRNA sequence were identified by oligonucleotide probe screening and sequence analysis (Supplemental Data Sheet Section I). Mapping the distribution of the paired-end sequences of each BAC clone relative to reference genomes suggested nearly equal and random coverage of the two unlinked 16S rRNA loci of these genomes (**Figure 1B**, Supplemental Data Presentation Figure

3). Interestingly, many Synechococcus A-like or B'-like BACs containing a 16S rRNA sequence were non-syntenous when compared with the respective genome sequence. Syntenous sequences are considered here to be those whose mate pairs were of the same orientation as the reference genome and were separated by a distance that was similar to the size of the DNA fragments used to construct the metagenomic library. Nonsyntenous sequences are considered here to be those whose mate pairs deviated from the insert length and/or orientation of the reference genome (See Supplemental Data Sheet Section II and reference, Klatt et al., 2011). Such cases were revealed by the paired-end sequences of the same clone matching regions near both 16S rRNA loci, as opposed to one. The percentage of normal- and anti-normal long orientations of paired-end sequences in these non-syntenous clones was greater than in randomly sampled BACs (**Figure 1C**, Supplemental Data Sheet Section II, Supplemental Data Table 5, and Supplemental Data Presentation Figure 3). This suggested that non-syntenous BACs may have been associated with genomic inversions (Rusch et al., 2007).

## MLSA Databases

BACs containing A-like and B'-like Synechococcus 16S rRNA genes were screened by PCR amplification of loci targeted for MLSA analysis (**Tables 1**, **2**; Supplemental Data Tables 1, 2). MLSA loci were selected on the basis of their being (i) singlecopy genes, (ii) under neutral or purifying selection, that (iii) represented a range of degrees of divergences and distances from the 16S rRNA locus and (iv) were not adjacent to potential mobile elements (see Supplemental Data Sheet Sections I, III).

The distributions of MLSA loci among the 267 A-like and 237 B'-like BACs containing the targeted 16S rRNA locus and at least one of the other MLSA loci are shown in Supplemental Data Sheet Section III and Supplemental Data Tables 6, 7. The majority of BACs (94 and 66% of Synechococcus A- and B' like BACs, respectively) exhibited locus combinations consistent with the gene order of the reference genomes. The number of BACs positive for various loci decreased as the distance separating them from the 16S rRNA locus increased, providing additional evidence of genome rearrangements in these lineages (Supplemental Data Sheet Section III and Supplemental Data Presentation Figure 4). There was no obvious pattern of decaying synteny with increasing phylogenetic divergence (see end of legends in **Figure 2** and Supplemental Data Presentation Figure 10 identifying syntenous and non-syntenous clones in the dataset).There were 71 Synechococcus A-like BACs that contained the total set of 7 selected loci (rbsK, PK [locus tag CYA\_2262], hisF, lepB, CHP [locus tag CYA\_2291], aroA and dnaG) in addition to the 16S rRNA/16S-23S rRNA internal transcribed spacer (ITS) region (Supplemental Data Table 6). However, only two Synechococcus B'-like BACs contained 7 of

TABLE 1 | Results from recombination and mutation rate and ratio analyses for A-like Synechococcus BACs.

*CI: 95% credibility interval.*

*a rbsK, CHP, PK, lepB, aroA.*

*bSum of all single gene lengths (number of nucleotides) and analysis pertains to a 73-sequence (145 sequence de-duplicated) 5-locus alignment except for hisF and dnaG which pertain to 7-locus alignment.*

*cWatterson's theta.*

*dMaximum at 4Ner(region); composite likelihood method (LDHat; pairwise) finite sites model (McVean et al., 2002, 2004).*

*eR/*θ *from LDHat is equivalent to the* ρ*/*θ *measure from Clonal Frame and define the ratio of the rates of recombination (R or* ρ*) to mutation (*θ*) illustrating how often recombination occurs relative to mutation (McVean et al., 2002, 2004; Didelot and Falush, 2007; Didelot and Maiden, 2010). For LDHat-pairwise analysis, R/*θ *range includes calculations from compositelikelihood, Rmin and Wakeleys moment method (McVean et al., 2002, 2004).*

*fRatio of the probabilities that a given site is altered through recombination (r) and mutation (m) illustrating how important the effect of recombination was in the diversification of the samples relative to mutation (Didelot and Falush, 2007; Didelot and Maiden, 2010).*

*gThere were no segregating sites in the sequence dataset.*

*hClonal Frame analysis is designed for concatenated multi-locus alignments (Didelot and Falush, 2007; Didelot and Maiden, 2010).*

*iavePWD: average pairwise distance (LDHat pairwise).*

*j varPWD: variance in pairwise distance (LDHat pairwise).*

*kRmin value (LDHat pairwise); minimum number of recombination events describing evidence for recombination in region (McVean et al., 2002, 2004) infinite sites model.*

*lPopulation scaled 4Ner (recombination estimate) by Wakeley Moment method (McVean et al., 2002, 2004).*

*mThe look up table for LDHat pairwise not generated from data but from existing lkgen table provided in LDHat distribution (lk\_n100\_t0.01; https://github.com/auton1/LDhat).*

*CI: 95% credibility interval.*

*aWatterson's theta.*

*bMaximum at 4Ner(region); composite likelihood method (LDHat; pairwise) finite sites model (McVean et al., 2002, 2004).*

*cR/*θ *from LDHat is equivalent to the* ρ*/*θ *measure from Clonal Frame and define the ratio of the rates of recombination (R or* ρ*) to mutation (*θ*) illustrating how often recombination occurs relative to mutation (McVean et al., 2002, 2004; Didelot and Falush, 2007; Didelot and Maiden, 2010). For LDHat-pairwise analysis, R/*θ *range includes calculations from compositelikelihood, Rmin and Wakeleys moment method (McVean et al., 2002, 2004).*

*dRatio of the probabilities that a given site is altered through recombination (r) and mutation (m) illustrating how important the effect of recombination was in the diversification of the samples relative to mutation (Didelot and Falush, 2007; Didelot and Maiden, 2010).*

*eClonal Frame analysis is designed for concatenated multi-locus alignments (Didelot and Falush, 2007; Didelot and Maiden, 2010).*

*f avePWD: average pairwise distance (LDHat pairwise).*

*g varPWD: variance in pairwise distance (LDHat pairwise).*

*hRmin value (LDHat pairwise); minimum number of recombination events describing evidence for recombination in region (McVean et al., 2002, 2004) infinite sites model.*

*iPopulation scaled 4Ner (recombination estimate) by Wakeley Moment method (McVean et al., 2002, 2004).*

*jThe look up table for LDHat pairwise* not *generated from data but from existing lkgen table provided in LDHat distribution (lk\_n100\_t0.01; https://github.com/auton1/LDhat).*

the targeted MSLA loci (Supplemental Data Table 7). In order to obtain a comparable sampling of Synechococcus B'-like and A-like BAC clones, the number of loci for the Synechococcus B'-like population was reduced to three protein-encoding loci (rbsK, aroA and pcrA) plus the 16S rRNA/16S-23S rRNA ITS sequence. The amount of nucleotide divergence at the 16S rRNA locus within A-like and B'-like BACs averaged 0.31% (±0.036% SE) and 0.71% (±0.064% SE), respectively, much less than that between BACS of the different lineages, which averaged 2.73% (±0.024% SE), which was somewhat lower than the 3.6% divergence reported for the Synechococcus strain A and B' genomes (Bhaya et al., 2007). We compared the rbsK sequences of A-like and B'-like BACs to rbsK sequences recovered by direct PCR amplification from the same mat samples and found that the BAC library sequences sampled broadly across B'-like and A-like lineages (Supplemental Data Sheet Section IV and Supplemental Data Presentation Figures 6, 7). All multilocus sequence datasets were assembled into sequence types (STs) comprised of individual BACs with identical sequences at all loci (Supplemental Data Sheet Section V and Supplemental Data Tables 8, 9).

# Recombination and Mutation within and Between A and B' Synechococcus Lineages

Evidence of recombination within the A- and B'-lineages was readily detected by various analyses of both single genes and multiple loci. The ability to detect recombination events was influenced by the number and length of sequences analyzed and by the number of genes included in the analysis. Furthermore, different methods detected different events. For instance, as shown in **Table 3**, RDP4 analysis of 7 loci in 71 A-like BACs detected 11 unique events. Clonal Frame detected 7 of the same events, and manual inspection of SNPs (using the definition from (Feil et al., 2000a), where >1 SNP within any single locus is indicative of recombination) suggested 24 additional unique putative events. Hence, 11–35 unique events were detected. However, in an analysis in which we more than doubled the number of sequences by removing 2 of the 7 genes (see below), RDP4 detected only 5 events; Clonal Frame detected 1 additional event and SNP analysis suggested 2 more events (range 6–8 events). Similarly, in an analysis of 4 genes in 72 B'-like BACs RDP4 detected 4 events, Clonal Frame detected an additional 9 events and SNP inspection suggested 11 additional unique events (range 13–24 events).

The majority of recombination events were detected in the rbsK locus in both lineages; to a lesser extent recombination events were also detected in PK, CHP, and lepB for A-like Synechococcus and in the 16S rRNA/ITS region and pcrA gene for B'-like Synechococcus (Supplemental Data Tables 10–13). Phylogenies forrbsK, constructed using RDP4 with sequence data from either side of the predicted breakpoints, were shown to be incongruent providing additional evidence of recombination at this locus in both lineages (Supplemental Data Presentation Figures 8, 9).

In contrast, evidence of recombination between the A-like and B'-like lineages was rare. Blast analyses of 7 MLSA loci on 71 A-like BACs and 4 MLSA loci on 72 B'-like BACs revealed only three examples, all involving rbsK loci on A-like (n =

TABLE 3 | The number of unique (u) and overlapping (ovl; supported by multiple methods) recombination events recorded by RDP4, Clonal Frame (CF), and single nucleotide polymorphism (SNP) analysis of MLSA datasets.

*All loci were considered in analysis.*

\**Overlapping RDP4 events with CF or SNP analyses were not recorded in CFu or SNPu columns, only the RDP4u column.*

*aSee Supplemental Data Table 10.*

*bSee Supplemental Data Table 11.*

*cSee Supplemental Data Presentation Figure 10.*

*dSee* Figure 2*.*

*eSum of "u" columns only.*

1) or B'-like (n = 2) BACs that were associated with the genome of the other lineage (Supplemental Data Sheet Section VI). Similarly, of the loci used to study both A-like and B'-like Synechococcus lineages, rbsK was the only gene for which RDP4 analysis detected evidence of recombination. Only 3 unique recombination events were identified among 123 combined Alike and B'-like rbsK sequences in the 511 nt overlapping region (all within B'-like BACs; Supplemental Data Table 14).

The importance of recombination relative to mutation for the evolution of A-like and B'-like Synechococcus lineages was estimated independently in two ways: (i) R/θ, which is the ratio at which recombination and mutation occur (McVean et al., 2002, 2004; Didelot and Falush, 2007) using LDHat analysis (reported as a range depending on LDHat method annotated in footnotes of **Tables 1, 2** and equivalent to ρ/θ in Clonal Frame analysis), and (ii) r/m, which is a ratio of the probabilities that a given site was altered through recombination or mutation (i.e., how important the effect of recombination was in the diversification of the sample relative to mutation in terms of the number of resulting nucleotide substitutions) as defined by Clonal Frame analysis (Didelot and Falush, 2007; Didelot and Maiden, 2010). For MLSA concatenated datasets, the ratios (R/θ, ρ/θ, and r/m) differed somewhat in the two lineages. Recombination events were generally less frequent than mutation events (ρ/θ ranging from 0.38 to 0.89 for A-like BACs and 0.24 to 0.54 for B'-like BACs, in Clonal frame analysis; R/θ ranging from 0.4 to 1.0 for A-like BACs and 0.44 to 1.06 for B'-like BACs in LDHat analysis; see **Tables 1, 2**). It was also clear from Clonal frame analysis that recombination had a greater impact than mutation on diversification of the lineage (r/m ratio ranging from 2.60 to 3.40 depending on number of loci for Synechococcus A-like BACs, and 2.42 for Synechococcus B'-like BACs; **Tables 1**, **2**).

# Ecotype Simulation Analyses of Concatenated A-like Synechococcus MLSA Datasets

We focus in the main text on the A-like Synechococcus lineage, which offered a greater number of genes for MLSA analyses (7 loci, hence termed MLSA7). Also, A-like Synechococcus BACs were recovered from both mat samples, thus enabling analysis of habitat association. Comparable results for a 4-locus MLSA analysis (MLSA4) of B'-like lineage variants, which were recovered only from the 60◦C sample, are presented in Supplemental Data Sheet Section VII. For clarity, PEs demarcated by Ecotype Simulation are named according to whether they are based on MLSA data (MLSA7 or MLSA4 followed by the PE designation; e.g., MLSA7 PEA1; MLSA4 PEB'4) or singlelocus data (specified by the gene analyzed followed by the PE designation; e.g., rbsK PEA2).

All individual-locus sequences sampled by BACs and multilocus concatenations constructed from them (STs) were analyzed using Ecotype Simulation (**Figure 2A**, **Table 4**). STs predicted by Ecotype Simulation to be members of the same putative ecotype in the MLSA7 analysis are given the same color, but members of different PEs are colored differently, and colors are maintained throughout all MLSA7 Figures. The MLSA7 concatenated phylogeny and the single-locus phylogeny for rbsK sequences sampled by the same BACs are compared in **Figures 2A,B** respectively, which also show PEs predicted by Ecotype Simulation. Phylogenetic comparisons with other loci (CHP, lepB, dnaG PK, aroA, and hisF) are discussed in Supplemental Data Sheet Section VIII and are shown in Supplemental Data Presentation Figure 14. Ecotype Simulation predicted 2–7 A-like PEs from individual BAC loci, depending on the locus analyzed. Among individual loci, the greatest number of PEs was predicted from the rbsK gene (**Figure 2B**), which had the greatest average evolutionary divergence among the loci studied (**Table 4**). Ecotype Simulation predicted 13 A-like MLSA7 PEs from concatenated MLSA sequence data (**Figure 2A**), though five of these PEs were based on a single ST sequence (i.e., singletons) and one was based on two occurrences of the same ST (i.e., a doubleton). MLSA7 PE clades contained from 1 to 28 BACs. Three of the 13 MLSA7 PE clades for Synechococcus A-like BACs contained a dominant variant (i.e., 2–14 clones that were identical at all MLSA7 loci) and singleton variants. In two cases, MLSA7 PEs A4 and A7, subdominant variants were also observed, along with singletons. Supplemental Data Table 23 shows how STs, dominant variants and subdominant variants were distributed among MLSA7 PEs.

Only one case of complete correspondence between rbsK PEs and MLSA7 PEs was observed (i.e., MLSA7 PE A9 = rbsK

shaded in gray and STs that contained a combination of sequences that were syntenous and non-syntenous are indicated by an asterisk; syntenous STs are not shaded and are not annotated by an asterisk. Shared SNP pattern between STs of two distinct PEs indicated by a bidirectional arrow colored according to gene (yellow = *lepB*). Bootstrap values are provided for major nodes. Reference genome indicated by "SynA"; Genbank accession number CP000239.

PEA 4) (green colored STs in **Figure 2**). Otherwise, the rbsK and MLSA7 phylogenies were mainly incongruent, providing further evidence that recombination complicates single-locus analyses of PEs (compare **Figures 2A,B**). For instance, some members of MLSA7 PE A4 were classified into three rbsK PEs (A1, A2, and A7; red colored sequences in **Figure 2**). Ecotype Simulation analysis of concatenated MLSA7 loci mainly split individual PEs demarcated by rbsK analysis into multiple MLSA7 PEs, as shown by dashed and solid black arrows in **Figure 2** (e.g., rbsK PE2 was split into MLSA PEs 3 (pink), 4 (red), and 5 (blue).

# Evidence of Habitat Association

In order to obtain a greater number of A-like BACs from samples collected at both temperatures sufficient to observe

TABLE 4 | Ecotype Simulation and eBURST output for 71 A-like Synechococcus BACs that contained all 7 loci assayed for the study.

*aMolecular resolution of hisF was not significant enough for Ecotype Simulation to effectively predict ecotypes.*

*bCalculated from concatenated sequence datasets only.*

*cNot enough 65*◦ *sequences to determine sample specificity of demarcated PEs.*

*dClonal complexes only apply to concatenated sequence datasets.*

habitat associations, we performed two separate MLSA analyses of 5 of the 7 loci (rbsK, lepB, aroA, CHP, and PK; MLSA5, see Supplemental Data Sheet Section VIII). One included 49 concatenated sequences with nearly equal representation of sequences from 60 and 65◦C (Supplemental Data Presentation Figure 15). The other included the full dataset of 145 concatenated sequences that could be compared at these 5 loci (Supplemental Data Presentation Figure 16). Fisher's exact tests of contingency on both sets of BACs from 60 and 65◦C yielded significant evidence of heterogeneity in habitat associations (i.e., collection temperatures; p-value = 0.001). Additionally, several MLSA5 PEs were comprised of multiple BACs obtained only from one of the two temperatures.

In order to investigate vertical distributions of alleles we pyrosequenced the rbsK locus in DNAs extracted from 80µmthick cryotome sections from different depths in the top green layer of the 63–65◦C mat, where A-like PEs predominate. Pyrosequencing analysis of the rbsK locus resulted in a phylogeny (Supplemental Data Sheet Section IX; Supplemental Data Presentation Figure 19) that closely resembled the rbsK phylogeny based on sequences sampled by BAC clones (compare **Figure 2B** and Supplemental Data Presentation Figure 19). In cases where pyrosequences matched only STs of the same PE or subclade, we could test the hypothesis that these PEs represent ecologically distinct populations. As shown in **Figure 3A**, the subclade representing MLSAs PE A4 and A5 (red and blue STs in **Figure 2A**), which were grouped together into rbsK PE A2 (**Figure 2B**), were shown to exhibit different vertical distributions (Fisher's exact test, ANOVA and G-tests p < 0.001). rbsK sequences representative of MLSA PE A4 declined with depth, whereas those representative of MLSA PE A5 increased with depth. This provides an example of how a single-locus analysis can lump populations that are ecologically distinct, giving the impression of a clade containing ecologically heterogeneous members. Additionally, the availability of genome sequences of isolates representative of MLSA7 PEs A5 and A9 (Olsen et al., 2015; and see Supplemental Data Sheet Section X) allowed us to associate rbsK and psaA sequences corresponding to the same PEs. In these cases we could demonstrate correspondence between rbsK distributions and the psaA distributions observed in single-locus analysis of the same samples by Becraft et al. (2015) (see also **Figures 3B,C**, where MLSA7 PE 5 = psaA PE A6 and MLSA7 PE 9 = psaA PE A14). Both of these PEs have been shown to be most prevalent toward the bottom portion of the upper green layer of the mat at this temperature (Becraft et al., 2015), and representative isolates have been shown to possess genes that enable harvesting of lower amounts of light, characteristic of in situ light conditions (Becraft et al., 2015; Nowack et al., 2015; Olsen et al., 2015; also see Supplemental Data Sheet Section X).

# eBURST Analyses

The eBURST algorithm assembles STs into clonal complexes (Maiden et al., 1998; Feil et al., 1999, 2000a, 2004) that are typically composed of a single predominant genotype (consensus group or here, dominant variant), plus variants that are identical to the predominant genotype at all but one locus (single-locus variants; Feil and Spratt, 2001; Feil et al., 2004). Single-locus variants can differ from the dominant variant by any number of SNPs as long as those SNPs occur in only one locus (Feil et al., 1999, 2000b, 2003). Only 3 clonal complexes were observed with eBURST analyses of the MLSA7 dataset (**Figure 4**, **Table 5**), and the correspondence between eBURST's clonal complexes and Ecotype Simulation's PEs was imprecise (see **Figure 4** and Supplemental Data Sheet Section VIII for details of these analyses and results). In the eBURST analyses of the MLSA5 datasets, Fisher's exact tests did not indicate differences among clonal complexes in their habitat associations (Supplemental Data Presentation Figure 17, p-value 0.943; Supplemental Data Presentation Figure 18, p-value 0.734).

to ST10.

#### TABLE 5 | eBURST analysis of clonal complexes for the 7-locus MLSA of A-like Synechococcus BACs.

*The number of BACs within an ST is in parentheses when greater than 1. Superscripts next to the allele number denote number of nucleotide differences compared to the consensus sequence. Corresponding to* Figure 4*.*

*DV, dominant variant.*

# Evidence of Historic and Recent Recombination Events

The presence of historic recombination events was further investigated after it was initially found through SNP analysis that identical SNP differences were observed in two different MLSA7 PEs. **Figure 5** shows SNP differences between variants that were grouped into MLSA7 PEs and clonal complexes sharing the same dominant variant. Most differences resulted either from the exclusion of STs containing SNPs in multiple loci from clonal complexes by eBURST, or from the exclusion of STs containing a large number of SNPs from phylogenetic clusters demarcated as PEs by Ecotype Simulation. Of particular interest, however, was the observation of identical SNP differences in the CHP locus in two different MLSA7 PE A4 STs relative to two different dominant variants of MLSA7 PE A7 (ST7 vs. ST6 and ST11 vs. ST1). This prompted us to consider the possibility that these SNP patterns were clade-specific. We investigated recombination events occurring on internal nodes (suggestive of historic recombination) using Clonal Frame. Clonal Frame predicted several internal nodes in both the Synechococcus A-like and B'-like lineages where there was a high probability of recombination having occurred in one or more loci (closed colored stars in **Figure 2A** and Supplemental Data Presentation Figure 10). This revealed a historical record of the impact recombination has had on the diversification of these organisms. Detailed manual inspection of SNP patterns in all clades allowed us to suggest additional putative recombination events in the phylogeny of these Synechococcus lineages (open colored stars in **Figure 2A** and Supplemental Data Presentation Figure 10; also see Supplemental Data Sheet Section XI and Supplemental Data Presentation Figure 21). For instance, all variants grouped into the clade containing MLSA7 PEs A2- A4 appear to have vertically inherited a recombinant CHP gene after an early recombination event (closed green star) and variants within MLSA7 PEs A2, A3, and a portion of PEA4, appear to have vertically inherited an rbsK sequence from a

subsequent recombination event (closed blue star). A third putative recombination event involving the rbsK locus (open blue star) may distinguish MLSA7 PEs A2 and A3 from MLSA7 PE A4.

These SNP analyses also revealed several examples of singleton MLSA7 PEs associated with recombination patterns involving a large number of SNPs [MLSA7 PEs A1 (ST4), A2 (ST27), A8 (ST12), and A11 (ST25); see **Figure 2A**]. These singleton-PE variants formed long-branches and were closely related to sister PEs comprised of multiple variants, but different in many SNPs from the dominant variant of that PE. For example, in **Figures 2A**, **5A**, singleton-PE A8 (ST12) was closely related to PE A7, but differed by 19 SNPs in the dnaG locus, suggesting that a recent recombination event might have caused enough sequence difference to prevent its accurate demarcation by Ecotype Simulation. There was only one case in which an identical SNP pattern shared by two STs (ST25 in singleton MLSA7 PE A11 and ST 34 in MLSA7 PE A9; lepB) was not clade-specific (**Figure 2** and Supplemental Data Presentation Figure 21 panel I).

Similar analysis was performed on the Synechococcus B'-like BACs with the interesting result that 2–3 recombination events were associated with some divergences, even though only 4 loci were included in the B'-like MLSA4 phylogeny (Supplemental Data Sheet Section VIII and Supplemental Data Presentation Figures S10, S12, S13).

# DISCUSSION

# Recombination within and Between Lineages of Synechococcus

There is abundant evidence that genomic rearrangements and recombination have influenced the evolution of the Synechococcus populations inhabiting the Mushroom Spring microbial mat. As in previous investigations of Sulfolobus islandicus (Cadillo-Quiroz et al., 2012), Wolbachia (Ellegaard et al., 2013), and Vibrio cyclitrophicus (Shapiro et al., 2012) mentioned above, recombination between the divergent A and B' lineages of Synechococcus appears to be less frequent than recombination within each lineage. Specifically, only approximately 2% of BACs exhibited evidence of recombination between A- and B'-like Synechococcus BACs (see Results, Supplemental Data Sheet Section VI and Supplemental Data Table 14). This corroborates our previous metagenomic analyses, in which we found a low percentage of metagenomic clones with end sequences stemming from both the A-like and B'-like lineages (Klatt et al., 2011) (Supplemental Data Sheet Section VI). In contrast, within-lineage recombination was greater, with 5.5% of A-like and 33.3% of B'-like Synechococcus BACs exhibiting evidence of recombination for approximately equal numbers of loci analyzed (MLSA5, 145 sequences and MLSA4, 72 sequences for Synechococcus A- and B'-like BACs respectively; **Table 3**).

Recombination events have occurred at various times in the history of the A-like and B'-like clades: some recombinant segments are fixed in large subclades (i.e., evidenced by clustering of SNPs in MLSA PEs 2, 3, and 4 in Supplemental Data Presentation Figures 21A,D,G or MLSA PEs 10, 11, 12, and 13 in Supplemental Data Presentation Figures 21C,F,I; for analogous B'-like fixed recombinant segments see Supplemental Data Presentation Figure 13), indicating a recombination event early in the clade's history followed by vertical inheritance of the recombined gene, while others are specific only to singleton (one doubleton) variants demarcated as PEs based on recent recombination in single loci (i.e., MLSA PEs 1, 2, 8, 10, 11, and 13 in **Figure 2A**). These results suggest that recent recombination events might cause overestimation of PEs by Ecotype Simulation. That is, a single recombination event, if involving a long enough fragment, can cause Ecotype Simulation to split a PE. In the case of eBURST, a single recombination event of any size cannot split a consensus group; however, two or more recombinants of any size will separate a recipient organism into a different consensus group.

The rates of recombination within the A and B' lineages were estimated to be less than or about the same as that of mutation (**Tables 1**, **2**). However, we determined that recombination had a somewhat greater impact than mutation on sequence change within these populations (**Tables 1**, **2**; Clonal Frame, r/m range 1.58–4.92). The impact of recombination (r/m) was estimated at 2.6 and 2.4, for the Synechococcus A-like and B'-like lineages, respectively. These r/m ratios are similar to values reported for other cyanobacteria (Vos and Didelot, 2009) and are in the same range as those reported for S. islandicus, Wolbachia, and Halorubrum (Papke et al., 2004; Whitaker et al., 2005; Ellegaard et al., 2013). The recent study on Synechococcus by Rosen et al. (2015), did not provide r/m, but their estimate of the ratio of recombination to mutation rates (ρ/θ) was close to 1, which fell within the R/θ range provided by our LDHat analysis (**Tables 1**, **2**). Why do such different interpretations arise from such similar results? Although, Rosen et al. interpreted from a very similar recombination/mutation ratio that recombination is "ubiquitous," we note that this ratio is not extraordinarily high when compared against the many groups of bacteria whose recombination rates were quantified by Vos and Didelot (2009). Furthermore, Rosen et al. have not considered the theoretical work of Haldane (1932), which shows that such low rates of recombination cannot prevent ecological diversification.

Given the rates of mutation estimated for various bacteria (Drake, 2009; Vos and Didelot, 2009), we may conclude that a typical gene in the Synechococcus populations is experiencing recombination at an extremely low absolute rate, so that recombination between nascent species is unlikely to hinder their divergence (Haldane, 1932; Vos, 2011; Wiedenbeck and Cohan, 2011). Low-frequency recombination may allow shuffling of niche-neutral alleles, which do not determine ecological niche, but low-frequency recombination of niche-determining alleles from other populations can be countered by natural selection (Haldane, 1932). Maladaptive niche-determining alleles from other populations will be kept at negligible levels, maintaining the integrity of populations' adaptations (Wiedenbeck and Cohan, 2011). Rates of recombination observed in Synechococcus should not lead us to dismiss a priori the possibility of ecological diversification within either clade.

We therefore disagree with Rosen et al.'s conclusion (Rosen et al., 2015) that frequent recombination has caused each of the A and B' clades of Synechococcus to be a "quasisexual population," where close relatives are prevented from diversifying into longstanding, ecologically homogeneous populations, or ecotypes. While Rosen et al. (2015) have shown a shuffling of alleles (really gene segments) within the A and B' clades, they have not shown that the recombination has randomized associations between niche-determining alleles, which are responsible for ecotype divergence. Moreover, we do not see how they could rule out ecological diversification without studying ecology. As we have shown in many previous papers, as well as in the present study, there is ample evidence for the coexistence of longstanding, ecologically divergent populations within both the A and B' clades, as we discuss below. Rosen et al. mentioned some of these studies, but did not use an ecological model in their work, as we have done here.

As with most (if not all) bacteria, recombination has been sufficient in Synechococcus to create new ecological diversity. Within the A lineage, for example, low-light PEs differ from high-light PEs in having acquired a gene cluster likely coding for an additional photoreceptor (Olsen et al., 2015). It does not take a recombination rate much higher than that of mutation to transfer a gene that is adaptive in a new niche. Only one adaptive transfer of a gene (or one mutation event); see Bantinaki et al. (2007) is required to found a new ecotype and thereby bring about ecological diversification. On the other hand, a much higher rate of recombination is required to prevent adaptive divergence between nascent populations; this requires recurrent inter-population transfer of genes poorly adapted to their new niche (Wiedenbeck and Cohan, 2011).

# Ecological Diversification within the A-like and B'-like Synechococcus Lineages

The Ecotype Simulation analysis of concatenated MLSA sequences hypothesized 13 and 29 PEs within the A-like and B'-like clades, respectively. This was greater than the number of PEs predicted from single-locus analyses of the most divergent locus, rbsK. This is likely because of greater resolution provided by MLSA7 (as well as MLSA4; see Supplemental Data Presentation Figure 10), as evidenced by the fact that MLSA7 split numerous PEs that were lumped by rbsK analysis. We were able to demonstrate the ecological distinctness of some of the A-like putative ecotypes, first with respect to the two habitat types investigated in this study (60 and 65◦C). That is, some hypothesized A-like ecotypes were sampled entirely from different habitats, while others were quantitatively different in their habitat associations. Also, pyrosequencing of PCR-amplified rbsK genes typical of different MLSA7 PEs demonstrated fine-scale distribution differences of MLSA7 PEs along vertical gradients. Additionally, the availability of genomes of isolates representative of MLSA7 PEs permitted us to demonstrate correspondence of rbsK distributions with distributions observed earlier in psaA analyses (Becraft et al., 2015). This led to demonstration of adaptation consistent with distribution patterns (Nowack et al., 2015) and genomic evidence of mechanisms likely to underlie these adaptations (Olsen et al., 2015). For instance, MLSA PEs 5 and 9 correspond to psaA PEs A6 and A14, whose distribution deep in the upper green layer of the mat (Becraft et al., 2015) is consistent with the presence of auxiliary photosystem genes that are likely to be associated with the ability to acclimate to low irradiance (Nowack et al., 2015; Olsen et al., 2015); and see Supplemental Data Sheet Section X). Becraft et al. (2011, 2015) also describe the existence and ecological distinction of B'-like Synechococcus ecotypes, as well as the ecological interchangeability of the members of each ecotype.

We conclude that both the A-like and B'-like Synechococcus lineages have clearly diversified into many younger species populations, whose niches correspond in part to differences along temperature and vertical gradients (Becraft et al., 2015). Our approach did not permit us to evaluate the role of positive selection within our dataset as our genes were under neutral/purifying selection. Evaluation of full genomes of strains of closely related species will allow us to further examine the roles that recombination, mutation, selection and neutral processes have played in the diversification of Synechococcus.

# Whither the Biological Species Concept?

That recombination has occurred more frequently within closely related clades than between them has previously been taken as support for the Biological Species Concept in bacteria and archaea (Cadillo-Quiroz et al., 2012; Shapiro et al., 2012; Ellegaard et al., 2013). However, the quintessential idea of the Biological Species Concept is that sexual isolation is required for irreversible divergence and ecological distinctness in the origin of species (Cohan, 2013). Support for the Biological Species Concept requires more than a correlation between sexual isolation and divergence among clades; one must demonstrate that reduced recombination was necessary for their divergence.

We argue that such evidence is lacking in the above-cited studies on Vibrio, Wolbachia, and Sulfolobus. As seen in the present work, Synechococcus lineages have been noted to have recombination levels similar to those reported for Halorubrum and Sulfolobus. Like these organisms, major Synechococcus clades that are highly divergent (3.6% diverged at the 16S rRNA locus and 18% diverged in average nucleotide identity (ANI), Bhaya et al., 2007) are ecologically distinct. However, our highresolution analysis of closest relatives within each clade has shown ecological speciation among organisms that have identical or nearly identical 16S rRNA sequences and have diverged only on the order of 1% ANI (Olsen et al., 2015; also see Supplemental Data Table 25) and show minimal sexual isolation. Similarly, in work on another system we have found ecotype differentiation among Bacillus clades that show very little sexual isolation (Majewski and Cohan, 1999; Connor et al., 2010). We hypothesize that analysis of fine-scale ecological divergence within Halorubrum or Sulfolobus may yield similar results (Papke et al., 2007; Cadillo-Quiroz et al., 2012). Recombination does not appear to have impeded ecological diversification among close relatives that are not sexually isolated. Thus, sexual isolation, in this context, is unnecessary for the origin of species (Wiedenbeck and Cohan, 2011).

Instead, sexual isolation is more likely to be the result of divergence among closely related groups than the cause. Not only do ecologically distinct populations tend to diverge over time into separate sequence clusters (Cohan and Perry, 2007), but ecological divergence may also cause occupancy of different microhabitats, reducing opportunities for genetic exchange (Shapiro et al., 2012); in addition, functional divergence may result in reduced fitness of recombinants (Retchless and Lawrence, 2010). For other mechanisms of sexual isolation, also see Wiedenbeck and Cohan (2011). Regarding evidence for the Biological Species Concept provided by previous bacterial and archaeal studies, we note that two of these studies have provided no direct evidence that the groups of sequences studied have distinct ecological properties (Cadillo-Quiroz et al., 2012; Ellegaard et al., 2013), and rely primarily on the coexistence of sequence clusters as evidence that they utilize different resources. A third study by Shapiro et al looking at genetic exchange between two recently diverged populations of Vibrio (Shapiro et al., 2012) reported similar patterns of sub-clustering within the Vibrio populations (their Figure S2), and we predict that these sub-clusters would show ecological diversification if the authors were to test this. We predict each subgroup would represent an ecotype, subject to genomic sweeps. The fourth study (Rosen et al., 2015) claimed that there cannot be discrete ecotypes within the Synechococcus A and B' clades, owing to recombination, yet the existence of such discrete ecotypes has been amply demonstrated here and in many earlier papers (Becraft et al., 2015; Nowack et al., 2015; Olsen et al., 2015).

Most importantly, an argument for the Biological Species Concept should provide evidence that ecological diversification was prevented among groups hypothesized to recombine more frequently. This would require a search for ecological diversification among the close relatives recombining at the highest rates. While previous studies have not attempted to search for ecological diversity among closest relatives (within the major clades); this has been a principal aim of our study.

# CONCLUSION

In conclusion, it appears that we have found some of the most newly divergent ecotypes of hot spring Synechococcus. They do recombine, and recombination has not prevented their divergence, as had been suggested by Doolittle and Zhaxybayeva (2009). Because a reduction in recombination does not appear necessary for ecological divergence, and because ecological patterning of bacterial and archaeal diversity is common, the Biological Species Concept does not appear appropriate for bacterial and archaeal speciation. While recombination has not hindered speciation, ongoing genomic and metagenomic analyses have begun to reveal evidence suggesting that (i) horizontal genetic transfer and recombination have likely led to ecological diversification by introducing niche-transcending adaptations to new lineages (Gogarten et al., 2002; Lawrence, 2002; Doolittle and Zhaxybayeva, 2009; Wiedenbeck and Cohan, 2011; Olsen et al., 2015), and (ii) mutation and selection have also played a role in the evolutionary ecology of Synechococcus ecological species (Olsen et al., 2015). Sub-clustering within major clades predicted to be ecologically distinct has been observed across a variety of taxa and habits (Smith et al., 2006; Zinser et al., 2006; Hunt et al., 2008; Connor et al., 2010; Denef et al., 2010; Mazard et al., 2012), suggesting that our observations may have general importance to the issue of species and speciation in bacteria and archaea. It is important to keep in mind that the molecular dimension of microbial species may be smaller than we have previously thought, and that it is important to understand theories underlying species and speciation.

# MATERIALS AND METHODS

Samples were collected on 2 October 2003 from the effluent channel of Mushroom Spring (44.5386◦N, 110.7979◦W), an alkaline siliceous hot spring in the Lower Geyser Basin of Yellowstone National Park, WY, at two temperatures, 60 and 65◦C. DNA extraction and construction of bacterial artificial chromosome (BAC) clones and oligonucleotide probe screening to identify BACs containing a Synechococcus A/B lineage-specific 16S rRNA gene are described in the Supplemental Data Sheet Section I. DNA from parallel samples was used for analyses of the 16S rRNA and ITS region on the BACs (Ward et al., 2006) and PCR-amplified single protein-encoding loci (rbsK, aroA, and apcAB) (Melendrez et al., 2011), as well as for construction of small-insert (2–12 kb) metagenomic libraries (Bhaya et al., 2007; Klatt et al., 2011).

# Locus Selection

From one hundred genes that were upstream and downstream of the two unlinked 16S rRNA genes in the Synechococcus strain A and B' genomes, loci were selected based on the following characteristics: (i) presence in both genomes, (ii) range of distances from the 16S rRNA locus, (iii) high degree of nucleotide divergence between Synechococcus strain A and B' homologs, (iv) high average nucleotide divergence and variance of metagenomic homologs from Synechococcus strain A and B' homologs, (v) not under positive evolutionary selection (dN/dS values < 1) and (vi) functionally useful for gene expression studies. In addition, loci adjacent to transposons or mobile elements were avoided due to the greater possibility of co-migration of adjacent loci. Loci used in this study are presented in Supplemental Data Tables 1, 2.

# PCR Amplification and Sequencing of BAC MLSA Loci

Primers for the separate amplification of Synechococcus A-like and B'-like BAC target genes (Supplemental Data Tables 1, 2) were designed and tested as described for amplification of rbsK, aroA, and apcAB gene segments by Melendrez et al. (2011), which also provides details of PCR conditions and gel-based size verification and purification of amplicons. Cycling conditions for the dnaG, pcrA, ispE, sufB, and argD genes were the same used for aroA and rbsK genes, and cycling conditions for the protein kinase (PK), lepB, and accC genes were the same as for the apcAB gene, except for the use of 40 cycles instead of 30 cycles. Cycling conditions for hisF and the conserved hypothetical protein (CHP [locus tag CYA\_2291]) gene were: an initial denaturing step at 94◦C (2 min) followed by 20 cycles of 94◦C (1 min), 60◦C (1 min), and 72◦C (1 min); then 20 cycles of 94◦C (1 min), 55◦C (1 min), and 72◦C (1 min) with a final extension at 72◦C for 10 min and storage at 4◦C. BAC clones exhibiting no PCR product for any gene were amplified two more times (directly from the BAC clone sample) to increase confidence that the gene was absent.

# Sequencing

Purified PCR products were sequenced using the forward and reverse primers for each gene described in Supplemental Data Tables 1, 2 and the BigDye v.3.1 cycle sequencing kit (Applied Biosystems) at the University of Nevada-Reno Sequence Center (Reno, NV). The sequences have been submitted to GenBank, via the BankIt submission tool; Accession numbers: pcrA; KT425377-KT425447, rbsK; HQ662694- HQ662843, HQ187926-187996, KT425519-KT425589, and KT426096-KT426240, 16S rRNA; KT425590-KT425660, PK (protein kinase); KT425661-KT425805, CHP (conserved hypothetical protein); KT425806-KT425950, lepB; KT425951- KT426095, aroA; HQ187854-HQ187925, KT425447-KT425518, and KT426241-426385, dnaG; KT426386-KT426455, hisF; KT426456-KT426525.

# Sequence Alignment and Phylogenetic Analysis

Bidirectional sequence data for each locus were analyzed using Sequencher v.4.8. Sequence data were then analyzed with NCBI-BLAST nr/nt (Altschul et al., 1990) to determine whether the top match was with genomic homologs of Synechococcus strains JA-3-3Ab [accession: CP000239.1], JA-2-3B'a(2-13) [accession: CP000240.1], or neither. Sequences used in MLSA were concatenated using a custom perl script (available from J.M. Wood) and uploaded in MEGA4. Alignments were made using the ClustalW algorithm in the MEGA4 software (Tamura et al., 2007) or the MUSCLE (Edgar, 2004) algorithm implemented in Seaview (Gouy et al., 2010). The Synechococcus strain A and B' genomes were included as references (Bhaya et al., 2007). MLSA including newly obtained isolates are described in Supplemental Data Sheet Section X. Alignments were analyzed using jModeltest2 (Darriba et al., 2012) to determine the best model fit for maximum likelihood tree construction. Maximum likelihood trees were constructed using PhyML (Guindon et al., 2010) as implemented in Seaview (Gouy et al., 2010) with aLRT support and edited using MEGA4 (Tamura et al., 2007) or Figtree (Rambaut, 2009). Estimates of average evolutionary divergence for single and multi-locus phylogenies were computed over all sequence pairs using the Maximum Composite Likelihood method in MEGA4 as previously described (Tamura et al., 2004, 2007; Melendrez et al., 2011).

# Ecotype Simulation and Demarcation

Concatenated and single-locus sequence alignments were analyzed using Ecotype Simulation to predict the number of PEs (n), rates of periodic selection (sigma), ecotype formation (omega), and 95% confidence intervals (CI) for all parameters at the best precision match between observed and simulated data for that sequence dataset (between 1.25x and 2x). Ecotypes were manually demarcated conservatively as previously described (Cohan and Perry, 2007; Koeppel et al., 2008; Becraft et al., 2011; Melendrez et al., 2011) (http://fcohan.web.wesleyan.edu/ ecosim/). Statistical tests were conducted using the R statistical package (http://cran.r-project.org/). Groups (rows) were defined as PEs that contained >5 sequences. The variables being tested were temperature (60 and 65◦C) and vertical (80µm intervals) distributions. Significance level was set at 0.05.

# Multi-Locus and eBURST Analyses

Alignments for each gene used in MLSA were organized into groups that were 100% identical at the nucleotide level using Sequencher v 4.8. Allele types were assigned for each unique sequence and were used to generate allelic profiles; BACs with identical allelic profiles at all loci were assigned as unique sequence types (STs), as described in Supplemental Data Tables 8, 9, 19, 20. Allelic profiles and their unique ST designation were uploaded into eBURST (Feil et al., 2004; Spratt et al., 2004) and population snapshots were generated to view A-like and B'-like Synechococcus diversity. Clonal complexes (CC) were defined as a consensus group of BACs and at least 3 single-locus variants, as suggested (Feil et al., 2004). Population snapshots to visualize CCs with less stringent criteria were generated by defining CCs as a consensus group with at least 2 single-locus variants and/or allowing for double-locus variants to be included in the CC and are defined as sub-clonal complexes (see Supplemental Data Sheet Section XII).

# Single Nucleotide Polymorphism Analysis

SNPs were analyzed by comparing sequences of PE clade variants and single-locus variants with that of the dominant variant of the same PE clade or consensus group of the same CC (often the same as the dominant variant) using the Perl program, Pigeon (https://github.com/sandain/pigeon). Pigeon reads in the FASTA file and compares each nucleotide in the dominant variant to the corresponding nucleotide in each PE variant (or single-locus variant), locating and reporting the position of SNPs.

# Detection of Recombination Signals

Five methods were utilized to detect recombination in singlegene and concatenated datasets; (i) outlier detection on phylogeny followed by BLAST (http://blast.ncbi.nlm.nih.gov/ Blast.cgi), (ii) RDP version 4 (RDP4) (Martin et al., 2015), (iii) Clonal Frame version 1.1 (Didelot and Falush, 2007), (iv) SNP analysis (PIGEON; https://github.com/sandain/pigeon), and (v) LDHat v2.2 (McVean et al., 2002).

In RDP, the window size was set to 30 as recommended (RDP4 Instruction Manual, available at http://web. cbio.uct.ac.za/~darren/RDP4Manual.pdf). In MaxChi and Chimera the 'variable window size' was used. In Siscan the window size was set to 200 bp with a step size of 20.

Recombination signals were considered present if they could be detected by at least 3 methods within the RDP4 package at significant p-values as suggested (RDP4 manual).

(iii) Clonal Frame identifies evidence of recombination in multilocus sequence datasets (Didelot and Falush, 2007). Extended multi-FASTA formatted alignments (XMFA) were generated using the MAUVE program (Darling et al., 2004, 2007) and given file extensions of.dat. DAT files were loaded into Clonal Frame and analyzed using default settings with 1,000,000 iterations after burn in, 250,000 burnin, sampling every 100 iterations and verbose mode (-v,-x, y,-z options). OUT files with phylogenies and statistics (mutation events, recombination events, rho/theta and r/m) were visualized using the cfgui.bat program included in the Clonal Frame installation. Recombination events identified along branches are noted with closed stars for MLSA datasets in **Figure 2A** and Supplemental Data Presentation Figure 10.

# Pyrosequencing

PCR primers and protocols described by Melendrez et al. (2011) or in Supplemental Data Table 1 were used to amplify rbsK sequences from Mushroom Spring samples collected at ∼63– 65◦C and sectioned along vertical gradients (∼80µm intervals using a cryotome). Samples were collected on 13 and 14 September 2008 and DNA was extracted as described in Becraft et al. (2011). Barcoding and Ti454-sequencing were completed at the J. Craig Venter Institute according to the GS FLX Titanium Series Rapid Library Preparation Method. DNA was

# REFERENCES

Allewalt, J. P., Bateson, M. M., Revsbech, N. P., Slack, K., and Ward, D. M. (2006). Effect of temperature and light on growth of and photosynthesis by Synechococcus isolates typical of those predominating in the octopus sheared using the Covaris S2 System, and qPCR was used to accurately estimate the number of molecules needed for emPCR. BAC clone MLSA sequences were trimmed at the 5′ and 3′ ends to match the average sequence length, which collapsed some SNP variation contained in the rbsK single-locus analysis and caused a decrease in molecular resolution. Homopolymer inserts were removed by aligning sequences to the rbsK reading frame of genomes from the A or B' Synechococcus isolates. High- frequency sequences (sequences totaling ≥50 identical copies in all combined samples) were then identified and PEs were demarcated from high-frequency sequences using Ecotype Simulation. Population percentages based on high-frequency sequences that matched STs of a single MLSA PE were calculated as described in Becraft et al. (2015). Sequence data for rbsK was deposited into the MG-RAST version 4 database (https:// metagenomics.anl.gov/) under accession numbers 4642340.3- 4642367.3.

# FUNDING

This research was supported by the National Science Foundation Frontiers in Integrative Biology Research Program (EF-0328698), the National Aeronautics and Space Administration Exobiology Program (NAG5-8807 and NX09AM87G), and the U.S. Department of Energy (DOE), Office of Biological and Environmental Research (BER), as part of BER's Genomic Science Program (GSP). This contribution originates from the GSP Foundational Scientific Focus Area (FSFA) at the Pacific Northwest National Laboratory (PNNL), contract no. 112443. It was also supported by Montana Agricultural Experiment Station project 911352.

# ACKNOWLEDGMENTS

We appreciate the assistance from National Park Service personnel at Yellowstone National Park. Special thanks to Darren Martin from the University of Cape Town, South Africa for assistance with the RDP4 software suite, Dr. Timothy McDermott from Montana State University, USA for use of equipment required for radioactive probing of BAC clones and the laboratory of Dr. Al Jesaitis from Montana State University, USA for use of the phosphoimager to identify BACs containing 16S rRNA loci.

# SUPPLEMENTAL DATA

The Supplemental Data for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fmicb. 2015.01540

spring microbial mat community of Yellowstone National Park. Appl. Environ. Microb. 72, 544–550. doi: 10.1128/AEM.72.1.544-550.2006

Altschul, S. F., Gish, W., Miller, W., Myers, E. W., and Lipman, D. J. (1990). Basic local alignment search tool. J. Mol. Biol. 215, 403–410. doi: 10.1016/S0022- 2836(05)80360-2

predominant populations and their functional potential. ISME J. 5, 1262–1278. doi: 10.1038/ismej.2011.73

lusitaniae revealed by multilocus sequence typing. PLoS ONE 3:e4002. doi: 10.1371/journal.pone.0004002

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

Copyright © 2016 Melendrez, Becraft, Wood, Olsen, Bryant, Heidelberg, Rusch, Cohan and Ward. 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.

# Characterization of a single mutation in TraQ in a strain of Escherichia coli partially resistant to Q**β** infection

### *Akiko Kashiwagi\*, Hikari Kitamura and Fumie Sano Tsushima*

Faculty of Agriculture and Life Science, Hirosaki University, Hirosaki, Japan

#### *Edited by:*

Hiroyuki Futamata, Shizuoka University, Japan

#### *Reviewed by:*

Evelien M. Adriaenssens, University of Pretoria, South Africa Hailan Piao, Washington State University in Tri-Cities, USA

#### *\*Correspondence:*

Akiko Kashiwagi, Faculty of Agriculture and Life Science, Hirosaki University, 3 Bunkyo-cho, Hirosaki, Aomori 036-8561, Japan e-mail<EMAIL_ADDRESS>
Bacteria and virulent bacteriophages are in a prey–predator relationship. Experimental models under simplified conditions with the presence of bacteria and bacteriophages have been used to elucidate the mechanisms that have enabled both prey and predator to coexist over long periods. In experimental coevolution conducted with Escherichia coli and the virulent RNA bacteriophage Qβ in serial transfer, both coexisted for at least for 54 days, during which time they continued to change genetically and phenotypically. By day 16, an E. coli strain partially resistant to Qβ appeared and caused an approximately 104-fold decrease in Qβ amplification. Whole-genome analysis of this strain suggested that a single mutation in TraQ was responsible for the partially resistant phenotype. TraQ interacts with propilin, encoded by the traA gene and a precursor of pilin, which is a component of the F pilus. The present study was performed to elucidate the mechanism underlying the coexistence of E. coli and Qβ by investigating how a mutation in TraQ altered the physiological state of E. coli, and thus the amplification of Qβ. Overexpression of wild-type TraQ in the partially resistant E. coli strain resulted in recovery of both TraA protein content, including propilin and pilin, and Qβ amplification to levels comparable to those observed in the susceptible strain. Intriguingly, overexpression of the mutant TraQ in the partially resistant strains also increased the levels of TraA protein and Qβ amplification, but these increases were smaller than those observed in the wild-type strain or the partially resistant strain expressing wildtypeTraQ.The results of this study represent an example of how E. coli can become partially resistant to RNA bacteriophage infection via changes in a protein involved in maturation of a receptor rather than in the receptor itself and of how E. coli can stably coexist with virulent RNA bacteriophages.

**Keywords: coevolution, prey–predator, experimental evolution, virulent phage, partial resistance**

#### **INTRODUCTION**

There have long been ecological and theoretical investigations regarding why predators do not eradicate their prey (Murdoch and Oaten, 1975; Anderson and May, 1978; Alexander, 1981; Berryman, 1992; Abrams, 2000; Briggs and Hoopes, 2004; Pettorelli et al., 2011; Brockhurst and Koskella, 2013). Bacteria and bacteriophages have been used as model systems of prey–predator interactions to elucidate the fundamental mechanism underlying this issue. Due to their short generation time, large population size, and ease of analyzing phenotypic and genomic changes, theoretical and experimental evolutionary research have been conducted extensively using these systems (Campbell, 1961; Levin et al., 1977; Lenski and Levin, 1985; Lenski, 1988; Schrag and Mittler, 1995; Sasaki, 2000; Paterson et al., 2010; Dennehy, 2012; Brockhurst and Koskella, 2013; Bull et al., 2014).

The existence of refuges for sensitive bacteria or the occurrence of endless arms races were suggested to be necessary to explain the coexistence of bacteria and bacteriophages (Chao et al., 1977; Lenski, 1988; Schrag and Mittler, 1995; Brockhurst et al., 2006). The numerical refuge in which density-dependent protection of susceptible cells from over-predation (Chao et al., 1977), spatial refuges such as wall populations on flasks in continuous culture or solid media used in serial passage (Chao et al., 1977; Lenski, 1988; Schrag and Mittler, 1995), and physiological refuges in which cells become transiently resistant or susceptible have been discussed (Lenski, 1988). In most of these previous studies, DNA bacteriophages, such as T2, T5, T7, λ, and 2, were used (Bohannan and Lenski, 2000; Buckling and Rainey, 2002; Dennehy, 2012). Although a great deal of knowledge has been accumulated regarding DNA bacteriophages, little is known about RNA bacteriophages in terms of stable coexistence.

In our previous study, *Escherichia coli* and the lytic RNA bacteriophage Qβ (Qβ) were shown to coevolve for at least 54 days, equivalent to 165 generations, under conditions of serial passage with shaking (Kashiwagi and Yomo, 2011). Qβ is a bacteriophage with a single-stranded RNA genome that specifically infects and lyses *E. coli* cells to release progeny phages (Van Duin and Tsareva, 2006). Phenotypic and genomic analyses indicated the coevolution of both *E. coli* and Qβ. In the course of experimental coevolution, partially resistant *E. coli* appeared in the 54th generation (16th day). Detailed analysis of the partially resistant *E. coli* is necessary to determine how *E. coli* and Qβ coexist in this experimental coevolution system. Genetic analysis revealed a single mutation in the *traQ* gene in the day-16 *E. coli* population (Kashiwagi and Yomo, 2011). It was reported that TraQ protein is a chaperone for insertion of propilin encoded by the *traA* gene (Moore et al., 1982; Kathir and Ippen-Ihler, 1991), and propilin was also reported to be unstable in *traQ*− cells (Maneewannakul et al., 1993). The 13 kDa propilin is processed by peptidase to a 7-kDa pilin and pilin proteins are assembled into filaments (i.e., the F pilus). Qβ adsorbs the F pilus of *E. coli* at the first step of infection (Van Duin and Tsareva, 2006), and the adsorption rate of Qβ on partially resistant cells estimated by first-order kinetics decreased markedly (Kashiwagi and Yomo, 2011). Therefore, the partially resistant phenotype of *E. coli* to Qβ infection observed in coevolution may be correlated with F pilus biosynthesis, especially TraA, and we focused on the relationships among mutation in TraQ, TraA content, and Qβ amplification.

Here, we report that a single amino acid change in TraQ was linked to the reduction of TraA content in the *E. coli* population. In addition, this decrease was recovered by supplying ancestral (wild-type) or mutant-type TraQ from an expression vector, and the ability of Qβ to amplify in the cell also recovered. These results represent one example of how *E. coli* can become partially resistant to RNA bacteriophage infection, which involves changes in a protein related to the maturation of a receptor, in this case the F pilus, rather than changes to the receptor itself. These results suggest that the mutation in TraQ may cause heterogeneity within the *E. coli* population, with a small number of cells supporting the phage population and a large number of cells supporting the *E. coli* population without Qβ infection, even though the *E. coli* cells were genetically identical.

### **MATERIALS AND METHODS**

#### **STRAINS, CULTURE MEDIA, AND PLASMID DNA**

*Escherichia coli* Anc(C), the partially Qβ infection-resistant mutant strain, M54(C) (Kashiwagi and Yomo, 2011), and DH1 *leuB*::(*gfpuv5-Kmr*) (hereafter called LKG; Kishimoto et al., 2010) were used to characterize the effects of S21P mutation in TraQ protein and a control F− strain. *E. coli* A/λ (Watanabe et al., 1979) was used as the host strain for titer assay. LB medium (10 g/L tryptone, 5 g/L yeast extract, 10 g/L NaCl; Nakalai Tesque, Kyoto, Japan) was used.

To construct TraQ with a Strep-tag II sequence (Schmidt and Skerra, 2007) at the C-terminus of the ancestral-type TraQ (TraQAnctag) and mutant-type TraQ (TraQS21Ptag), the *traQ* gene was amplified by PCR with Anc(C) and M54(C) genome as the template, the primers traQ\_XbaI and traQ\_strep\_HindIII, and Phusion® High-Fidelity DNA polymerase (New England Biolabs, Ipswich, MA, USA). The oligonucleotide DNAs used in this study are listed in **Table 1**. The *traQ*Anctag and *traQ*S21Ptag genes with *Xba*I and *Hin*dIII sites at both ends were subcloned into the *Xba*I/*Hin*dIII sites of pASK-IBA3plus (IBA Biologics GmbH, Goettingen, Germany). The resulting plasmid DNAs were designated as pASK-traQAnctag and pASK-traQS21Ptag, respectively.

#### **Table 1 | Oligonucleotide DNA sequence.**

#### **DNA SEQUENCING OF THE** *traQ* **GENE OF ANCESTRAL AND PARTIALLY RESISTANT** *E. coli*

To determine the *traQ* gene sequneces of Anc(C) and M54(C), the *traQ* region in 10 colonies each of Anc(C) and M54(C) was amplified by PCR with the primers F\_4f\_2 and F\_4r\_2 and Phusion® High-Fidelity DNA polymerase (New England BioLabs), and PCR products were directly sequenced by the dideoxynucleotide chain termination sequencing method (Sanger et al., 1977).

#### **ESTIMATION OF Qβ AMPLIFICATION**

Anc(C)/pASK-traQAnctag, Anc(C)/pASK-traQS21Ptag, Anc(C)/ pASK-IBA3plus, M54(C)/pASK-traQAnctag, M54(C)/pASKtraQS21Ptag, and M54(C)/pASK-IBA3plus were cultured in 5 mL of LB with 100 μg/mL ampicillin overnight and 50 μL of the culture was inoculated into 5 mL of LB with 100 μg/mL ampicillin for approximately 2 h. Aliquots of 1 mL of the culture were transferred into 4 mL of LB medium with 100 μg/mL ampicillin and 100 nM doxycycline-HCl (Dox) and cultured for a further 2 h. Qβ was added to infect the bacterial cells and the free phage was separated immediately or 4 h after infection by centrifugation at 13400 × *g* for 1 min. The free phage in the supernatant was diluted, the number of plaque forming units per milliliter was determined (PFU/mL), and the amplification ratio was calculated as *x* = N4/N0, where N4 and N0 represent 4 h after and initial (0 h) free phage density, respectively. The titer assay was conducted according to the standard method described previously (Carlson, 2005).

#### **WESTERN BLOTTING ANALYSIS**

A polyclonal antibody to TraA protein raised against the keyhole limpet hemocyanin-conjugated peptide (CDLMASGNTTVKAT-FGKDSS) was obtained from Sigma-Aldrich Japan (Tokyo, Japan).

Cell preculture was conducted as described in the Section "Estimation of Qβ Amplification." For induction with 100 nM Dox, Dox was added to the 2-h culture and cells were cultured for a further 5.75 h. Without Dox induction, the cells were cultured for 7.25 h in LB medium. Proteins from the cells obtained from 0.1 mL of culture with OD600 = 2.0 were subjected to SDS-PAGE using Any kDTM Mini-PROTEAN® TGXTM precast gels (Bio-Rad Laboratories, Hercules, CA, USA), and the TraA protein was determined by Western blotting analysis with anti-TraA antibody and HRP-conjugated goat anti-rabbit polyclonal antibody (Santa Cruz Biotechnology, Inc., Santa Cruz, CA, USA) as the primary and secondary antibodies, respectively, diluted with Can Get Signal® immunoreaction enhancer solution (Toyobo Co. Ltd., Osaka, Japan). TraQ with a Strep-tag II at the C-terminus was detected using Precision ProteinTM Strep-Tactin-HRP conjugate (Bio-Rad Laboratories). The signals were detected with Chemi-Lumi One L (Nakalai Tesque).

#### **NORTHERN HYBRIDIZATION ANALYSIS**

To compare the *traA* mRNA contents by Northern hybridization, Anc(C), M54(C), and LKG with pASK-IBA3plus were cultured in LB medium and total RNA was extracted from the cells in logarithmic phase using the SV Total RNA isolation system (Promega, Madison, WI, USA) according to the manufacturer's instructions. Total RNA from cells obtained from 0.16 mL of culture with OD600 = 0.27 was subjected to Northern hybridization using a digoxigenin-labeled single-stranded DNA probe, i.e., traA1 for *traA* and 16SrRNA\_1 for 16SrRNA, as an indicator of the amount of total RNA used (Moran et al., 1995). The signals were detected with CDP-Star (GE Healthcare UK Ltd., Little Chalfont, UK). DynaMarker® Prestain Marker for RNA High (BioDynamics Laboratory Inc., Tokyo, Japan) was used to obtain the standard curve for RNA length.

#### **DETERMINATION OF 5- - AND 3- -TERMINAL SEQUENCES OF** *traA* **mRNA**

To determine the 5- -terminal sequence of *traA* mRNA, cDNA was synthesized using total RNA of Anc(C) as the template, Super-Script® III reverse transcriptase (Life Technologies, Carlsbad, CA, USA), and traA\_r primer. The 5- -phosphorylated DNA linker 5PpACYC\_rev was ligated at the 3- -terminus of the first strand cDNA with T4 RNA ligase 1 (New England Biolabs). PCR was performed using the resultant cDNA as the template, PrimeSTAR® HS DNA polymerase (Takara Bio Inc., Shiga, Japan), and the primers pACYC\_rev2 and traA\_r2. To determine the 3- -terminal sequence of *traA* mRNA, universal miRNA cloning linker (5- rAppCTGTAGGCACCATCAAT-NH2-3- ; New England Biolabs) was ligated with the 3- -terminus of total RNA of Anc(C) using T4 RNA ligase 1 (New England Biolabs). The first strand cDNA was synthesized using the primer Linker\_r and SuperScript® III reverse transcriptase (Life Technologies), and then purified cDNA was subjected to PCR using the primers Linker\_r and traA\_f. Three and two bands of PCR products for 5- - and 3- -terminal sequence determination, respectively, were sliced from the gel and subcloned using a Zero Blunt® TOPO® PCR Cloning Kit for Sequencing (Life Technologies). Three to six clones were randomly picked and sequenced by the dideoxynucleotide chain termination sequencing method.

#### **STATISTICAL ANALYSIS**

Amplification ratios were compared by one-way ANOVA with the *post hoc* Tukey's test (Zar, 2010). In all analyses, values of log10 (N4/N0) of each strain were used for statistical analysis. The Studentized range, *q*, is shown in the text. In all analyses, *P* < 0.01 was taken to indicate statistical significance.

#### **RESULTS**

#### **RECOVERY OF Qβ AMPLIFICATION IN RESISTANT** *E. coli* **BY SUPPLYING TraQ**

We first analyzed the *traQ* gene sequences from 10 single colonies derived from the coevolved *E. coli* population to confirm that the majority harbored the T61C mutation. We picked 10 single colonies from each of Anc(C) and M54(C) populations, which were the initial and day-16 *E. coli* populations in the coevolution experiment (Kashiwagi and Yomo, 2011). All 10 colonies of Anc(C) had T and all 10 colonies of M54(C) had C at position 61, and this T61C mutation resulted in S21P in the TraQ protein.

We analyzed the amplification ratio of Qβ on Anc(C) and M54(C) that harbored only the vector (pASK-IBA3plus) to determine the extent of reduction in the amplification ratio of Qβ on M54(C). Anc(C) and M54(C) expressed inherent ancestraltype TraQ protein and mutant-type TraQ protein from the F plasmid, respectively. The amplification ratio of Qβ was calculated as described in the Section "Materials and Methods." The amplification ratios of Qβ on Anc(C)/pASK-IBA3plus and M54(C)/pASK-IBA3plus were 3.5 × 104 and 1.8, respectively (**Figure 1**). Although the amplification ratio on M54(C) was much lower than that on Anc(C), Qβ could undergo amplification on M54(C), indicating that M54(C) had a partially rather than fully resistant phenotype.

To analyze whether ancestral-type TraQ expression in M54(C) could compensate for Qβ amplification in this strain, ancestraltype TraQ with the Strep-tag II sequence at the C-terminus