The expansion of a pan-genome can be examined by plotting the number of genomes considered against the total number of genes observed. The associations on the plot can then be mathematically evaluated through fitting regression models (Tettelin et al., 2008), according to which it can be classified as “open” or “closed” (see Materials and Methods for details). While the size of an open pan-genome would increase unboundedly with the inclusion of new genomes, the size of a closed pan-genome would reach a plateau after a certain number of sample genomes were included. Previously, such analyses have only been performed at the species level (Tettelin et al., 2008). In order to compare the pan-genomes of different taxonomic ranks, we constructed regression models for five groups of Epsilonproteobacteria that correspond to three different levels of taxonomy (Table S1): two groups representing strains of same species (intra-species), two representing species of same genera (intra-genus), and one representing different species of the entire class (intra-class or inter-species). It is worth mentioning that groups of higher taxonomic rank encompassed those of lower rank. For example, the inter-species group (Class) included all the different species in the intra-genus groups (Genus), and the intra-genus groups included a representative strain from each of the intra-species groups (Species). Hence, our analyses compared the extent of pan-genome expansion at all three taxonomic levels.

We used a two-step process to model the expansion of the pan-genomes. First, complete or random permutations were carried out with step-wise additions of new genomes, and a median was taken on the size of pan-genomes after each step. Second, the median counts was extrapolated using two different models: power-law regression (Tettelin et al., 2008) and exponential regression (Tettelin et al., 2005). The resulting extrapolations were normalized by the median genome sizes in their respective sets to assist the visualization and comparison of the fitted curves (Figure 1).

The power-law regression showed that all five pan-genomes in our dataset are open, with an average γ parameter of 0.29, 0.53, and 0.66, respectively, for intra-species, intra-genus, and inter-species sets (Materials and Methods). The extrapolated curves of intra-species (red), intra-genus (blue), and inter-species (green) pan-genomes followed distinct slopes: while the intra-species curves were the shallowest, the inter-species curves were the steepest. In contrast, the curves within the same category (the two intra-species or the two intra-genus groups) presented similar slopes (Figure 1A and Table S2). While both intra-species curves (H. pylori and C. jejuni) were slightly shallower than the pan-genome curves of E. coli (Rasko et al., 2008) and SAR11 (Grote et al., 2012) and steeper than that of Streptococcus agalactiae (Tettelin et al., 2005), the two intra-genus curves (Helicobacter and Campylobacter) were much steeper than all evaluated intra-species curves (Figure 1A).

The power-law model fit the available data with high R² values of more than 0.98 for all extrapolations, while the exponential model had on average lower R² values, especially when the number of considered genomes is relatively small (Table S2). To further evaluate this, we fitted both models to a varying number of data points to monitor how the number of available genomes may affect the accuracy of pan-genome modeling (Figures 1C,D). For example, when using a range of 6–10 considered genomes in the modeling of H. pylori (Figure 1C), the exponential curves (dotted lines) are more spread out than the power-law curves (solid lines). Similarly, modeling of the entire Epsilonproteobacteria also showed that the power-law model is more stable than the exponential regression model and less influenced by the availability of genomic data (Figure 1D). Besides this, results in Figures 1C,D also showed that the number of available genomes are far from saturating the power-law model. In other words, new genes are still been discovered with the addition of each genome. Therefore, the diversity of epsilonproteobacterial pan-genomes is yet to be fully explored with additional genomic sequences.

The pan-genome of all the examined Epsilonproteobacteria contains 16,349 clusters of non-redundant protein coding genes (Materials and Methods). Among them, 289 clusters (1.7%) represented Conserved Single Copy Genes (CSCGs) that appear only once in every examined genome. These CSCG clusters contain 11,271 genes, which account for 16% of the more than 70,000 genes in all analyzed genomes. At the other end of the spectrum are 10,944 clusters (67%) that occur only in a single genome, accounting for about 15% of the total genes. We classified the total genes in the pan-genomes into three sets based on their occurrence: (1) the “core” genes, which are universally conserved in all considered genomes, include both CSCGs and conserved genes of multiple copies per genome, (2) the “singleton” genes are specific to single genomes, and (3) the “flexible” genes are found in more than one, but not all genomes.

The functional distribution of core, flexible, and singleton genes was examined for the intra-genus groups of Helicobacter, Campylobacter, and Sulfurimonas, as well as the inter-species group. Figure 2 shows the classifications based on the Clusters of Orthologous Groups (COG) database (Tatusov et al., 2003). The genes that cannot be classified into any existing COG clusters were grouped into the “unmapped” category. For all groups, the majority of the singleton genes (56–79%) and a large fraction of the flexible genes (38–48%) were unmapped or poorly characterized, while only a small fraction of the core genes (13–27%) had no clear functional assignment. The core, which represents an indispensable part of all genes in Epsilonproteobacteria, contained mainly housekeeping genes that encode the central machinery of a cell, such as translation, protein modification and turnover, replication and repair, cell wall biogenesis, as well as co-enzyme metabolism. In contrast, the flexible and singleton genes that could be mapped to a functional category mainly encoded functions in signal transduction, inorganic ion transport, and energy production.

While both considered to be non-core (Medini et al., 2005), the flexible and singleton genes presented slightly different functional distributions. A larger fraction of the flexible genes (66–78%) were mapped to existing COG families than the singleton genes (33–55%) despite both being poorly characterized. Additionally, a slightly larger fraction of the flexible genes (10–14%) encoded functions in energy production and inorganic ion transport than the core (5–12%) and singleton (3–7%) genes. Therefore, these processes were more likely shared among subgroups of Epsilonproteobacteria than being common to all or unique to a single species. In the following two sections, we further investigate the composition and evolution of pan-genomes by analyzing the core and flexible genes. First, phylogeny of the core genes was reconstructed and used to infer the evolutionary history of Epsilonproteobacteria. Second, an analysis was carried out to identify lineage-specific and niche-specific signals within the flexible genome of selected species.

The core genome of Epsilonproteobacteria contains 289 CSCGs and 40 conserved genes with multiple copies per genome. Specifically, the CSCGs accounted for 15% of all the genes in an average epsilonproteobacterial genome. The fact that these genes are universally present in a single copy makes them useful markers for inferring the phylogenetic relationships within the class (Figure 3A), as well as evaluating the phylogenetic position of the Epsilonproteobacteria as a whole (Figure 4).

The phylogenomic analysis (Figure 3A) is overall consistent with a phylogenetic tree based on 16S rRNA sequences (Figure 3B), in which species of the same taxonomic groups were tightly clustered. According to both trees, the deepest lineages within the Epsilonproteobacteria are represented by Nautilia profundicola (Smith et al., 2008) and Caminibacter mediatlanticus (Voordeckers et al., 2005), both of which are moderate thermophiles and obligate anaerobes isolated from deep-sea hydrothermal vents. The other species from either vent or non-vent environments emerged later in evolution, which paralleled the emergence of host-associated species. Despite these similarities, the CSCG-based tree (Figure 3A) was overall better resolved compared to the 16S rRNA gene tree (Figure 3B), and it indicated slightly different branching patterns of certain clades, such as the precise phylogenetic positions of the genera Arcobacter and Helicobacter.

Besides assisting in the evolutionary reconstruction within the Epsilonproteobacteria, the CSCGs can also help in evaluating the phylogenetic affiliation of this class as a whole (Figure 4). We used a set of 37 phylogenomic markers that are universally present in single copies in a set of fully sequenced genomes to build the bacterial tree. These markers included 31 universal genes published in a previous study (Wu and Eisen, 2008), as well as six additional genes that were identified through a global search of epsilonproteobacterial CSCGs using Hidden Markov Models (HMMs). The concatenated protein tree revealed a close affiliation of the Epsilonproteobacteria with the Aquificae, as well as provided evidence that the Epsilonproteobacteria represent a distinct clade that is separated from Proteobacteria. The global analyses also positioned Deltaproteobacteria into a distal branch. With the exception of Hippea maritima (Anderson et al., 2011), all other examined deltaproteobacterial genomes grouped with the phyla Thermodesulfobacteria (Anderson et al., 2012; Elkins et al., 2013), Acidobacteria (Ward et al., 2009; Challacombe et al., 2011; Rawat et al., 2012), and Nitrospirae (Lücker et al., 2010; Fujimura et al., 2012).

While a subset of the core genes supported the phylogenomic reconstruction of Epsilonproteobacteria, they only covered a small fraction of an average genome. The remaining genes (85%) were either shared among a subset (flexible genes), or they only occurred in one genome (singleton genes) of currently sequenced Epsilonproteobacteria. Therefore, analyses of these genes can provide a more complete picture of evolution of this class. Since the singleton genes were largely unmapped to COG families (Figure 2) and because the definition of singleton genes is subject to change, i.e., what appears to be a singleton today may be reclassified into a multi-gene cluster when new genome data becomes available in the future, we decided to focus our study on the flexible genes.

We hypothesized that the flexible genome carries two types of information: (1) lineage specific signals that trace the evolutionary heritage of a single lineage, and (2) niche specific signals that represent the adaptation of different species to similar environments. We tested this hypothesis by comparing the genome of a free-living species from a particular environment to genomes of species within the same genus, but inhabiting different environments, as well as to genomes of other species inhabiting a similar environment, but belonging to different genera. We chose Sulfurimonas autotrophica for our studies, which was the only free-living species that matched our criteria. We compared the genome of S. autotrophica with two subsets (Figure 5), one containing its phylogenetic neighbors of the same genus that were isolated from different environments (S. gotlandica, and S. denitrificans), and the other containing organisms with similar habitats, but belonging to different genera (Sulfurovum sp. NBC37-1, and Nitratiruptor sp. SB155-2) (Figure 3).

Our analyses successfully identified genes that were either uniquely shared among the Sulfurimonas species, i.e., genus specific, or among the vent isolates, i.e., habitat specific (Figure 5). This confirmed our hypothesis that both lineage-specific and niche-specific evolutionary events were recorded in the flexible genome. Overall, there were almost three times as many genes specifically inherited within the genus of Sulfurimonas (195 unique genes, Table S3) than shared among the vent isolates (67 unique genes, Table S4). Further, the Sulfurimonas-specific genes were distributed evenly across the genome of S. autotrophica, while the vent-specific genes clustered tightly at certain locations of the genome. The vent-specific gene clusters encode three major functions, including glycogen production and utilization, inorganic ion and efflux transport, and energy production and conversion. We specifically focused on two vent-specific protein complexes related to energy production and conversion (Hyf-like and PYDH complex) as examples to further elucidate the evolution of vent-specific genes in the flexible genome.

The vent-specific hyf-like operon (hyfBCEFGI) encodes a putative hydrogenase-4 complex (Figure 5B). This operon is conserved among all the vent species, while the genomic neighborhood of the operon is completely unrelated between different species (Figure 5B). Besides the hyf-like operon, all examined genomes of vent-inhabiting Epsilonproteobacteria also encode a homolog of formate dehydrogenase H (Fdh-H). This combination resembles the formate hydrogenlyase (FHL-2) complex of E. coli, which oxidizes formic acid to carbon dioxide and molecular hydrogen (Andrews et al., 1997). The pyruvate dehydrogenase complex (PYDH) is a multi-enzyme complex composed of three different enzymes: a decarboxylase (E1p), a dihydrolipoamide acyltransferase (E2p), and a dihydrolipoamide dehydrogenase (LPD), among which the E1p carries out the initial step of an enzymatic reaction that converts pyruvate to acetyl-CoA, NADH and CO₂ (Neveling et al., 1998). The E1p can be present in two forms, one containing multiple copies of a single subunit (type-I) and the other containing multiple copies of two subunits, alpha and beta (type-II). These two forms are not evolutionarily related, and either or both forms may be present in the same organism (Schreiner et al., 2005). Our analyses confirmed the presence of PYDH in seven epsilonproteobacterial genomes: two encode the type-I form and belong to the genus Arcobacter, while the other five encode the type-II form and are from species inhabiting deep-sea hydrothermal vents (Table S5). We constructed phylogenetic trees of the E1p enzymes in order to further investigate the evolution of the two forms of PYDH in Epsilonproteobacteria (Figure 6). According to the phylogenetic trees, the type-I form, which existed solely in the genus Arcobacter, was closely related to the PYDHs encoded in Gammaproteobacteria (Figure 6A), whereas the type-II form, which existed solely in the vent species, grouped with the enzymes of Sulfurihydrogenibium azorense, which belongs to the Aquificales, and of Geobacter species (Figure 6B).

A total of 39 published genomes were collected for this study (Table S1), the majority of which came from human or animal associated species of clinical relevance, with 16 genomes representing different strains of two widely studied pathogens, Helicobacter pylori and Campylobacter jejuni. Despite this bias toward clinical species, nine complete genomes in our dataset represented free-living species from multiple genera that were isolated from a variety of environments, including deep-sea hydrothermal vents, coastal marine sediments, oil fields, and salt marshes. Additionally, we also included near-complete draft genomes from two free-living species, Caminibacter mediatlanticus TB-2 and Sulfurimonas gotlandica GD1, in order to increase the number of genomes from non-pathogenic species.

In order to identify non-redundant gene clusters, we performed pairwise comparison on the collected genomes to identify bi-directional best hits of the encoded proteins. The determination of bi-directional best hits was based on the BLAST software (Altschul et al., 1990, 1997) using three criteria: (1) the e-value should be better than 0.001; (2) the alignment should cover at least 70% of the sequences; and (3) the sequence similarity should be better than 30%. Here the sequence similarity cutoff may seem low, but we used it to accommodate the orthologs in distantly related species. We showed that for the majority of the multi-gene clusters we identified, at least 90% of genes in these clusters encode proteins of the same Pfam family assignment, COG family assignment and functional annotation based on exact word matches (Figure 7). This is significant considering the many different existing ways to name identical functions. Figure 7 only presented a lower estimation of the accuracy of our approach in obtaining coherent orthologous clusters.

We performed permutations on five genome sets: one inter-species, two intra-genus, as well as two intra-species sets (Table S1). At the inter-species level, we examined all of the 25 distinct species in our genome set of Epsilonproteobacteria; at the intra-genus level, we looked at the different species within the genera Helicobacter (6 species) and Campylobacter (6 species), disregarding strain variations over the same species by selecting one representative strain each for H. pylori and C. jejuni; at the intra-species level, we compared the different strains of H. pylori (10 strains) and C. jejuni (6 strains). All possible permutations were explored for the intra-species and intra-genus datasets. However, the inter-species dataset contains too many genomes to be fully permuted. Alternatively, we created 25,000 random permutations by randomly selecting genomes one after another over the whole dataset until all genomes has been visited.

Through each permutation, a data vector is produced to record the increased number of unique genes in the pan-genome with the step-wise addition of new genomes. Then, medians and the 25th and 75th percentile values were calculated over all permutations for each data point of the vectors. Finally, these counts were normalized by the median size of genomes, respectively (Table S2), so that pan-genomes of different datasets can be compared with one another (Figure 1).

The power law regression was performed following the approach described in Tettelin et al. (2008). The regression function n = σN^γ was used to model the median sizes of the pan-genomes generated from all permutations, where n is the total number of non-orthologous genes in the pan-genome, N is the number of genomes considered, and σ and γ are free parameters (Table S2). When 0 < γ < 1, the pan-genome is considered open because it is an unbounded function over the number of genomes. When γ < 0, the pan-genome is considered closed since it approaches a constant as more genomes are considered.

Instead of fitting the total number of genes in pan-genomes, the exponential regression model fits the number of new genes per added genome, which were implemented with an exponential decay function F_s(N) = K_s exp[−N/T_s] + tg(θ), where F_s(N) is the number of new genes with the addition of each new genome, tg(θ) is the same number when N approaches infinity, and K_s, T_s, and tg(θ) are free parameters (Tettelin et al., 2005). Based on the estimated parameters in the exponential decay function, the pan-genomes were modeled with the formula Pan(N)=D+∑j=2N{ Ksexp[ −j/Ts ]+tg(θ) }, where D is the median number of genes per sequenced genome at each dataset (Table S2).

We adapted the protocol by Wu et al. for phylogenomic reconstructions (Wu and Eisen, 2008). In a first step, the individual clusters of CSCG-encoded proteins were aligned using MUSCLE (Edgar, 2004), and HMMs were built for each cluster using hmmbuild from the HMMER package (Eddy, 2011). Then, the models were used as queries to search against other genomes and the resulting alignments were trimmed adapting scripts from AMPHORA (Wu and Eisen, 2008). In a next step the trimmed alignments were concatenated with one another into a master alignment, which was further refined using Gblocks (Talavera and Castresana, 2007) to remove the less conserved columns. Finally, the refined master alignment was used as the input for PhyML (Guindon et al., 2010) for phylogenetic reconstruction.

The CSCG tree of Epsilonproteobacteria (Figure 3A) included six additional draft or complete genomes that were published after our initial steps of data collection. These included Sulfurospirillum barnesii SES-3, Uncultured Sulfuricurvum sp. RIFRC-1, Arcobacter butzleri ED-1 (Toh et al., 2011), Arcobacter sp. L (Toh et al., 2011), Sulfurovum sp. AR (Park et al., 2012), as well as the single-cell genomes of Thiovulum sp. ES (Marshall et al., 2012). To accommodate the incompleteness of draft genomes, we selected a subset of the CSCG-encoded proteins that occurred once in every draft genomes, and used only these as markers for tree construction. As a result, 194 of the CSCG-encoded proteins were used in the above procedure to construct the local phylogeny for Epsilonproteobacteria.

The global bacterial phylogeny was constructed with 37 globally conserved single copy markers (Figure 4). In addition to the 31 applied in the AMPHORA package (Wu and Eisen, 2008), we identified six additional phylogenetic markers using the HMM of core proteins: DNA gyrase subunit B (gyrB), Tryptophanyl-tRNA synthetase (TrpRS), SSU ribosomal protein S12p (S23e), LSU ribosomal protein L17p, SSU ribosomal protein S4p (S9e), and SSU ribosomal protein S15p (S13e). Among these new marker genes, GyrB (Kasai et al., 2000; Holmes et al., 2004; Peeters and Willems, 2011) and TrpRS (Rajendran et al., 2008) have been used in previous studies to determine the phylogeny of selected taxonomic groups, and the rest are ribosomal proteins.

The global bacterial tree in Figure 4 was rooted using mid-point rooting. The 16S and CSCG trees in Figure 3 were rooted based on the relative positions of different epsilonproteobacterial species at the global bacterial tree and using all other bacteria as an outgroup. As indicated with a black arrow in Figure 4, the root of Epsilonproteobacteria is located between Nautiliales and the other examined lineages.

We collected representative sequences of the type-I and type-II forms of PYDH based on annotations in the UniProt database (Bairoch et al., 2005). The protein phylogeny was reconstructed with PhyML (Guindon et al., 2010) using default settings. The tree of type-I form PYDH (Figure 6A) is rooted using mid-point rooting, and the tree of type-II form PYDH (Figure 6B) is rooted with the Eukaryotes as an outgroup.

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.