Chlamydia trachomatis In Vivo to In Vitro Transition Reveals Mechanisms of Phase Variation and Down-Regulation of Virulence Factors

Abstract

Research on the obligate intracellular bacterium Chlamydia trachomatis demands culture in cell-lines, but the adaptive process behind the in vivo to in vitro transition is not understood. We assessed the genomic and transcriptomic dynamics underlying C. trachomatis in vitro adaptation of strains representing the three disease groups (ocular, epithelial-genital and lymphogranuloma venereum) propagated in epithelial cells over multiple passages. We found genetic features potentially underlying phase variation mechanisms mediating the regulation of a lipid A biosynthesis enzyme (CT533/LpxC), and the functionality of the cytotoxin (CT166) through an ON/OFF mechanism. We detected inactivating mutations in CT713/porB, a scenario suggesting metabolic adaptation to the available carbon source. CT135 was inactivated in a tropism-specific manner, with CT135-negative clones emerging for all epithelial-genital populations (but not for LGV and ocular populations) and rapidly increasing in frequency (~23% mutants per 10 passages). RNA-sequencing analyses revealed that a deletion event involving CT135 impacted the expression of multiple virulence factors, namely effectors known to play a role in the C. trachomatis host-cell invasion or subversion (e.g., CT456/Tarp, CT694, CT875/TepP and CT868/ChlaDub1). This reflects a scenario of attenuation of C. trachomatis virulence in vitro, which may take place independently or in a cumulative fashion with the also observed down-regulation of plasmid-related virulence factors. This issue may be relevant on behalf of the recent advances in Chlamydia mutagenesis and transformation where culture propagation for selecting mutants/transformants is mandatory. Finally, there was an increase in the growth rate for all strains, reflecting gradual fitness enhancement over time. In general, these data shed light on the adaptive process underlying the C. trachomatis in vivo to in vitro transition, and indicates that it would be prudent to restrict culture propagation to minimal passages and check the status of the CT135 genotype in order to avoid the selection of CT135-negative mutants, likely originating less virulent strains.

Fig 1. Phase variation mediated by variable homopolymeric tracts.

Panel A. The graph shows the evolution throughout passaging of the percentage of sequence reads with different ‘A’ counts in the homopolymeric tract upstream from CT533/lpxC for the strain E/CS1025/11. The poly(A) tract corresponds to poly(T) in the annotated leading strand. Panel B. Schematic view of the putative promoter region of CT533/lpxC. The predicted transcription start site [126] is labeled by +1. The variable poly(A) tract (in bold) falls between the predicted -35 and -10 hexamers (underlined). BLAST analyses revealed the existence of variable number of ‘A’ counts in C. trachomatis genomes, and also that the nucleotide indicated with an arrow is deleted exclusively in all LGV strains. Panel C. The graph shows the percentage of sequence reads with different ‘G’ counts in the variable homopolymeric tract of CT166 found in the initial populations of the epithelial-genital strains. “G” counts of nine correspond to an “ON” protein. Panel D. Schematic view of the four positions (numbers 1 to 4) relative to a gene at which contingency loci (e.g., homopolymeric tracts) can cause phase variation (adapted from van der Woude MW and Bäumler AJ, Clin Microbiol Rev 17:581–611, 2004 [152]). Whereas positions 1 and 2 are associated with transcription initiation and position 4 with translation (ON/OFF), the mechanism regarding the position 3 is not completely disclosed. We found heterogeneity in length within homopolymeric tracts located in positions 2 (for CT533/lpxC) (blue), 3 (for CT043/slc1 and the operon CT134-CT135), and 4 (for CT166—cytotoxin) (red).

Fig 2. Mutational scenario throughout experimental evolution.

Panel A. Chromosomal location of the genomic alterations observed during the in vitro passaging. The chromosomal position of each mutation (scale adjusted and given by the locus name) and the type of mutation event (inactivating events represented in red) are shown for each strain (see also S3 Table for details). Inactivating SNPs or indels refer to events leading to protein truncation (regardless the length of the resulting protein). For the strain D/CS637/11, the CT135 inactivating event involved the entire gene deletion between direct repeats (Fig 3). Panel B. Dynamics of the emergence and spread of mutations and their frequency in the evolving bacterial populations. For each time-point (passages 5–7, 10, 20, 30, 50 and 100), circular graphs show the frequency of the mutations in the bacterial population, where each color represents a different mutated locus. The number of bacterial generations was estimated taking into account the minimum and maximum values of the mean doubling time of the strains analyzed at each time-point, and assuming a conservative approach by considering 15 hours of exponential phase per bacterial life-cycle (i.e, per passage). Loci designations are based on genome annotation of the D/UW3 strain (GenBank accession number NC_000117).

Fig 3. Schematic representation of the CT135 deletion in the serovar D strain.

The inactivating event of CT135 involved the complete gene deletion between direct repeats (in blue) and the putative formation of a fusion gene enrolling the two flanking genes (CT134 and CT136). The underlying mechanism likely relied in one of three major pathways: A—intermolecular crossing over between direct repeats followed by recombination (yielding both a tandem duplication and a deletion); B—looping out in between direct repeats followed by recombination; and C—DNA polymerase slippage during DNA replication [–91]. The figure also shows the position of all CT135 frameshift mutations (labeled by Ψ) reported here and elsewhere [46, 48, 131, 149, 150, 164], demonstrating that the strains evolved towards CT135 inactivation regardless the “genetic pathway” that drove that inactivation The bilobal hydrophobic domains that putatively enable the insertion of the CT134 and CT135 proteins into the inclusion membrane [33] are highlighted in grey.

Fig 4. Impact of in vitro passaging on the C. trachomatis growth kinetics.

Panels A and B. Comparison of the growth rates and doubling times between ancestral (grey) and evolved populations (black). The percentage values above the bars correspond to the growth rate increment of the evolved population relatively to the ancestral. Panel C. Comparison of the one-step growth curve between D/CS637/11 CT135-positive and CT135-negative strains. Cells grown in the same conditions were infected at a MOI of 1, and cell scrapings were collected over time after infection for analysis of inclusion-forming units (IFUs). The black line represents the evolved CT135-negative D/CS637/11 strain, whereas the grey line represents the ancestor CT135-positive strain. The shaded area indicates the time points chosen for RNA-seq differential expression comparative analyses.

Fig 5. CT135 mRNA decay analysis.

Panel A shows the comparison of the relative amount of transcripts at 4 h post-infection (pi) and after 10 min of transcriptional blockage with rifampicin (10 μg/ml) between the ancestral (grey) and the evolved populations (black). The assay was performed for all strains with emergent CT135-negative clones (i.e., all epithelial-genital isolates) and for the strain L2b/CS19/08 (control). The number of transcripts was quantified by independent RT-qPCR targeting the two genes of the operon CT134-CT135 (see methods for details), except for serovar D strain as the evolved population lacks CT135. Data was normalized against the number of C. trachomatis genomes quantified on the corresponding DNA samples. In order to facilitate the comparative analysis, the normalized value before rifampicin treatment (4h pi) was arbitrarily set to 1. Panel B shows the mRNA half-life times calculated based on the fit of an exponential decay between the quantified values at 4h pi and the values calculated 10 minutes after the transcriptional arrest.

Fig 6. Comparative analysis of global gene expression (RNA-seq) between D/CT135-positive and D/CT135-negative populations.

Panels A-B. Comparison of gene expression between biological replicates for the D/CT135-positive (A) and D/CT135-negative (B) populations. Pearson correlation coefficients are shown. Panel C. Comparison of gene expression between the D/CT135-negative and D/CT135-positive populations. The red points mark genes and the non-coding RNA for which the fold change of expression exceeds two-fold and the FDR-corrected P-values were below 0.05. For panels A to C, axes are log₁₀-transformed normalized expression levels (FPKM). Panel D. Volcano plot of –log₂ fold change (D/CT135-positive versus D/CT135-negative) versus –log₁₀ adjusted P-values. In order to better fit the scale to data, corrected P-values ≤10⁻³ were set as 10⁻³. Points in red indicate genes and the non-coding RNA for which the fold change of expression exceeds two-fold and the FDR-corrected P-values were below 0.05.