Genome assemblies were downloaded from the RefSeq database, and from Genbank for the assemblies absent from RefSeq (Table 1; O'Leary et al., 2016). All genome assemblies from the phylum Chlamydiae (NCBI taxonomy ID 204428) were included in this study with the exception of the Chlamydia genus. Since several Chlamydia species are well-studied pathogens with many sequenced genomes, a single representative genome assembly was considered for each of the 13 Chlamydia species out of a total of 326 genome assemblies (June 2017).

Many genomes included in this work are metagenomics bins that may be of unequal quality. Indeed, metagenomics bins can be an incomplete representation of the organism's genomes, can be a mixture of multiple genomes or include chimeric DNA sequences. To screen for contaminants, the completeness and redundancy of each genome was thus evaluated based on the identification of 104 nearly universal (and generally single copy) bacterial markers with checkM (Parks et al., 2015). In addition, translated predicted coding sequences (CDS) of each genome were compared to RefSeq database using PLAST, a fast sequence similarity search tool, with the following parameters: -M BLOSUM62 -s 45 -seeds-use-ratio 60 -G 11 -E 1 -F F -max-hit-per-query 100 -max-hsp-per-hit 1 (Van Nguyen and Lavenier, 2009). PLAST was favored over BLAST for its lower memory usage. The taxonomy of the first PLAST hit (excluding hits against RefSeq sequences classified as belonging to the same species) was investigated with help of the NCBI taxonomy database (Federhen, 2012).

Orthologs were identified using OrthoFinder version 0.4 (Emms and Kelly, 2015). A reference phylogeny was built based on single copy orthologs conserved in at least 55 of the 60 studied genomes (Table 1). Core single copy orthologs were aligned using mafft version 7.058b (Katoh and Standley, 2013). The concatenated alignment was used for the reconstruction of the species tree using FastTree 2.1.9 with double precision (Price et al., 2010). Circular plots were drawn using Circos version 0.69 (Krzywinski et al., 2009), whereas other plots were made with R (R Core Team, 2016).

All genomes were annotated using GhostKOALA (Kanehisa et al., 2016) and Interproscan version 5.23–62.0 (Mitchell et al., 2015). Orthologs of the type III and type IV secretion systems, known effectors, division proteins, membrane proteins, the Euo master regulator, respiratory chain complexes and enzymes involved in menaquinone biosynthesis and glycogen metabolism were identified based on OrthoFinder grouping into groups of orthologous proteins (see Table S1 for the detailed reference locus list). The ATP:ADP antiporters homologs were identified by InterProScan annotation (interpro accession: IPR004667). The flagellar apparatus subunits, proteins involved in peptidoglycan, purine, pyrimidine and ubiquinone biosynthesis as well as the glycolysis pathway, the citrate cycle (TCA), and the pentose phosphate pathway (PPP) were identified based on KEGG ortholog (KO) annotation. Reported counts for each pathway/module are the non-redundant number of identified KO. Predicted coding sequences were compared to COG database (version 2014) (Galperin et al., 2015) using BLASTP version 2.3.0+(Camacho et al., 2009) with an e-value cutoff of 1e⁻⁵, a minimal query coverage of 50% and a minimal identity of 20%.

Proteins conserved between pairs of genomes were identified by pairwise protein sequences comparisons with BLASTP version 2.3.0+ (Camacho et al., 2009). The identity of all reciprocal best blast hit (RBBH) with a minimal e-value of 10⁻⁵ and 50% of both query and hit coverage were retained for the comparative analyses. 16S and 23S rRNA sequences were identified using barrnap (https://github.com/tseemann/barrnap). HMM profiles were built from 21 reference sequences for nine previously identified phylogenetically informative markers (Pillonel et al., 2015). These profiles were used to identify markers in all chlamydial genomes included in this study using hmmsearch v3.1 (Eddy, 2011). For each marker, a phylogenetic tree was reconstructed based on the best hmmsearch hit in each genome. Bitscore cut-offs were defined individually for each marker based on the density distribution of hmmsearch bit scores and gene tree topologies: DnaA (260), SucA (750), the hypothetical protein 325 (350), FabI (400), RpoN (300), FtsK (750), PepF (470), Adk (175), and HemL (360). Pairwise amino acid sequence identities were calculated based on pairwise Needleman-Wunsch global alignments made with Needle from the EMBOSS package version 6.6 (Rice et al., 2000). Gaps were not considered in identity calculations. The probability distributions of pairwise amino acid identity values were estimated using the statistical package R with kernel density estimations (R Core Team, 2016). Density plots were made using the ggplot2 package (Wickham, 2016).

In total, 60 genomes (including one outgroup) were compared in this study: 21 poorly characterized metagenomics bins classified as Chlamydiae were retrieved from Genbank. Most (n = 17) originate from a single study of groundwater microbial ecosystems (Anantharaman et al., 2016). Others were sequenced as part of the microbial community of a biologically active filter of a water treatment plant (n = 1) (Pinto et al., 2016), as part of an investigation of bacterial communities in estuary sediments (n = 1) (Baker et al., 2015), and from an experimental bioreactor used to treat contaminated goldmine water (n = 2) (Kantor et al., 2015). In addition, three single-cell amplified partial genomes from marine Chlamydiae (Collingro et al., 2017) and 25 previously published chlamydial genomes (Table 1) were studied. The latter contained a majority (n = 20) of draft assemblies. The verrucomicrobial genome of Akkermansia muciniphila was included as an outgroup.

While most Chlamydiae share a similar life cycle, their genetic diversity is considerable. Standard procedures for the classification of Chlamydiae rely on a small set of reference genes or protein sequences, including the 16S rRNA sequence (Pillonel et al., 2015; Greub, 2017). Given that metagenomic datasets frequently lack rRNA operons, pairwise comparisons of all protein sequences were undertaken to better characterize the genetic diversity of newly sequenced genomes. The protein sequences of the 59 chlamydial genomes included in this comparative analysis were clustered into 10,162 orthologous groups (Table S2, Figure S1A). Among those groups, 4,553 are restricted to a single genome (44.8%, Figure S1B). Only three single copy orthologs were conserved in all chlamydial genomes (the two ribosomal proteins S8 and L6, and the tRNA threonylcarbamoyladenosine biosynthesis protein TsaE), but 108 orthologous groups were conserved in more than 90% of the genomes (55 out of 60, including the outgroup genome.The list of 108 orthologous groups is reported in Table S3). The concatenated alignment of those 108 protein sequences was used to build a reference phylogeny of all representative genomes of the phylum Chlamydiae (Figure 1A).

The GC content of uncultured strains varies widely from 26.2 to 49.3% (Table 1, Figure 1B). Twelve of the metagenomics bins are predicted to be more than ninety percent complete, based on the identification of 104 nearly universal bacterial proteins (Table 1, Table S4, Figure 1B). Those nearly complete assemblies range from 1.42 to 2.83 Mb in size. The largest metagenomics bin (Chlamydia sp. 32–24), is phylogenetically related to Neochlamydia spp. The smallest one, Chlamydiae bacterium RIFCSPLOWO2_01_FULL_28_7, only comprises about 700 kb that were predicted to represent approximately 58% of the genome. The most incomplete genomes are the three genomes from marine Chlamydiae (Figure 1, groups 2, 6, and 7). Those genomes are amplified from a single cell. Such approach frequently yield incomplete genomes (Collingro et al., 2017). Chlamydiales bacterium SCGC AG-110-M15 missed 55 of the 108 core proteins used to build the species phylogeny (50.9% of the dataset, Figure 1C).

Despite the metagenomics approach to sequence Chlamydiae from more complex samples, only few genomes exhibited signs of contamination. The strongest evidence was present in the S. epinepheli genome where 32 of the 104 nearly universal proteins were present in more than one copy, indicating that it might be an admixture of two genomes (Figure 1C). Nevertheless, the 32 duplicated S. epinepheli markers were all monophyletic in phylogenetic reconstructions including all chlamydial homologs and their closest homologs in the RefSeq database, indicating that the assembly might be an admixture of two closely related strains. In order to further evaluate if the assemblies contained only fragments of Chlamydiae genomes, all CDS were compared to the RefSeq database. A majority of CDSs of each analyzed genome exhibited a best non-identity PLAST hit (excluding hits against RefSeq sequences classified as belonging to the same species) against chlamydial sequences present in the RefSeq database (Figure 1C). Assemblies reconstructed from metagenomics and single cell data exhibit between 55.77 and 75.93% of best PLAST hits against chlamydial sequences, which is higher than for S. negevensis (45.01%). S. negevensis is the only representative of the family Simkaniaceae currently available in RefSeq. The higher proportion of PLAST hits against non-chlamydial sequences (32.48%) and without any hits (22.51%) as compared to representative of other new families absent from RefSeq might be related to the much larger size of S. negevensis genome (Figure 1B). Several of the most fragmented assemblies harbored contigs without any best hits against chlamydial sequences (Figure 1C). Nevertheless, chlamydial hits were systematically distributed along the whole assembly, with only few small contigs exhibiting no best hits against chlamydial sequences, as can be seen on Figure 2 for the assembly RIFCSPHIGHO2_01_FULL_44_39. Overall, these results indicate that all assemblies contain mostly chlamydial DNA sequences, supporting the reliability of the analysis.

We recently proposed a scheme for the classification of chlamydial genomes at various taxonomic levels (species, genus and family). This scheme was used here to update the classification of all unclassified chlamydial genomes (Figure 3). Briefly, the strain HS-T3 likely belongs to a new genus in the family Parachlamydiaceae (Figure S2). Chlamydia sp. 38–26 is a new Neochlamydia species (Figure S3) and Chlamydia sp. 32–24 is a Rubidus massiliensis strain (Figure S4). There are four Pr. amoebophila strains (Figure S5) and three strains of the same Neochlamydia species (S13, TUME1 and EPS4, Figure S6).

Regarding the new deeply branching assemblies, Group 3 and Group 4 encompass each two candidate species from two different genera (Figures S7, S8). The seven assemblies of group 5 are part of the same candidate species (Figure S9). Altogether, the new genomes can be classified into 12 new candidate species (Figure 3). Half of them are only represented by <90% complete genomic data (Figure 1B, Table 1).

Most genome assemblies sequenced using a metagenomics approach lack rRNA operons: 15 lack 16S rRNA, whereas 17 lack 23S rRNA (Figure S2). This happens frequently as rRNA operons are too highly conserved to be assigned reliably to specific genome assemblies. In addition, most genomes lack one or several of the nine proposed taxonomic markers (Figures S2–S9; Pillonel et al., 2015). This is a major limitation for any scheme relying on a limited set of genes. Therefore, in order to evaluate the genus and family level diversity of the Chlamydiae phylum, the identity of reciprocal best blast hits (RBBH) was used to evaluate the divergence of chlamydial lineages at various taxonomic depths.

Figure 4 shows the distribution of RBBH identity values between pairs of genomes. As expected, the distribution of sequence identity values shifts continuously toward lower values with increased phylogenetic distance (Figure 4A). W. chondrophila and the different genera classified in the Parachlamydiaceae family show similar levels of divergence (Figure 4B). This is not the case with S. negevensis, exhibiting a median identity lower than 50% with P. amoebophila (Figure 4B). The median of the median pairwise identities among members of the Chlamydiales and Parachlamydiales (Figure 3) is of 44.44%, whereas Chlamydiaceae and the two deep branching clades 1 and 2 exhibit a median identity of respectively, 39.81 and 41.14% (using one representative genome per species, see detailed values in Figure S10). The high divergence of clade 1 and 2 is even more obvious on Figure 4A, with a clear shift of the distribution of the RBBH identities peaking at about 35% identity between C. trachomatis and clade 1 and 2. Altogether, the phylum exhibits at least three highly divergent clades, and may now include seven new candidate family-level lineages and 11 new genera (Figure 3, Figure S11).

Gupta and colleagues recently disputed the 16S rRNA identity cutoffs proposed to delineate families of the phylum Chlamydiae, arguing that the number of chlamydial families was inflated by the high cutoff in use (Gupta et al., 2015). Nevertheless, different bacterial clades may exhibit variable rates of sequences evolution (Kuo and Ochman, 2009). In addition, the 16S rRNA gene may not evolve at a constant rate across the entire tree of life and may not necessarily be a reliable indicator of whole genome relatedness (Konstantinidis and Tiedje, 2005; Kuo and Ochman, 2009). In order to put the Chlamydiae phylum into a broader perspective, we compared the 16S rRNA pairwise sequence identity and the conservation of protein sequences with data from another diverse group of intracellular bacteria, the order Rickettsiales. Representatives of the Rhizobiales order were also included for comparison (Table S5).

Figure 5A shows the relationships of 16S rRNA and median RBBH identities for the three considered orders. For both measures, the Chlamydiales and Rickettsiales orders exhibit higher sequence divergence than the Rhizobiales order. The sequence divergence of different families belonging to the same order (see Tables S5, S6 for the detailed classification used here) reveals that the 16S rRNA gene of Rickettsiales families is more divergent than in the Chlamydiales families (Figure 5B). On the other hand, the median pairwise protein identity is also relatively low for most Chlamydiales (Figures S12A–D). Given the skewed shape of RBBH identity distributions (Figure 4), the median identity may not be the most appropriate summary metrics to estimate genome divergence. An alternative would be the maximum of the density distribution of RBBH identity rather than the median (Figures 5C,D). Using the maximum of the estimated distribution of amino acid identities, the Chlamydiae cloud shifted toward lower values, indicating that the distribution of RBBH identities in pairs of chlamydial genomes is more skewed than for Rickettsiales and Rhizobiales. This still holds true when comparing Parachlamydiales only (excluding Chlamydiaceae and Similichlamydiaceae).

The difference between Rickettsiales and Chlamydiales could be due to the fact that the Chlamydiales dataset includes incomplete genomes. Missing data might lead to spurious reciprocal best blast hits, leading to an undervaluation of whole genome relatedness. Nevertheless, the trend is still visible when only complete chlamydial genomes are considered (Figures S12E,F). In conclusion, chlamydial families exhibit higher conservation of 16S rRNA but a skewed distribution of whole genome RBBH identity toward values that are as low, or even lower than genomes exhibiting lower 16S rRNA identities.

Several genomes included in this analysis are only a partial representation of the true chlamydial genome, which precludes any in-depth genome content comparisons. Nevertheless, the identification of even partial molecular machineries and metabolic pathways can provide first insights into the ancestrality and conservation of genomic features that may be essential to their unique lifestyle. The recent sequencing of S. epinepheli, the most diverging Chlamydiae described thus far, revealed the presence of a type III secretion system (T3SS) and the conservation of several key virulence factors (Taylor-Brown et al., 2017). T3SS components were identified in all new chlamydial clades (Figure 6, column T3SS), including the most incomplete ones. Several described T3SS effectors such as Mip and Nue are also widely distributed within the whole phylum (Figure S13). The entire clade that includes the Simkaniaceae, the Rhabdochlamydiaceae and several new lineages lack any homolog of Type II secretion system effector protease CPAF (Figure S13). In C. trachomatis, CPAF mutants display impaired generation of infectious elementary bodies (Snavely et al., 2014) and cpaf is therefore an essential gene. Our observations suggest that CPAF is dispensable in some chlamydial lineages, and its role in other Chlamydia-related organisms where it is conserved remains to be confirmed.

As obligate intracellular bacteria, Chlamydiae acquire essential nutrients from the cytosol of their host, including nucleotides and sugars such as ATP and glucose-6-phosphate. All but one partial genome encode at least two ADP/ATP antiporter homologs. These transporters are involved in nucleotide, ATP or nicotinamide adenine dinucleotide uptake from the host (Tjaden et al., 1999; Greub and Raoult, 2003; Haferkamp et al., 2004; Fisher et al., 2013). Furthermore, most genomes also encode a homolog of the glucose-6-phosphate transporter UhpC.

The genomes of Simkania negevensis and several Parachlamydiaceae strains encode an apparently complete type IV secretion system (T4SS) located either on the chromosome or on a plasmid (Bertelli et al., 2016) that may be involved in conjugative DNA transfer (Greub et al., 2004). Several strains of Parachlamydia acanthamoebae encode remnants of the T4SS in their genome (Greub et al., 2009; Bertelli et al., 2016). No traces of T4SS genes could be found in any other chlamydial lineages (Figure 6). The T4SS of S. negevensis, Protochlamydia naegleriophila KNic and R. massiliensis are located on a plasmid, which may have facilitated its loss or acquisition in the various chlamydial lineages. Assignment of plasmid sequences to specific organisms within the context of a metagenomic dataset is extremely difficult. Chlamydial plasmids may thus be missing in some assemblies, including for the genome of Chlamydia sp. 32–34.

It was recently suggested that some marine Chlamydiae may possess a flagellum regulated by a chemotaxis system (Collingro et al., 2017). The identification of several flagellar and chemotaxis genes in highly divergent chlamydial genomes retrieved by single cell-genomics was interpreted as an evidence that motility was an ancient feature of the Chlamydiae phylum (Collingro et al., 2017). Although most new uncultivated chlamydial clades were also sequenced from aquatic samples, none exhibited homologs of the flagellar genes identified in the genome assemblies SCGC AG-110-P3, AB-751-O23 and AG-110-M15 (Figure 6, group 2, 6, and 7), questioning whether this flagellar system is as old as Collingro et al. hypothesized. Two proteins of the type III secretion system that are homologous to flagellar proteins, fliQ/sctS (K02420) and fliN/sctQ (K02417), were sometimes annotated as part of the flagellum (Figure 6) but did not branch with SCGC Chlamydiae homologs in phylogenetic reconstructions (see for example Figure S14). A chemotaxis system is also encoded in several chlamydial genomes lacking flagellar genes including several Parachlamydiaceae (P. phocaeensis, P. naegleriophila, and P. acanthamoebae) and the two Criblamydiaceae genomes currently available. This system was thus proposed to be involved in signal transduction (Collingro et al., 2011). In the present analysis, chemotaxis genes were identified in the genome of Chlamydia sp. 38–26, which is closely related to Neochlamydia sp. strains that do not harbor a chemotaxis system. Overall, our data suggests that all Chlamydiae exhibit a similar obligate intracellular lifestyle, and that most are non-motile.

All described Chlamydiae exhibit a similar biphasic developmental cycle, including the most deeply branching clades (Abdelrahman and Belland, 2005; Steigen et al., 2013; Seth-Smith et al., 2017). The transcription factor Euo targets more than 100 promoters in W. chondrophila, and may be a key player in cell cycle regulation (Domman and Horn, 2015; de Barsy et al., 2016). It is conserved throughout the whole phylum, including S. epinepheli and “Ca. Novochlamydiaceae” (Figure 6), with exception of three of the most incomplete genomes probably due to the partial data available. Interestingly, the genome assembly of S. epinepheli lacks all proteins involved in peptidoglycan biosynthesis (Figure 6) and all but three proteins involved in chlamydial division (Figure S13). Remnants of those proteins could not be found in the raw assembly using tBLASTn. Although the genome of this strain is partially covered, peptidoglycan synthesis genes are located in various regions along the chromosome, which reinforces the likely absence of these genes from S. epinepheli. Indeed, this bacterium divides through a budding process (Seth-Smith et al., 2017), which suggests that the family Parilichlamydiaceae may present a significantly different division mechanism than other chlamydial lineages.

Menaquinones are electron shuttles that allow the transfer of electrons between membrane-bound protein complexes in the electron transport chain (Meganathan, 2001). The Chlamydiaceae synthesize menaquinone through the alternative futalosine pathway (Hiratsuka et al., 2008; Barta et al., 2014). Other clades, including the deep branching “Ca. Novochlamydiaceae” (Figure 6), encode the traditional pathway, supporting the hypothesis that the futalosine pathway was horizontally acquired by the Chlamydiaceae lineage (Barta et al., 2014; Subtil et al., 2014). No menaquinone biosynthesis enzymes could be identified in the S. epinepheli genome (Figure 6). A homolog of MenB was identified (Figure 6) but it exhibits more similarity with enoyl-CoA hydratase and did not cluster with other chlamydial sequences in phylogenetic reconstructions, thus likely having another enzymatic role than MenB (Figure S15). Most genomes encode the necessary components of a minimal respiratory chain except the seven genomes of “Ca. Arenachlamydiaceae” (Figure 6). Those 1.4 Mb genomes are predicted to be nearly complete but do not encode the menaquinone biosynthesis pathway, nor any component of the respiratory chain, suggesting that they use alternative respiration systems.

Overall, enzymes of the respiratory chain exhibit variable patterns of parallel losses across the entire phylum. The NADH-quinone oxidoreductase (complex I) is for instance restricted to the Parachlamydiaceae-Waddliaceae-Criblamydiaceae clade, and was lost at least three times independently in subclades (Figure 6). Few proteins homologous to complex I were identified in the low GC genomes recovered from estuary sediments and groundwater (“Ca. Limichlamydiaceae”). They encode a putative NADP-reducing hydrogenase exhibiting 44.5–56% amino acid identity to the hnd operon of Desulfovibrio fructosivorans (Malki et al., 1995; Figure S16). The F-type ATPase is also restricted to Parachlamydiaceae-Waddliaceae-Criblamydiaceae clade, as opposed to the V-type ATPase that is conserved in all sequenced chlamydial genomes.

The citrate cycle (TCA) is incomplete in all Chlamydiaceae due to the absence of three enzymes; citrate synthase (GltA), aconitase (Acn) and isocitrate dehydrogenase (Icd) (Omsland et al., 2014). Several Chlamydiae such as the three Neochlamydia sp. and the seven “Ca. Arenachlamydiaceae” genomes (repeatedly missing the same set of enzymes) exhibit an even less complete TCA (Figure 6). Indeed, no homologs of TCA enzymes could be identified in the S. epinepheli genome assembly. On the other hand, several homologs of the pentose phosphate pathway (PPP), glycolysis and glycogen metabolism were identified in all assemblies that are more than 90% complete (Figure 6).

The various Parachlamydiaceae species exhibit significant differences in their ability to synthetize amino acids, with reduced biosynthetic abilities in Protochlamydia amoebophila strains and Neochlamydia strains and extended metabolic abilities for Parachlamydia strains (Figure S13). The 3 closely related Neochlamydia genomes exhibit genomes larger than 2.5 Mbp, but overall highly reduced metabolic capacities (Figures S17, S18). On the other hand, they encode large numbers of mobile genetic elements and repeat proteins (Figure S17). S. epinepheli has the most reduced metabolic capacities, but its genome is only about 70% complete. While the metabolism of amino acids, nucleotides, cofactors and vitamins was not investigated in detail, we generally observed that all genomes had limited predicted biosynthetic capabilities for all those compounds and that marked differences exist within the phylum (Figure S13), as already observed in previous comparative analyses (Bertelli et al., 2010, 2015; Collingro et al., 2011; Omsland et al., 2014).

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.