High-throughput single-cell sequencing identifies photoheterotrophs and chemoautotrophs in freshwater bacterioplankton

Abstract

Recent discoveries suggest that photoheterotrophs (rhodopsin-containing bacteria (RBs) and aerobic anoxygenic phototrophs (AAPs)) and chemoautotrophs may be significant for marine and freshwater ecosystem productivity. However, their abundance and taxonomic identities remain largely unknown. We used a combination of single-cell and metagenomic DNA sequencing to study the predominant photoheterotrophs and chemoautotrophs inhabiting the euphotic zone of temperate, physicochemically diverse freshwater lakes. Multi-locus sequencing of 712 single amplified genomes, generated by fluorescence-activated cell sorting and whole genome multiple displacement amplification, showed that most of the cosmopolitan freshwater clusters contain photoheterotrophs. These comprised at least 10-23% of bacterioplankton, and RBs were the dominant fraction. Our data demonstrate that Actinobacteria, including clusters acI, Luna and acSTL, are the predominant freshwater RBs. We significantly broaden the known taxonomic range of freshwater RBs, to include Alpha-, Beta-, Gamma- and Deltaproteobacteria, Verrucomicrobia and Sphingobacteria. By sequencing single cells, we found evidence for inter-phyla horizontal gene transfer and recombination of rhodopsin genes and identified specific taxonomic groups involved in these evolutionary processes. Our data suggest that members of the ubiquitous betaproteobacteria Polynucleobacter spp. are the dominant AAPs in temperate freshwater lakes. Furthermore, the RuBisCO (ribulose 1,5-bisphosphate carboxylase/oxygenase) gene was found in several single cells of Betaproteobacteria, Bacteroidetes and Gammaproteobacteria, suggesting that chemoautotrophs may be more prevalent among aerobic bacterioplankton than previously thought. This study demonstrates the power of single-cell DNA sequencing addressing previously unresolved questions about the metabolic potential and evolutionary histories of uncultured microorganisms, which dominate most natural environments.

Figure 1

Principal coordinates analysis of weighted UniFrac pairwise distances between 16S rRNA gene sequences from the studied environmental samples and the HNA and LNA cell fractions. A neighbor-joining tree (Jukes–Cantor substitution model) including all 16S rRNA gene sequences from SAGs served as the input data for the Fast UniFrac analysis. The archaeon Nitrosopumilus maritimus (CP000866) was used as an outgroup.

Figure 2

The relative frequency of photoheterotrophs in freshwater bacterioplankton, as determined by single-cell approach and metagenomic sequencing. The frequency of photoheterotrophs among SAGs was determined as the fraction of 16S rRNA-positive SAGs from which rhodopsin or pufM gene was recovered. The frequency of photoheterotrophs in metagenomes was determined as the ratio of either rhodopsin or pufM to recA+radA. M, metagenomics; SCG, single-cell genomics; Dam, Damariscotta Lake.

Figure 3

Maximum-likelihood tree of 119 rhodopsin proteins from single cells: (a) general tree; (b) subtree of actinorhodopsin-like proteins and (c) subtree of proteorhodopsin-like proteins. In all, 120 amino-acid positions were used in the tree construction. Bootstrap values ⩾50 are displayed. SAGs obtained from the HNA and LNA cells are indicated in regular and italic fonts, respectively. The taxonomic identity of rhodopsin-containing SAGs, based on their 16S rRNA gene phylogeny, is provided next to the SAG name (for detailed phylogeny of the 16S rRNA genes, see Supplementary Figure S2). (d) Recombination analysis of rhodopsin genes. The dual multiple change-point model that considers the spatial variation of tree topologies and the substitution process parameters was applied in a Bayesian framework using reversible jump Markov chain Monte Carlo sampling to approximate the joint posterior distribution of all model parameters. Parameters of transition:transversion (κ) and expected divergence (μ) and spatial variation of tree topologies are indicated. Each one of the breakpoints shown in the tree topologies together with κ and μ parameters indicate a putative recombination event. Recombination was not detected within the actinorhodopsin and proteorhodopsin clusters.

Figure 4

Maximum-likelihood phylogenetic analysis of pufM and RuBisCO genes and the corresponding 16S rRNA sequences from single cells (limited to Betaproteobacteria owing to space constrains). The taxonomic identity of the pufM- and RuBisCO-containing SAGs is indicated next to the SAG name. BchlY-containing SAGs are indicated in the phylogenetic tree of 16S rRNA gene (for detailed phylogeny of the 16S rRNA genes, see Supplementary Figure S2). Bootstrap (1000 replicate) values ⩾50 are displayed. In the case of pufM gene, the analysis was conducted on nucleotide sequence alignment (250 nucleotide positions). For the RuBisCO gene, the analysis was based on protein sequence alignment (amino-acid positions 100–262).