The Chlamydiales Pangenome Revisited: Structural Stability and Functional Coherence

The entire publicly available set of 37 genome sequences from the bacterial order Chlamydiales has been subjected to comparative analysis in order to reveal the salient features of this pangenome and its evolutionary history. Over 2,000 protein families are detected across multiple species, with a distribution consistent to other studied pangenomes. Of these, there are 180 protein families with multiple members, 312 families with exactly 37 members corresponding to core genes, 428 families with peripheral genes with varying taxonomic distribution and finally 1,125 smaller families. The fact that, even for smaller genomes of Chlamydiales, core genes represent over a quarter of the average protein complement, signifies a certain degree of structural stability, given the wide range of phylogenetic relationships within the group. In addition, the propagation of a corpus of manually curated annotations within the discovered core families reveals key functional properties, reflecting a coherent repertoire of cellular capabilities for Chlamydiales. We further investigate over 2,000 genes without homologs in the pangenome and discover two new protein sequence domains. Our results, supported by the genome-based phylogeny for this group, are fully consistent with previous analyses and current knowledge, and point to future research directions towards a better understanding of the structural and functional properties of Chlamydiales.


Introduction
Members of the order Chlamydiales are obligate intracellular bacteria, characterized by a unique developmental cycle and are important pathogens of humans and animals resulting in a wide range of diseases, including several zoonoses [1][2][3]. The order Chlamydiales, separated from other eubacteria by forming a deep branch in ribosomal RNA-based phylogenetic trees, has been enriched by new lineages. Beside the family Chlamydiaceae, in which important chlamydial pathogens are grouped, new families, such as Parachlamydiaceae, Simkaniaceae and Waddliaceae, have been recognized to accommodate newly discovered pathogenic and non-pathogenic chlamydial organisms [4][5][6].
Since the release of the first chlamydial genome sequence from Chlamydia trachomatis (serovar D) [7], new genomes are being sequenced, thus offering insights into the genome organization and functional capacity of the corresponding species [8]. Besides its crucial importance for applied research in medical and veterinary microbiology [9], this corpus of genomic information is also key to understanding the evolutionary position of various chlamydial species (or strains) and the inference of the internal phylogeny of this distinct taxon [8,[10][11][12].
As the intracellular lifestyle imposes constraints on gene content and metabolic capabilities, the Chlamydiales might represent one of the best datasets for the development of pangenome analysis methods [13]. Additional challenges are the wide variety of chlamydial genome sizes with unequal rates of reduction, and a repertoire of less characterized proteins than other bacterial groups whose pangenomes have been analyzed, e.g., Streptococcus or Salmonella [14,15].
Previously, we have used the genome of Chlamydia trachomatis [7] as a case study for annotation transfer quality [16]. Using a novel encoding scheme and a scoring function called TABS for transitive annotation-based scale [16], our main finding regarding annotation was that, despite a number of inconsistencies, automated annotation pipelines performed remarkably well when benchmarked against a manually curated annotation corpus [16]. These results are important for the quantification of reproducibility and consistency in genome-wide annotation [17].
In this work, we explore the entire set of the Chlamydiales pangenome with a broad collection of genome sequences publicly available to date (31 Chlamydiaceae and six other Chlamydiales genomes), twice as many as in a similar recent analysis [18]. Importantly, our pangenome analysis pipeline Genes 2012, 3 incorporates recently sequenced genomes of key Chlamydiaceae species not previously reported, thus augmenting our understanding from previous findings [18,19].
We focus on key aspects of pangenome analysis and explore multiple facets of the Chlamydiales gene content in terms of protein-coding genes and families. We also provide certain key findings that might illuminate the evolutionary history of this group as well as interesting sequence motifs not widely shared within this order. Beyond the confirmation of the recent analysis of the Chlamydiales as mentioned above [18], we also use this group to expand on methods for pangenome analysis [13,20] by proposing a pangenome analysis pipeline. Our results are consistent with wider studies of pangenomes [21] and provide additional knowledge for Chlamydiales. In conclusion, pangenome analysis offers an opportunity for the study of bacterial genome evolution, the development of relevant methods and the understanding of genome structure and proteome function on a large scale.

Data Collection
All protein sequence data from 37 genomes were compiled into a single data collection (February-July 2011), including the most recent published Chlamydiales genomes. In total, 43,736 protein-coding genes were extracted from public databases corresponding to the entire set of 37 genome sequences from the bacterial order Chlamydiales currently available (Table 1). Sequence data were codified following the style of the COGENT database [22], for easy identification both by programs and human users (Supplement S1). The above notation is followed throughout this work. The COGENT scheme encodes genus and species names into a four-character identifier prefix string, followed by a code for the strain name, its version (in this collection all versions are considered as version 1 and optionally hidden) and finally for proteins the relative order of the sequence within the genome [23] (Table 1). We have also recorded the date of publication for the corresponding genome (or the release date where no publication was available) (Supplement S2).  CTRA-SWE-01  875  Total 43,736 The first column signifies the inclusion order into the genome collection and does not reflect any other relationship. The second column lists the species and strain name, the third column the COGENT-style identifier and the last column the number of protein-coding genes.

Sequence Comparison
All protein sequence data were masked using CAST with default parameters (threshold = 40), to exclude compositionally biased regions [24]. In total, 6,906 such regions were filtered out, provided for further study (Supplement S3).
The masked sequences were used as queries against the genome corpus, in an all-against-all mode with BLAST (blastall, e-value threshold 10 í6 ) [25,26]; in total, more than 40,000 BLAST searches were performed and 1,709,325 significant similarities below threshold were obtained (Supplement S4).

Clustering and Annotation
The similarity pairwise list (from Supplement S4) was submitted to MCL sequence clustering [27], with default parameters (e.g., inflation value 2.0); clusters were incrementally assigned to an integer identifier. Clusters are sorted by their size (number of members in a cluster, Supplement S5); thus, the largest clusters have smallest-integer identifiers (see Results and Appendix- Table 1). This approach has also been used successfully elsewhere [28] as a method of choice.
Annotation transfer based on the first chlamydial genome ever sequenced was implemented through the direct matching of the lead sequences to a previously highly curated dataset for Chlamydia trachomatis D/UW-3/CX [7].
The annotation qualifiers used in the manually curated corpus [16] are: ENZYME (for enzymes with EC number assignments), FUNCTION (for other protein functions), SIMILAR-TO (for those sequences with a similarity to a protein of known function but no specific assignment) and DOMAIN (for the existence of a known, named protein sequence domain) [16] (Appendix- Table 1).
Sequence matching of the original dataset to the data collection presented here was performed by MagicMatch [29], which was the first scheme to implement the MD5 checksum for protein sequence identification, an approach later propagated in all major database resources.

Analysis of Unique Genes
All unique genes, i.e., more than 2,000 genes with no similarity within the pangenome, were searched against the non-redundant protein sequence database (nrdb: 15,052,178 entries) [30]. Results from this search were evaluated manually and key similarities were extracted for further investigation (Supplement S6).

Genome Trees
Genome-based trees were calculated using phylogenetic profile distance [31,32]. Similarity values were measured by the shared number of genes represented by phylogenetic profiles, symmetrified by the minimum shared value, normalized by minimum self-similarity and turned into distance values as previously described [32,33] (Supplement S7).

Sequence Alignments
Multiple sequence alignments were performed and visualized by JalView [34]. Novel motifs reported in this work are provided below and in FASTA format (Supplements S8,9).

Data Availability
Per genome contributions to the pangenome are also provided (Supplement 10). All sequence data and results (in 10 Supplements) have been made available at datadryad.org, under the identifier [35].

General Characteristics of the Chlamydiales Pangenome
The Chlamydiales collection herein contains over 40,000 protein-coding genes in total, with ~1,200 genes/genome on average, with significant deviations (Table 1). We take the view to present the two extreme tails of this data collection in detail following the clustering step for the identification of protein families within the pangenome and comment on the intermediate cases. In other words, we primarily focus on the two classes of the most interesting clusters, (i) those containing the core genes and (ii) those corresponding to "unique" genes, without significant similarities within the pangenome, thus singleton clusters. The functional characterization of the entire complement as well as further issues listed in the discussion for future research are clearly beyond the scope of this critical review.

Protein Families
In total, the clustering has yielded 5,554 clusters corresponding to protein families. For practical purposes, we define a protein family as one that contains at least three genes: in that sense, there are 294 cases, which do not detect themselves in this comparison (typically because of either short length, abnormal composition, or both), 2,038 unique genes (singletons) and 1,177 doublets. The remaining 2,045 clusters represent protein families with three or more members, distributed across 37 genomes ( Figure 1).

Figure 1.
Pangenome protein family size distribution. Cluster size is displayed on the x-axis (bins until 50 are all shown; above 50, bins are shown for each ten counts, labels for every five bin sizes); absolute frequency of clusters is shown on the left y-axis (bars, green curve); cumulative count of clusters is shown on the right y-axis (orange curve). Families are defined as those clusters with at least three members (see text); all cluster frequencies are shown here for completeness. The bimodal nature of the distribution can be seen between the peak at low cluster sizes and 37; above 37 there are multi-member and multi-species protein families (see text).
It is evident that the protein family size distribution follows, as expected, the shape of other pangenome analyses, with a clear bimodal distribution, with one peak at low-count families which has been called the "accessory pool" and another peak at the limit of the genomes under consideration, which has been called the "extended core" [21]. The so-called "character genes" (which we prefer to define as "peripheral", as opposed to "core" genes) exhibit, by definition, a heterogeneous distribution across genomes (and between peaks) and present an additional challenge for further interpretation (Figure 1). Genes 2012, 3 The peak at exactly 37 with 312 counts, i.e., 312 families with exactly 37 members, corresponds to the number of 37 genomes analyzed across the pangenome. Beyond that peak, there are 180 protein families with more than 37 members (clusters 1-180) (Supplement S5), of which ten contain more than 100 members and are discussed below.

Figure 2.
Top ten multi-member families within the pangenome. Genomes (with full COGENT-like codes) are shown on the x-axis, sorted by total protein-coding gene count (see also Table 1). Absolute cumulative counts of multi-member families are shown on the y-axis (displayed in the figure legend from left to right and then top to bottom, e.g., ABC transporter permeases, POMPs, type III secretion system ATPases, etc. according to size, see text), color coded according to figure legend.
Following those, there are another four families with more than 110 members each: the EF-Tu/EF-G/LepA family (119 members), the oligopeptide binding protein family OppA (114 members), the GroEL family (111 members) and finally the Ile-Leu-Val (ILV)-tRNA synthetases (111 members). These are followed by two families with more than 100 members, namely the Dihydrolipoamide acetyltransferase E2 component/Dihydrolipoamide succinyltransferase (110 members) and the 3-oxoacyl-[acyl-carrier protein] reductase families (109 members) ( Figure 2). A significant number of multi-member families contain proteins of known function (Supplement S5). Interestingly, families containing only homologues from S. negevensis, W. chondrophila, P. acanthamoebae and Protochlamydia amoebophila are 172 in total, remarkably close to the 171 clusters of "orthologous" proteins in this group of species reported recently [18].

Core Genes
At the other end of the bimodal distribution, there are 312 families with 37 genes each, reflecting the number of genomes analyzed. However, there are eight clusters here with duplicates per genome (clusters 224, 460: S. negevensis; 254, 276, 420: P. acanthamoebae; 255, 272: P. amoebophila; 429: W. chondrophila 203) (two of which of unknown functional roles, Appendix- Table 1). Thus, there are exactly 304 protein families with 37 genes each represented once in each genome, which can be truly called "core" genes, most of which have some source of annotation (Appendix- Table 1). These represent just over a quarter of the average chlamydial genome (304/1182 = 26%).
Annotations transferred from the manually curated seed annotation corpus of C. trachomatis reveal a wide range of functional roles for this core set, as expected (Appendix -Table 1). Indeed, 227 families of the core set can be assigned to a functional role, according to the annotation qualifiers originally used (see Experimental Section). Only an additional 77 cases in this set do not contain any annotation (Appendix- Table 1). It can be argued, therefore, that this level of characterization of 75% (227/304) across 37 genomes signifies a functional coherence that is consistent with our current knowledge of this taxonomic order. This list is provided for further investigation by the community; it is worth pointing out that it encompasses basic cellular roles in genetic information processing (e.g., cluster 184), including transcription (e.g., cluster 187) and translation (e.g., clusters 242-243), metabolic transformations (e.g., cluster 182 or 196), transport systems (e.g., clusters 193-195) and other key processes (e.g., cluster 192). It is interesting to note that apart from complements represented by ribosomal proteins or aminoacyl-tRNA synthetases, other systems are also coherently detected, for example the NifU [39]/NifS [40] genes (clusters 221-222).

Peripheral Genes
In the midst of the two extremes (viz. peaks) of the bimodal family size distribution, there exists a wide variety of cases with an anomalous and clearly heterogeneous pattern. There are 428 families with more than ten and less than 37 members (not shown, available in Supplement S5). Their hererogeneous composition is reflected by the fact that 217 of the 428 families (just over 50%) do not contain a homolog outside the Chlamydiaceae, i.e., across the larger genomes mentioned above. Within this group, however, there is a significant variation of family phylogenetic distribution (not shown) that needs to be explored in future research.

Unique Genes
In total, there are 2,038 unique genes represented by singleton clusters, thus not falling into families within the pangenome. The content of genomes with unique genes varies significantly, from 0 to 796 (S. negevensis), with 55 unique genes on average. In percentage points, this varies from obviously 0 to 32% of the genome (S. negevensis), with an average of just over 3% per genome ( Figure 3). . Correlation between genome size and unique genes. Genome size is given as the number of protein-coding genes (shown on the x-axis) against the count of unique genes (number of unique genes without homologs within the pangenome, shown on the y-axis; y-axis is displayed on logarithmic scale). The six points on the upper right part of the graph are evidently those genomes with largest gene counts, all outside the Chlamydiaceae family (see Table 1 and text). The pattern observed is primarily due to the sampling of taxonomic space of the Chlamydiales and will vary as more genomes from this group become available.
The densest part of the phylogeny exhibits no unique genes-17 genomes, including most of the C. trachomatis and C. psittaci strains, C. pneumoniae CWL029 and C. abortus LLG (Figure 3, missing points corresponding to 17 genomes with zero value on the y-coordinate, available in Supplement S5). Twenty genomes have unique genes, of which six genomes have less than 10 such genes and one with 15 unique genes (Figure 3), all from the above group, or less than 2% of their genome entries. Another five genomes with a handful of unique genes are C. pneumoniae AR39 (33/3%), TW-183 (43/4%) and LPCoLN (60/5%) as well as C. felis (27/3%) and C. caviae (29/3%). The remaining eight genomes contain the majority of unique genes, 1818 in number or 89% of total, ranging from 66 (W. chondrophila WSU, 3% of genome) to 796 genes (S. negevensis, 32% of genome). This is not entirely a biological effect, rather a sampling artifact arising from the deeper sequencing of the C. trachomatis/C. pneumoniae group (see below).
The six outliers which form a different group above (upper right, Figure 3) are all species with large genomes (ca. 2,000 protein-coding genes or more): the two W. chondrophila strains (3-4%), the two P. acanthamoebae strains (4-8%), P. amoebophila (20%), and S. negevensis (32%), listed here according to the absolute number of their unique genes per genome. In relative terms, however, two species namely C. muridarum (36/4%), and C. pecorum (76/8%) contain a significant number of unique genes given their relatively small genome size (both less than 1,000 protein-coding genes).

Properties of Unique Genes
The genes considered as singletons in this analysis are 2,038 as mentioned above. Of those, a number of short genes might fall into pangenome families (not shown) but do not seriously affect the overall assessment (e.g., case CCAV-GPI-01-000824 in Supplement S6). This is an artifact of sensitivity for the two different searches, first against the 40,000 or so genes of the pangenome and second against the entire nrdb database of more than 15 million sequences. While a full analysis of the unique gene complement of the Chlamydiales is under progress, it is interesting to report on a number of findings pertinent to this work.
One such domain is an enigmatic, short and highly conserved motif containing the triplet Pro-Cys-Tyr (PCY), present in the C. pneumoniae AR39 CP0988 protein. This protein is 52 residues long and does not exhibit significant similarities to any other protein in the Chlamydiales pangenome. However, it does show similarity to a set of short proteins (<100 residues long) from various species, including Acinetobacter, Brucella, Clostridium, Coxiella, Curvibacter, Eubacterium, Parvimonas, Rhizobium, Ruminococcus, Selenomonas, Streptomyces, other longer proteins from Chloroflexi, Heliobacterium, Lactobacillus, the C-terminus of a Propionibacterium protein (HL046PA2) and an uncultured Acidobacteria bacterium HF4000_26D02, and importantly, to a number of longer plant proteins from Nicotiana tabacum, Pinus koraiensis, Solanum demissum (middle of protein) and Vitis vinifera (N-terminus, total length 1,193 residues) (Supplement S8). This conserved region with this peculiar phylogenetic distribution has not been characterized previously to our knowledge, and can be considered a genuine novel domain of unknown function (Figure 4). It remains unclear whether the domain has been universally lost from the Chlamydiales pangenome or acquired from C. pneumoniae through horizontal transfer.
Another interesting example of a unique protein is the P. amoebophila pc0506. This 82-residue-long uncharacterized protein is evidently absent from the core pangenome and yet it exhibits significant similarity to four Verrucomicrobia proteins from Verrucomicrobium spinosum, Chthoniobacter flavus, Pedosphaera parvula and Coraliomargarita akajimensis, in this order of similarity, ranging from 53% down to 44% sequence identity ( Figure 5). The above mentioned proteins reportedly belong the leucyl aminopeptidase superfamily (Supplement S9). The functional significance of this biochemical role for P. amoebophila is not yet understood. Yet, the strong mutual similarity of this protein family with Verrucomicrobial and P. amoebophila members (no other member in the entire pangenome) can be placed within the general controversy of the connection of Chlamydiales with the so-called PVC group [41,42] (see below). The domain was discovered following five iterations with PSI-BLAST with CP0988 as query sequence (GI:16752158), until convergence and an e-value threshold 0.005. In total 70 sequences were recovered; redundancy was removed at 95% with Jalview [34], resulting in 32 sequences shown here. The length of the domain is just 30 residues; boxes signify sequence identity at 50% or above (darker color: more conserved). GI labels are provided, along with sequence coordinates on the left of the alignment (see text for more details and discussion). In all, it appears that properties encoded from most unique genes, apart from their unusual phylogenetic distribution, represent accessory functional roles that provide additional versatility to the largest genomes in the group, possibly related to their extra functional capabilities. Two exceptions with seemingly central functions are wcw_0805, with similarity to the 50S L34 ribosomal protein family and wcw_861, with similarity to 6-pyruvoyl tetrahydrobiopterin synthases, both from W. chondrophila WSU (not shown).

Protein Family Contributions from Genome Projects
As mentioned above, we have tracked the original publication (and/or release) data for the genomes under consideration, in terms of novel families detected per genome sequence (Supplement S10). By mapping the protein families which appear first in this ranking order, we can thus estimate the relative "novelty" or contribution of previously unseen protein families within the chlamydial pangenome and the typical "pangenome saturation curve" (Figure 6). Figure 6. Protein family contributions from genome projects. Genome codes are sorted according to their original publication date (and/or release date, x-axis); absolute number of "novel" protein families within the pangenome are given (left y-axis, blue curve and square symbols); cumulative sum of protein families (up to 5,260, excluding those without self-hits, see text) is also shown, defined as a "pangenome saturation curve" (right y-axis, green curve and square symbols).
As expected, and discussed above (Figure 3), for the densest part of the group, little or no contributions have been provided. Apart from the larger genomes, which have added hundreds of new gene types [19], the more distant members of the group with small genomes, for instance C. caviae or C. pecorum, have also contributed a significant number (80 and 76, respectively-Supplement S10).

Genome Phylogeny
Finally, we have reconstructed the genome phylogeny of the pangenome based on the sharing of phylogenetic profile patterns based on the above analysis (see Experimental Section). Evidently, the pangenome is stratified according to the known, established phylogeny patterns [10] (Figure 7). The genome tree is another concise way to visualize the "novelty" components of the various species and strains that have been sequenced, exemplified above in various contexts, e.g., number of unique genes (Figure 3) or the tracking of the relative contributions of novel protein families ( Figure 6). A future aspect of this work will be to infer the history of the pangenome using methods of ancestral state reconstruction [43]. The evolutionary history of the Chlamydiales as reflected by the genome tree might also shed light on the ongoing controversy about their status within the tree of life [41]. The genome tree accurately reflects the current taxonomy of Chlamydiales [4,44], with a couple of notable exceptions namely the clustering of C. abortus with C. psittaci, the closer relationship of C. felis with the former two species against C. caviae-in agreement with previous findings [6,44] but not with other proposals [8]-as well as the distinct relationship of C. pecorum at the root of Chlamydiacae and not as a sister group of C. pneumoniae [4,44]. The resulting phylogenetic tree using genome-wide phylogenetic profile sharing patterns can also act as an internal control of the pangenome analysis, since all the closely related strains sequenced are grouped together with very high accuracy (Figure 7).

Discussion
Our results suggest that the Chlamydiales pangenome reflects a certain degree of structural stability, as core genes represent over a quarter of an average genome, as well as functional coherence, in the sense that most functional properties of these genes are consistent with current knowledge. Unlike various claims in the recent literature, it turns out that, at least in the case of a highly constrained pangenome of intracellular pathogens, there is an unexpected degree of stability, given the wide range of phylogenetic relationships within this particular taxon.
It is thus shown that for the smallest of genomes (<900 protein-coding genes), over a third of their gene content is shared with larger genomes (>2,000 genes), decorated by a broader element of so-called "character", or peripheral, genes. This distribution, which in turn is influenced by the sampling of phylogeny and other factors, requires further investigation, being beyond the scope of this work.
It should also be pointed out that the Chlamydiales pangenome exhibits general characteristics of distribution not dissimilar to other recent pangenome analyses, including those of the Salmonella pangenome with 45 strains [15], the Streptococcus pneumoniae pangenome with 44 strains [14] and the Campylobacter pangenome with 96 strains [28], suggesting the conservation of a core pangenome within and across bacterial taxa that have been sampled adequately. In the case of Salmonella, tracking the contributions of new strains to the entire core set and the pangenome suggests a slight expansion with more sampling and a stable core, reminiscent of the Chlamydiales, with one third of the pangenome represented in the core set [15]. A slightly less stable pattern is detected in the Streptococcus pneumoniae group [14], possibly due to a wider diversity in that sample, yet with a similar pattern of core set saturation. Interestingly, an attempt for ancestral reconstruction in the S. pneumoniae/S. mitis complex suggests that there is a dual process of genome expansion and reduction in the different paths leading to the genomes of contemporary strains [14]. A more comprehensive analysis of the Campylobacter pangenome with 96 strains [28], using a combination of experimental and theoretical work, also points to the same direction: Within the two species groups examined, the core gene set overlap reaches 80%, supporting earlier findings for the related Helicobacter pylori strains [45].

Conclusions
We have thus examined the salient features of the Chlamydiales pangenome, introducing a pangenome analysis pipeline and certain definitions that facilitate the discovery of core and peripheral genes, the identification of unique genes with various origins as well as the detection of novel protein sequence domains. We expect that analogous efforts will lead to rigorous standards for pangenome analysis in the future. Future research opportunities abound, for example: ancestral reconstruction [43], syntenic patterns of genome structure (e.g., [28,45]), the (presently limited) enrichment with expression data, the evolutionary histories of 'peripheral' genes (as discussed above), the connection of Chlamydiales with plants [46][47][48][49][50], the position of the Chlamydiales in the tree of life, and the connection with the PVC superphylum [41,42,50]. Wider challenges that go beyond the above pangenome-specific issues might include a more detailed annotation of the entire dynamic range of family distribution [21], the characterization of protein function in a wider context including comparative metabolic reconstructions [19], the evolution of mobile elements [51], the deeper understanding of the physiological and pathological properties [52,53] of the strains that have been sequenced and the connection with other pangenomes [28]. Genes 2012, 3