Transcriptional profiling of Chlamydia trachomatis and its host in an ex vivo endocervical primary cell culture system using dual RNA sequencing

Abstract

Chlamydia trachomatis (Ct) is an obligate intracellular bacterium that causes significant ocular and urogenital morbidity worldwide. Understanding host-pathogen interactions is challenging but dual RNA sequencing offers simultaneous transcriptome data for comprehensive interrogations into these interactions. While transcriptional profiling of both Ct and host-derived immortalized cells has been performed, this study used dual RNA sequencing to examine host-pathogen interactions in ex vivo human primary endocervical stromal cells infected with Ct strain E/Bour. At 1-hour post-infection (1hpi), 168 differentially expressed host genes (DEGs) were identified, 40% of which were non-coding RNAs, novel proteins, or pseudogenes. Pathway analysis revealed no significant enrichment at this stage, indicating a quiescent host response. At 24hpi, 212 DEGs were identified, with strong upregulation of interferon-stimulated genes and activation of the cGAS-STING and RLR pathways, despite the absence of detectable type I interferons. Pro-inflammatory and leukocyte recruitment genes were also highly expressed, suggesting an immunoreactive phenotype at this later stage. Ct transcriptomics identified 331 early and 903 mid-infection genes. Inclusion-membrane genes peaked at 1hpi, while hemolysin-like and polymorphic membrane protein genes were upregulated at 24hpi. Enrichment analysis identified pathways related to catalytic activity, host modulation, and bacterial survival. This study demonstrates distinct temporal dynamics in Ct-host interactions, including early host immune quiescence and robust mid-infection activation of innate immunity in contradistinction to previous host and Ct findings in immortalized cell lines. The findings emphasize the utility of ex vivo human primary cell culture for investigating Ct pathogenesis using clinically relevant Ct strains and provide a foundation for future exploration of uncharacterized genes and pathways critical to Ct infection.

Introduction

Chlamydia trachomatis (Ct) is an obligate intracellular bacterium and the leading cause of sexually transmitted infections (STI) worldwide (Sexually transmitted infections (STIs)). Ct STIs are frequently asymptomatic, and if left untreated in the female genital tract, can lead to pelvic inflammatory disease, ectopic pregnancy, and infertility (Land et al., 2010; Rowley et al., 2019; Sexually Transmitted Disease Surveillance, 2023). Ct exhibits a biphasic developmental cycle that alternates between an infectious, non-replicative elementary body (EB) and a replicative, metabolically active reticulate body (RB). However, if RBs are exposed to host environmental stressors such as antibiotics, RBs can halt replication and enter a reversible persistent state that prioritizes cell functions required for long-term survival (Panzetta et al., 2018).

Research on the dynamics of Ct gene expression have increasingly substantiated how Ct modulates host-pathogen interactions that drive pathogenesis. Time course experiments have previously identified transcriptional changes throughout Ct development, categorizing them into temporal groups: Early (1–3 hours) to mid-early (>3–8 hours), mid (>8–24 hours), and late (>24 hours) gene expression (Belland et al., 2003b; Nicholson et al., 2003; Wurihan et al., 2024). One in vitro modelling study evaluated gene activity during persistence by monitoring Ct transcriptional responses to interferon (IFN)γ and identified key pathways that were differentially regulated (Belland et al., 2003a). For instance, Ct genes involved in phospholipid utilization and protein translation were found to be upregulated (Belland et al., 2003a). A rigorous understanding of the mechanisms that maintain and regulate the Ct transcriptome within the context of the host response can greatly enhance the development of diagnostic tools and preventive strategies for Ct.

Dual RNA-sequencing (dual RNA-seq) provides the simultaneous interrogation of both the host and pathogen transcriptomes (Westermann et al., 2017). Despite the potential of this technology, its application in chlamydial pathogenesis has thus far been limited to the study of Ct infection in immortalized cells (Humphrys et al., 2013; Hayward et al., 2021). Immortalized cells are unlikely to represent the actual in vivo response to pathogen infection. Building on our previous study of the human primary endocervical cell model of Ct infection, we found that human primary endocervical stromal (ECS) cells, unlike primary endocervical epithelial cells, elicited a variety of responses including pro- and anti-inflammatory mediators that varied depending on the infecting Ct ocular or urogenital strains (Jolly et al., 2019). We therefore applied dual RNA-seq to investigate the host-pathogen dynamics of Ct infection of ECS cells with urogenital reference strain E/Bour and mock infection of ECS at the early (i.e., 1 hour post infection [1hpi]) and mid (24hpi) timepoints of development. In the present study, we found unique host-pathogen responses compared to the literature that can be used to further elucidate Ct pathogenic mechanisms.

Materials and methods

C. trachomatis infection of primary cells

ECS cells were prepared from biopsied cancer-free de-identified hysterectomy tissue as previously described (Jolly et al., 2019). Briefly, tissue plugs from a single donor were plated the same day as the tissue biopsy on fresh culture dishes, grown in RPMI 1640 + 10% FBS Media, and fibroblasts were collected 1 to 2 weeks later, as they migrated from the tissue plug during this period. Confirmation of the cell type was based on fibronectin-specific antibody staining (Sheep anti-human fibronectin antibody, R&D systems; Pittsburgh, PA) of the cells. ECS cells used in this experiment were at passage number 2-3. ECS cells were plated in two 6-well tissue culture dishes (9.6 cm² per well) in RPMI 1640 + 10% FBS Media (Lonza, Allendale, New Jersey) in the absence of antibiotics/antifungal agents. One plate contained two triplicate experiments for each condition (infected and mock-infected) at one timepoint. The timepoints of 1hpi and 24hpi were tested. Cells were grown to 80% confluency as confirmed by brightfield microscopy and then used to calculate the volume of inoculum required to achieve the target multiplicity of infection (MOI) of 1. The cells were mock-infected or infected with Ct reference strain E/Bour to achieve an MOI of 1 and centrifuged for 1 hour at 1500 rpm, 37°C, to facilitate bacterial contact with the cells. The infected cells were further incubated for one hour or 24 hours at 37°C with 5% CO₂. At 24hpi, immediately prior to harvesting, the cells were visualized using a Nikon T1-E real time microscope with 5% CO₂ under high-resolution imaging to confirm an MOI of 1 similar to our previous study (Jolly et al., 2019). Cells were washed twice with DPBS (ThermoFisher, Waltham, MA), harvested in TRI reagent (Zymo Research Corp, Irvine, CA) and immediately subjected to RNA extraction.

Total RNA isolation, rRNA and polyA depletion, and library preparation and sequencing

Total RNA was extracted with Zymo Research Direct-zol kit (Zymo Research Corp, Irvine CA) following manufacturer’s instructions. RNA samples were frozen in -80°C for storage and then transferred to the sequencing facility (Functional Genomics Laboratory, UC Berkeley, Berkeley, CA). RNA concentrations of total RNA samples were measured using Qubit™ RNA HS Assay Kit (Invitrogen; Carlsbad, CA), and RNA integrity was measured using a Bioanalyzer 2100 Expert and the Agilent RNA 6000 Pico assay (Santa Clara, CA). All samples had an RNA integrity number greater than 7.0. To enrich samples for mRNA, all total RNA samples underwent ribosomal RNA depletion using Illumina Human and bacteria Ribo-Zero rRNA Removal Kits (Human/Mouse/Rat, and Bacteria; San Diego, CA). At this point, rRNA depleted RNA from mock-infected samples were sent for the library preparation phase of the procedure. rRNA depleted RNA from infected samples were split into two portions. One portion was sent directly for the library preparation phase of the procedure while the other portion was enriched for prokaryotic mRNA using the KAPA mRNA Capture Kit (Roche, San Francisco, CA) prior to sequencing. Illumina mRNA sequencing libraries were prepared using the KAPA mRNA HyperPrep Kit (Roche Sequencing Solutions, Pleasanton, CA) as per manufacturer’s instructions. Prior to sequencing, the quality of the libraries was verified on an AATI Fragment Analyzer (Agilent). Individually indexed RNA libraries were pooled and sequenced with Illumina NovaSeq 6000 using 150 bp paired-end reads, yielding an approximate 10 Gb per sample. Sequencing data were demultiplexed and analyzed independently. A simplified workflow from RNA extraction to sequencing is shown in Supplementary Figure S1.

RNA sequence read processing, quantification and gene expression analysis

Host transcript quantification was performed using the Nextflow rnaseq package (Patel et al., 2024) with default settings. To quantify the bacterial reads, quality control on sequenced reads was performed using Fastp (Chen, 2023). Reads were subsequently depleted of contaminating host reads using Scrubby (Steinig, 2024). With Scrubby, two host read removal methods were implemented: Kraken2, with taxa set as Chordata, and Minimap2 targeting the human genome (GRCh38). Host-depleted reads were subsequently run through the Nextflow rnaseq package (Patel et al., 2024) for read quantification using the Ct E/bour genome assembly (NC_020971.1 and NC_020947.1) as reference genome. Ct E/bour gene identifiers were compared and converted to Ct D/UW-3/CX gene identifiers where necessary.

Differential gene expression analysis of host and bacterial reads was performed using the DESeq2 (Love et al., 2014) package on the iDEP (Ge et al., 2018) online platform with default settings. This included normalization to account for differences in library size across all samples. Genes with zero counts across all samples were filtered prior to analysis. To select differentially expressed genes for downstream analysis, a cut-off of a false discovery rate (FDR) adjusted p-values of ≤ 0.05 and log₂FC ≥ 2 was used, except otherwise stated. To identify enriched biological pathways, a reactome analysis was performed using gProfiler (Raudvere et al., 2019), STRING (Szklarczyk et al., 2018) and ShinyGO (Ge et al., 2020).

Results

Host transcriptome of primary cells during C. trachomatis infection

Differentially expressed genes (DEGs) were defined as genes meeting the minimum cut-off of log₂ fold change (LFC) ≥2 and a false discovery rate FDR of ≤0.05 when mock-infected cells were compared to Ct reference strain E/Bour infected cells. Using this threshold, 168 genes were found to be differentially expressed in ECS cells at 1hpi: 107 genes were upregulated, while 61 genes were downregulated (Figure 1a; Supplementary Table S1a). At an LFC of 23, Small Nucleolar RNA C/D Box 3D (SNORD3D) was the highest upregulated gene. Other top upregulated genes include Acetylserotonin N-Methyltransferase-Like (ASMTL), Small Integral Membrane Protein 11B (SMIM11B) gene and ENSG00000271741, a novel transcript of ZMYM6-ZMYM6NB readthrough. Top downregulated genes include those encoding a novel protein gene (ENSG00000272617) and a novel transcript (ENSG00000284952) in addition to a long non-coding RNA (ENSG00000267520) (Figure 1b). A substantial number of DEGs identified at 1hpi (n=68, ~40%) were either long intergenic non-protein coding RNAs (lincRNA), novel protein genes and transcripts, or pseudogenes (Supplementary Table S1a).

Figure 1

Gene set enrichment analysis of the host reactome was performed on the DEGs present at 1hpi. For this analysis, pathways with an FDR ≤0.05 were considered significant. This analysis identified no significantly enriched pathway at 1hpi. To assess whether our strict log₂FC ≥2.0 for DEGs requirement limited downstream interpretation, we relaxed the cutoff slightly down to 1.5. Re-analysis showed that there was a significant downregulation of signaling especially by Rho GTPases, Miro GTPases and RHOBTB3 (Supplementary Table S2), suggesting that suppression of intracellular and cell-to-cell signaling was crucial to entry and establishment of Ct infection in ECS cells. Genes involved in these pathways include ARHGAP11B, HNRNPC, PPP2R5D, ARHGEF11, RHOBTB3 and CDC42EP3. A similar re-analysis of upregulated genes at 1hpi failed to yield any enriched pathways. This lack of significantly enriched pathways may be due to the abundance of long intergenic non-protein coding RNA, novel proteins and transcripts, or pseudogenes (Supplementary Table S1a).

At 24hpi, 212 DEGs were identified, of which 164 genes were upregulated while 48 genes were downregulated (Figure 2a; Supplementary Table S1b). Six of the top 10 upregulated DEGs at 24hpi were IFN stimulated genes participating in the type I IFN pathway (Figure 1b). When we compared all DEGs at 1hpi and 24hpi, only eight were present at both timepoints, and seven of them had similar magnitudes of up- or down-regulation at both timepoints. These genes include ENSG00000243679, ANKRD33, SMIM10, HNRNPA1L3, ARHGAP11B, ENSG00000273619, SNORA67 and POC1B-GALNT4 (Supplementary Table S3).

Figure 2

Gene set enrichment analysis of upregulated DEGs at 24hpi showed that the most enriched biological process (BP) term was “response to virus” with 28 genes. This occurred principally through genes participating in the cellular response to type I IFN (Figure 3; Supplementary Table S4a). A separate pathway enrichment analysis performed on all upregulated DEGs using the Reactome database showed IFN alpha(α)/beta(β) signaling as the most enriched pathway with the majority of the DEGs being IFN-stimulated genes (ISGs) (Supplementary Table S5a). Surprisingly, neither IFN-α nor IFN-β were expressed in these cells either at 1hpi or 24hpi (Supplementary Table S5b), suggesting that the upregulation of ISGs occurred in the absence of IFN-α and IFN-β. We therefore investigated the possible activation of the cyclic GMP–AMP synthase (cGAS)–stimulator of interferon genes (STING) and retinoic acid-inducible gene I (RIG-I)-like receptors (RLRs)-MAV pathways for host recognition of cytosolic chlamydial DNA and RNA, respectively. Seven key genes participating in the upstream processes of the cGAS-STING and RLR-MAV pathways were significantly upregulated without the production of detectable IFN-α and IFN-β mRNA (Figure 4). These genes were DDX58, IFIH1, IRF7, and the OAS family members OAS1–3 and OASL. In total, we identified 19 upregulated DEGs and an additional eight upregulated genes with LFC of <2.0 (FDR ≤ 0.05) participating in the cGAS-STING and RLR-MAV pathways (Figure 4).

Figure 3

Figure 4

The BP term “response to molecule (lipopolysaccharide) of bacterial origin” was significantly enriched at 24hpi (Supplementary Table S4a). Of the 411 genes that are known to as participate in this term, 241 were expressed in ECS cells at 24hpi of which 68 met the FDR ≥0.05 criterion. Of these 68 DEGs, 10 were downregulated while 62 were upregulated. By further applying the LFC of ≥2.0, we identified 15 genes that were differentially expressed, and all of them were significantly upregulated. These genes included CCL5, CMPK2, CXCL10, CXCL2, PTGS2, IL1B, FOS, CXCL8, LGALS9, IL24, CCL2, CSF2, PLSCR1, NCF2 and IRAK2 (Supplementary Table S4b). Of note, five of these upregulated genes (CCL2, CCL5, CSF2, CXCL10, CXCL2) are involved in leukocyte recruitment and activation.

C. trachomatis transcriptome during infection of primary endocervical cells

Analyses of the Ct bacterial transcriptome at 1hpi and 24hpi provides an opportunity to understand host-pathogen interactions within the same host cell infection system. Differential gene expression analyses at these two timepoints in primary cells can aid in discerning mechanisms of Ct pathogenesis that are more relevant for human infection.

To adequately capture Ct transcripts, two depletion methods were utilized for the RNA samples obtained from each timepoint: a) rRNA depletion only; used to study both Ct and host transcriptome profiles; and b) rRNA depletion with polyA depletion; used to study Ct transcriptome profiles only. At 1hpi, the rRNA depletion method generated an average of 24,191 Ct reads per sample (total = 145,143), while rRNA and polyA depletion generated an average of 39,083 Ct reads per sample (total = 234,500) (Figure 5a). At 24 hpi Ct reads from the rRNA depletion method averaged 1,429,227 reads per sample (total = 7,146,133), while rRNA and polyA depletion averaged 2,533,238 reads per sample (total = 12,666,188) (Figure 5b). A Euclidean distance between variance stabilized read count for all samples was computed, and the resulting clustering showed a distinct clustering between 1hpi and 24hpi Ct reads (Figure 5c). Taken together, a total of 379,643 and 19,812,321 Ct reads were available at 1hpi and 24hpi respectively for differential expression analysis. At 1hpi, Ct genome coverage ranged from 4.7x to 14.8x (mean: 9.1x; standard deviation: 3.5), while at 24hpi, coverage ranged from 243.8x to 1419.6x (mean: 712.6x; standard deviation: 343.9).

Figure 5

Using a minimum count per million of 0.5 in at least 5 samples, a total of 338 Ct genes were identified as expressed at 1hpi representing about a third of the known Ct transcriptome. To identify genes with high expression levels, the minimum count per million threshold was increased to 100. With this cut-off, 184 genes were found to be highly expressed at 1hpi (Supplementary Table S7a). The top 10 highly expressed early genes were ompA, encoding the most abundant outer membrane protein; npt1, encoding an ATP/ADP exchange transporter involved in host energy parasitism; CT147, encoding an inclusion membrane protein likely involved in avoiding fusion with host lysosomes; CT529, encoding a putative inclusion membrane protein; cpoS (CT229), encoding an inclusion membrane protein important in host trafficking system modification; omcB, encoding the conserved outer membrane complex protein; rpoB and rpoC, encoding DNA-directed RNA polymerase subunits; fusA, encoding elongation factor G; and ahpC, encoding an antioxidant thiol-specific peroxidase. A complete list of the highly expressed early genes is shown in Supplementary Table S7a.

Using the same thresholds as 1hpi, 914 genes were expressed at 24hpi, of which 632 were found to be highly expressed. Not surprisingly, four early genes were also among the top 10 highly expressed genes at 24hpi as they are primarily constitutively expressed (ompA, rpoB, rpoC and fusA). Other members of the top 10 highly expressed genes were nrdA, a gene encoding a ribonucleoside-diphosphate reductase subunit alpha that catalyzes the biosynthesis of deoxyribonucleotides; glgA, a glycogen synthase gene; CT691, a TIGR00153 family protein gene; rpsA, a gene encoding the 30S ribosomal protein S1; ndk: a nucleoside-diphosphate kinase gene; and CT708 that encodes the DEAD/DEAH box helicase (Supplementary Table S7b).

When both timepoints were compared with the cut-off criteria of FDR ≤0.05, Log₂FC ≥2.0, and CPM ≥0.5, 462 DEGs were identified (Figure 5d). These included 71 genes that were highly expressed at 1hpi, and 393 genes at 24hpi. Top DEGs at 1hpi include CT867, a gene encoding an effector protein with deubiquitinase enzyme activity; cpoS (CT229); and the inclusion (inc) gene CT249 (Figure 5e). Of note, while only 11 tRNA genes were differentially expressed at 1hpi relative to 24hpi, they made up 50% (10/20) of the top 20 DEGs. These tRNA are tRNA-Ser (3 copies), tRNA-Glu, tRNA-Leu (2 copies), tRNA-Val, tRNA-Arg, tRNA-Gln, tRNA-Asp, and tRNA-Phe (Supplementary Table S6). At 24hpi, the top DEGs include recJ, a DNA specific exonuclease gene; pmpC, a polymorphic outer membrane protein (Pmp) gene; and CT708 (Figure 5e). We also identified 60 DEGs that were hypothetical proteins (Supplementary Table S6).

Gene set enrichment analysis of Ct DEGs with higher expression at 1hpi showed “Secreted and Transmembrane” as the most enriched term (Figure 6a; Supplementary Table S8a). DEGs in this term included seven inclusion membrane protein genes (incB, incC, incD, incE, cpoS, CT227, and CT228). A manual search of the literature showed a total of 30 inc genes were differentially expressed in this study, 21 (68%) of which had a significantly higher expression level at 1hpi relative to 24hpi (Table 1). Non-inclusion proteins (n=9) known to be secreted into the host cytoplasm were found to be differentially expressed in this study. Of note, only three of these had a significantly higher expression at 1hpi (CT473: a hypothetical protein gene; CT867 and CT868: deubiquitinase genes) (Table 1).

Figure 6

Gene set enrichment analysis of Ct DEGs with higher expression at 24hpi revealed “catalytic activity” as the most enriched term (Figures 6b, d). Other enriched terms include “ion binding,” “carbohydrate metabolic process,” and “modulation of process of other organisms” (Supplementary Table S8b). As “modulation of process of other organisms” was the highest enriched term with 100% of its associated genes differentially expressed, we investigated this process further. Of the 11 DEGs identified in this process, eight were polymorphic membrane protein (pmp) genes: pmpB to pmpI. The other three were hemolysin family protein genes (CT256, CT257 and CT423).

Discussion

An understanding of the complex host-pathogen interactions between the human host and Ct, especially in the early stages of infection, has the potential to aid in identifying biomarkers and develop diagnostic tools and therapeutic targets for preventive strategies. Previous attempts to study transcriptome-wide host-Ct interactions utilized either microarray or RNA-seq technologies and lab-adapted immortalized cells such as HeLa229 and HEp2 (Belland et al., 2003b; Humphrys et al., 2013). While immortalized cells are excellent models and provide an opportunity for hypothesis generation/testing, they are unlikely to represent the actual in vivo response to the pathogen. To date, no transcriptome-wide analysis of primary human cell responses to Ct infection and the simultaneous responses of Ct to those cells have been reported. This study is therefore the first to describe the transcriptional profiling of Ct—in particular strain E that represents the most prevalent strain responsible for STIs in women—and host cells in an ex vivo human primary cell culture system. We show that clinically relevant data that are essential for interrogating Ct pathogenesis are obtainable using this model.

Previous work describing the early host transcriptome following Ct infection reported differential regulation of 13 host genes in HeLa cells infected with Ct strain L₂/434/Bu at 2hpi using microarray technology (Xia et al., 2003), and 622 host genes in Hep2 cells infected with E/Bour at 1hpi based on dual RNA-seq (Humphrys et al., 2013). In the current study, 168 host genes were differentially expressed in E/Bour infected ECS cells at 1hpi, 107 of which were upregulated while 67 were downregulated, representing 3.7x fewer DEGs compared to the prior study (Humphrys et al., 2013). We suspect that the lower number of DEGs in our study points to a less aggressive primary cell response to Ct infection compared to immortalized cells. This would be consistent with the stronger upregulation of extracellular matrix remodeling in immortalized cells (Humphrys et al., 2013) compared to primary cells.

The top host upregulated gene at 1hpi was SNORD3D that encodes a member of the box C/D small nucleolar ribonucleoproteins (snoRNPs). These proteins catalyze site-directed modifications, typically the 2′-O-methylation of ribosomal RNA (Baldini et al., 2021). While their role in host immune responses are not well known, a recent study showed that increased box C/D snoRNPs corresponded to a repressed immune function in the Caenorhabditis elegans model of Pseudomonas aeruginosa infection (Tjahjono et al., 2022). In Ct, SNORD3D upregulation suggests that the organism may influence host RNA processing and ribosome biogenesis to suppress host innate immune responses and thereby evade detection. Box C/D snoRNPs maintain mitochondrial dynamics that include function and surveillance (Tjahjono et al., 2022). Ct manipulates mitochondrial fusion to regulate host cell metabolism and acquire ATP and control mitochondrial fission to prevent apoptosis (Yang et al., 2022), both of which are essential for Ct growth/reproduction and survival, respectively. However, the lack of specific studies on the role of box C/D snoRNPs in chlamydial pathogenesis highlights the need for further research.

In addition to multiple novel protein genes and pseudogenes, we found several host DEGs at 1hpi that were lincRNAs. The biological role of most lincRNAs are largely unknown, however, some have been shown to participate in either negative or positive modulation of host responses to bacterial infection (Agliano et al., 2019; Arunima et al., 2023). Recent data show that lincRNAs may promote Ct replication in vitro (reviewed in (Arunima et al., 2023)). However, none of the nine lincRNAs identified in the current study have been studied in the context of Ct pathogenesis. The absence of functional characterization for many of these lincRNAs emphasizes the need for future studies to determine their precise contributions to Ct infection dynamics, particularly their potential as biomarkers or therapeutic targets.

At 24hpi, the host transcriptome primarily represented an immunoreactive phenotype accompanied by upregulated intracellular signaling. Alongside pro-inflammatory cytokines and adhesion molecules, the immune response was characterized by strong upregulation of genes linked to IFNα and IFNβ and, to a lesser extent, the IFNγ signaling responses. Since none of the IFNs were differentially expressed in infected cells, it is possible that Ct-related CDNs and/or RNA/DNA may be sufficient to trigger an IFN signaling cascade with a bias towards type I IFNs. The role of type I IFN signaling in the clearance of bacteria is enigmatic as both IFNα and IFNβ seem to have opposite effects (McNab et al., 2015; Boxx and Cheng, 2016). A recent study showed that type I IFNs did not prevent infection or replication of Ct in primary synovial fibroblasts (Di Pietro et al., 2020) similar to our findings. A separate study showed that type I IFNs may work in concert with the type II IFN response to inhibit Ct growth in epithelial cells (Ishihara et al., 2005). However, depending on the levels of type I IFNs, B cells may be inhibited along with reducing macrophage responses to IFNγ thereby facilitating infection (McNab et al., 2015). While it is evident that Ct does not completely evade immune responses, it remains unclear if the induction of a strong type I IFN response by ECS cells at 24hpi is sufficient to clear infection, facilitate infection, or encourage persistence by regulating indoleamine 2,3-dioxygenase (IDO), the induction of which depletes tryptophan, an essential amino acid required for Ct growth and replication (Puccetti, 2007; Trinchieri, 2010; Bommana et al., 2021).

The cGAS-STING and RLRs-MAV pathways are important for the recognition of cytosolic chlamydial DNA and RNA, respectively, by host cells (Wen and Li, 2020) and typically end in the production of type I IFN genes and autophagy. The Ct ybbP gene has been shown to be crucial in the synthesis of cyclic di-AMP, a major ligand in the cGAS-STING pathway (Barker et al., 2013; Meier et al., 2023), and was upregulated 16-fold at 24hpi (Supplementary Table S6). However, neither IFNα nor IFNβ were expressed at 24hpi by the primary cells. Surprisingly, a strong stimulation of multiple type I IFN-stimulated genes (ISGs) was observed. There is evidence for cooperation between the cGAS-STING and RLRs-MAV pathways for ISG production (Zevini et al., 2017). Data from the current study suggests that at 24hpi, the main response of primary cells to Ct was the induction of ISGs through the cGAS-STING and RLR pathways (Figure 4). The lack of detectable type I IFN gene expression is intriguing and might suggest a deliberate evasion strategy. Ct may inhibit upstream IFN production to prevent a full-scale immune response, while allowing ISG activation to maintain host cell viability and metabolic activity. This selective modulation could also help Ct evade inflammatory cell death pathways that are typically triggered by robust IFN responses.

A recent study using a murine L929 immortalized fibroblast cell line infected with Ct strain L₂ reported the expression of over 700 Ct genes as early as 1hpi, suggesting widespread early transcriptional activity (Wurihan et al., 2024). In contrast, our study identified 338 Ct transcripts at 1hpi in ECS cells infected with the E/Bour strain. Several key methodological differences likely contribute to this discrepancy. Notably, the previous study employed a high MOI of 50, whereas we employed an MOI of 1 to better reflect physiologically relevant infection burden. While a high MOI may not accelerate developmental timing (Chiarelli et al., 2020), it may influence the number of infected cells resulting in enhanced transcript detection. Additionally, differences in host cell type (i.e., immortalized versus primary cells) as well as Ct strain further influence early gene expression profiles as is evident in our study. Together, these factors highlight the importance of infection context in interpreting pathogen transcriptional dynamics.

Inclusion membrane proteins are a known group of Ct proteins that mediate Ct-host interactions (reviewed in (Bugalhão and Mota, 2019)). For instance, incD redirects host lipids away from their intended host destinations and towards the chlamydial inclusion (Derré et al., 2011; Agaisse and Derré, 2014), while cpoS (CT229) controls host cell vesicular trafficking (Bugalhão and Mota, 2019). In the current study, we report the upregulation of 21 known inc genes at 1hpi compared to 24hpi in primary cells (Table 1). In contrast, only 9 inc genes had higher expression levels at 24hpi compared to 1hpi. While the functions of the majority of inc genes identified in this study are still unknown, we suspect that they are crucial to a successful infection event and subsequent establishment of the Ct inclusion in primary host cells.

This study also underscores the complementary roles of upregulated hemolysin-like proteins (encoded by CT256, CT257, CT423) and PmpB to PmpI in the pathogenesis of Ct during the mid-infection stage (24hpi). Given their known roles in other bacteria, hemolysin-like proteins likely facilitate Ct survival by modifying the host intracellular environment, potentially through membrane disruption, nutrient acquisition, and suppression of immune responses (Aroian and van der Goot, 2007; Gonzalez et al., 2011; Wiles and Mulvey, 2013; Mondal and Chattopadhyay, 2022). In parallel, PmpB to PmpI are critical for mediating host-pathogen interactions such as adhesion, intracellular trafficking, and immune evasion (Vasilevsky et al., 2016; Pickering et al., 2017; Jury et al., 2023). Their concurrent expression suggests a coordinated strategy where hemolysin-like proteins create conditions favorable for Pmp-mediated processes. For instance, hemolysins may alter autophagy or inhibit lysosomal degradation, enabling Pmps to enhance Ct adhesion, endosomal trafficking and inclusion stability as well as other interactions with host cellular pathways. Both groups appear to contribute to Ct’s ability to evade host immune defenses while ensuring access to resources needed for replication. Functional studies investigating the co-regulation and physical interactions of these proteins could provide insights into their mechanistic roles and reveal potential targets for therapeutic intervention.

We also identified 60 DEGs that were hypothetical proteins with unknown functions. This observation suggests that the transcriptome of Ct at the early infection stage is currently under-characterized thus hampering an adequate understanding of their roles in Ct pathogenesis. Describing the role of these genes and their proteins may provide more insights into the natural course of Ct pathogenesis.

This study offers novel insights into host and Ct transcriptional dynamics in ECS cells, although some limitations should be noted. Similar to human immortalized cells that are derived from a single donor, our ECS cells were obtained from a single donor. Further transcriptomic analyses in ECS cells from additional donors are warranted to expand on our findings and capture any inter-individual variations. Further, while we previously performed extensive head-to-head comparisons between Ct infected primary and immortalized cells that informed our experimental approach and framework for interpreting our findings (Jolly et al., 2019), we did not do so here. These types of experiments could be performed to directly compare specific cell-type responses. Finally, protein-level quantitative assays on cell culture supernatants would be valuable to further refine and substantiate the transcriptomic data.

In conclusion, our study highlights the utility of ex vivo human primary cell cultures and dual RNA-seq in uncovering unique host-pathogen dynamics during Ct infection. Early infection showed a subdued host response alongside elevated expression of Ct inclusion genes, while mid-infection revealed robust host immune activation, including upregulated interferon-stimulated genes and cGAS-STING and RLR pathway activation. On the bacterial side, mid-stage expression emphasized hemolysin-like proteins and Pmps. These findings underscore the value of physiologically relevant ex vivo models for studying Ct pathogenesis that include the use of clinically relevant Ct strains and identifying novel therapeutic targets. Future research should explore the roles of uncharacterized genes and validate these results in vivo to deepen our understanding of Ct biology, pathogenesis, and host responses.

References

• 1

  AgaisseH.DerréI. (2014). Expression of the effector protein IncD in Chlamydia trachomatis mediates recruitment of the lipid transfer protein CERT and the endoplasmic reticulum-resident protein VAPB to the inclusion membrane. Infect. Immun.82, 2037–2047. doi: 10.1128/IAI.01530-14

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AglianoF.FitzgeraldK. A.VellaA. T.RathinamV. A.MedvedevA. E. (2019). Long non-coding RNA lincRNA-EPS inhibits host defense against listeria monocytogenes infection. Front. Cell Infect. Microbiol.9, 481. doi: 10.3389/fcimb.2019.00481

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  AroianR.van der GootF. G. (2007). Pore-forming toxins and cellular non-immune defenses (CNIDs). Curr. Opin. Microbiol.10, 57–61. doi: 10.1016/j.mib.2006.12.008

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  ArunimaA.van SchaikE. J.SamuelJ. E. (2023). The emerging roles of long non-coding RNA in host immune response and intracellular bacterial infections. Front. Cell Infect. Microbiol.13, 1160198. doi: 10.3389/fcimb.2023.1160198

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  (2023). Sexually Transmitted Disease Surveillance. Available online at: https://www.cdc.gov/sti-statistics/?CDC_AAref_Val=https://www.cdc.gov/std/statistics/2021/default.htm (Accessed July 12, 2023).

  • Google Scholar
  • View reference in article

• 6

  Sexually transmitted infections (STIs). Available online at: https://www.who.int/news-room/fact-sheets/detail/sexually-transmitted-infections-(stis) (Accessed January 15, 2025).

  • Google Scholar
  • View reference in article

• 7

  BaldiniL.CharpentierB.LabialleS. (2021). Emerging data on the diversity of molecular mechanisms involving C/D snoRNAs. Non-Coding RNA7 (2). doi: 10.3390/ncrna7020030

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  BarkerJ. R.KoestlerB. J.CarpenterV. K.BurdetteD. L.WatersC. M.VanceR. E.et al. (2013). STING-dependent recognition of cyclic di-AMP mediates type I interferon responses during Chlamydia trachomatis infection. mBio4, e00018–e00013. doi: 10.1128/mBio.00018-13

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  BellandR. J.NelsonD. E.VirokD.CraneD. D.HoganD.SturdevantD.et al. (2003a). Transcriptome analysis of chlamydial growth during IFN-gamma-mediated persistence and reactivation. Proc. Natl. Acad. Sci. U S A.100, 15971–15976. doi: 10.1073/pnas.2535394100

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  BellandR. J.ZhongG.CraneD. D.HoganD.SturdevantD.SharmaJ.et al. (2003b). Genomic transcriptional profiling of the developmental cycle of Chlamydia trachomatis. Proc. Natl. Acad. Sci. U S A.100, 8478–8483. doi: 10.1073/pnas.1331135100

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  BommanaS.SomboonnaN.RichardsG.TarazkarM.DeanD. (2021). Tryptophan Operon Diversity Reveals Evolutionary Trends among Geographically Disparate Chlamydia trachomatis Ocular and Urogenital Strains Affecting Tryptophan Repressor and Synthase Function. mBio12 (3), 10.1128/mbio.00605-21. doi: 10.1128/mbio.00605-21

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  BoxxG. M.ChengG. (2016). The roles of type I interferon in bacterial infection. Cell Host Microbe19, 760–769. doi: 10.1016/j.chom.2016.05.016

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  BugalhãoJ. N.MotaL. J. (2019). The multiple functions of the numerous Chlamydia trachomatis secreted proteins: the tip of the iceberg. Microb. Cell Graz Austria.6, 414–449. doi: 10.15698/mic2019.09.691

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  ChenS. (2023). Ultrafast one-pass FASTQ data preprocessing, quality control, and deduplication using fastp. iMeta2, e107. doi: 10.1002/imt2.v2.2

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  ChiarelliT. J.GrieshaberN. A.OmslandA.RemienC. H.GrieshaberS. S. (2020). Single-inclusion kinetics of chlamydia trachomatis development. mSystems5 (5). doi: 10.1128/mSystems.00689-20

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  DerréI.SwissR.AgaisseH. (2011). The lipid transfer protein CERT interacts with the Chlamydia inclusion protein IncD and participates to ER-Chlamydia inclusion membrane contact sites. PloS Pathog.7, e1002092. doi: 10.1371/journal.ppat.1002092

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  Di PietroM.FilardoS.FrascaF.ScagnolariC.ManeraM.SessaV.et al. (2020). Interferon-γ Possesses anti-microbial and immunomodulatory activity on a chlamydia trachomatis infection model of primary human synovial fibroblasts. Microorganisms8, 235. doi: 10.3390/microorganisms8020235

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  GeS. X.JungD.YaoR. (2020). ShinyGO: a graphical gene-set enrichment tool for animals and plants. Bioinforma Oxf Engl.36, 2628–2629. doi: 10.1093/bioinformatics/btz931

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  GeS. X.SonE. W.YaoR. (2018). iDEP: an integrated web application for differential expression and pathway analysis of RNA-Seq data. BMC Bioinf.19, 534. doi: 10.1186/s12859-018-2486-6

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  GonzalezM. R.BischofbergerM.FrêcheB.HoS.PartonR. G.van der GootF. G. (2011). Pore-forming toxins induce multiple cellular responses promoting survival. Cell Microbiol.13, 1026–1043. doi: 10.1111/j.1462-5822.2011.01600.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  HaywardR. J.HumphrysM. S.HustonW. M.MyersG. S. A. (2021). Dual RNA-seq analysis of in vitro infection multiplicity and RNA depletion methods in Chlamydia-infected epithelial cells. Sci. Rep.11, 10399. doi: 10.1038/s41598-021-89921-x

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  HumphrysM. S.CreasyT.SunY.ShettyA. C.ChibucosM. C.DrabekE. F.et al. (2013). Simultaneous transcriptional profiling of bacteria and their host cells. PloS One8, e80597. doi: 10.1371/journal.pone.0080597

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  IshiharaT.AgaM.HinoK.UshioC.TaniguchiM.IwakiK.et al. (2005). Inhibition of chlamydia trachomatis growth by human interferon-alpha: mechanisms and synergistic effect with interferon-gamma and tumor necrosis factor-alpha. BioMed. Res. Tokyo Jpn.26, 179–185. doi: 10.2220/biomedres.26.179

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  JollyA. L.RauS.ChadhaA. K.AbdulraheemE. A.DeanD. (2019). Stromal fibroblasts drive host inflammatory responses that are dependent on chlamydia trachomatis strain type and likely influence disease outcomes. mBio10, e00225–e00219. doi: 10.1128/mBio.00225-19

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  JuryB.FlemingC.HustonW. M.LuuL. D. W. (2023). Molecular pathogenesis of Chlamydia trachomatis. Front. Cell Infect. Microbiol.13. doi: 10.3389/fcimb.2023.1281823

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  LandJ. A.Van Bergen JE a.M.MorréS. A.PostmaM. J. (2010). Epidemiology of Chlamydia trachomatis infection in women and the cost-effectiveness of screening. Hum. Reprod. Update.16, 189–204. doi: 10.1093/humupd/dmp035

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  LoveM. I.HuberW.AndersS. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol.15, 550. doi: 10.1186/s13059-014-0550-8

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  McNabF.Mayer-BarberK.SherA.WackA.O’GarraA. (2015). Type I interferons in infectious disease. Nat. Rev. Immunol.15, 87–103. doi: 10.1038/nri3787

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  MeierK.JachmannL. H.TürközG.SaitM. R. B.PérezL.KeppO.et al. (2023). The Chlamydia effector CpoS modulates the inclusion microenvironment and restricts the interferon response by acting on Rab35. mBio14, e03190–e03122. doi: 10.1128/mbio.03190-22

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  MondalA. K.ChattopadhyayK. (2022). Chapter Seven—Structures and functions of the membrane-damaging pore-forming proteins In R. Donev (Ed.). Adv. Protein Chem. Structural Biol. 128, 241–288. doi: 10.1016/bs.apcsb.2021.07.001

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  NicholsonT. L.OlingerL.ChongK.SchoolnikG.StephensR. S. (2003). Global stage-specific gene regulation during the developmental cycle of Chlamydia trachomatis. J. Bacteriol.185, 3179–3189. doi: 10.1128/JB.185.10.3179-3189.2003

  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  PanzettaM. E.ValdiviaR. H.SakaH. A. (2018). Chlamydia Persistence: A Survival Strategy to Evade Antimicrobial Effects in-vitro and in-vivo. Front. Microbiol.9, 3101. doi: 10.3389/fmicb.2018.03101

  • CrossRef
  • Google Scholar

• 33

  PatelH.EwelsP.ManningJ.GarciaM. U.PeltzerA.HammarénR.et al. (2024). nf-core/rnaseq: nf-core/rnaseq v3.18.0 - Lithium Lynx. Zenodo. doi: 10.5281/zenodo.14537300

  • CrossRef
  • Google Scholar

• 34

  PickeringH.TengA.FaalN.JoofH.MakaloP.CassamaE.et al. (2017). Genome-wide profiling of humoral immunity and pathogen genes under selection identifies immune evasion tactics of Chlamydia trachomatis during ocular infection. Sci. Rep.7, 9634. doi: 10.1038/s41598-017-09193-2

  • CrossRef
  • Google Scholar

• 35

  PuccettiP. (2007). On watching the watchers: IDO and type I/II IFN. Eur. J. Immunol.37, 876–879. doi: 10.1002/eji.200737184

  • CrossRef
  • Google Scholar

• 36

  RaudvereU.KolbergL.KuzminI.ArakT.AdlerP.PetersonH.et al. (2019). g:Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res.47, W191–W198. doi: 10.1093/nar/gkz369

  • CrossRef
  • Google Scholar

• 37

  RowleyJ.Vander HoornS.KorenrompE.LowN.UnemoM.Abu-RaddadL. J.et al. (2019). Chlamydia, gonorrhoea, trichomoniasis and syphilis: global prevalence and incidence estimates, 2016. Bull. World Health Organ97, 548–562P. doi: 10.2471/BLT.18.228486

  • CrossRef
  • Google Scholar

• 38

  SteinigE. (2024). esteinig/scrubby. Available online at: https://github.com/esteinig/scrubby (Accessed June 5, 2024).

  • Google Scholar

• 39

  SzklarczykD.GableA. L.LyonD.JungeA.WyderS.Huerta-CepasJ.et al. (2018). STRING v11: protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res.47(D1), D607–D613. doi: 10.1093/nar/gky1131

  • CrossRef
  • Google Scholar

• 40

  TangD.ChenM.HuangX.ZhangG.ZengL.ZhangG.et al. (2023). SRplot: A free online platform for data visualization and graphing. PloS One18, e0294236. doi: 10.1371/journal.pone.0294236

  • CrossRef
  • Google Scholar

• 41

  TjahjonoE.RevtovichA. V.KirienkoN. V. (2022). Box C/D small nucleolar ribonucleoproteins regulate mitochondrial surveillance and innate immunity. PloS Genet.18, e1010103. doi: 10.1371/journal.pgen.1010103

  • CrossRef
  • Google Scholar

• 42

  TrinchieriG. (2010). Type I interferon: friend or foe? J. Exp. Med.207, 2053–2063. doi: 10.1084/jem.20101664

  • CrossRef
  • Google Scholar

• 43

  VasilevskyS.StojanovM.GreubG.BaudD. (2016). Chlamydial polymorphic membrane proteins: regulation, function and potential vaccine candidates. Virulence7, 11–22. doi: 10.1080/21505594.2015.1111509

  • CrossRef
  • Google Scholar

• 44

  WenY.LiZ. (2020). The STING pathway in response to chlamydial infection. Microb. Pathog.140, 103950. doi: 10.1016/j.micpath.2019.103950

  • CrossRef
  • Google Scholar

• 45

  WestermannA. J.BarquistL.VogelJ. (2017). Resolving host-pathogen interactions by dual RNA-seq. PloS Pathog.13, e1006033. doi: 10.1371/journal.ppat.1006033

  • CrossRef
  • Google Scholar

• 46

  WilesT. J.MulveyM. A. (2013). The RTX pore-forming toxin α-hemolysin of uropathogenic Escherichia coli: progress and perspectives. Future Microbiol.8, 73–84. doi: 10.2217/fmb.12.131

  • CrossRef
  • Google Scholar

• 47

  WurihanW.WangY.YeungS.ZouY.LaiZ.FondellJ. D.et al. (2024). Expression activation of over 70% of Chlamydia trachomatis genes during the first hour of infection. Infect. Immun.92 (3), e0053923. doi: 10.1128/iai.00539-23

  • CrossRef
  • Google Scholar

• 48

  XiaM.BumgarnerR. E.LampeM. F.StammW. E. (2003). Chlamydia trachomatis infection alters host cell transcription in diverse cellular pathways. J. Infect. Dis.187(3), 424–434. doi: 10.1086/367962

  • CrossRef
  • Google Scholar

• 49

  YangY.LeiW.ZhaoL.WenY.LiZ. (2022). Insights into mitochondrial dynamics in chlamydial infection. Front. Cell Infect. Microbiol.12, 835181. doi: 10.3389/fcimb.2022.835181

  • CrossRef
  • Google Scholar

• 50

  ZeviniA.OlagnierD.HiscottJ. (2017). Crosstalk between cytoplasmic RIG-I and STING sensing pathways. Trends Immunol.38, 194–205. doi: 10.1016/j.it.2016.12.004

  • CrossRef
  • Google Scholar