Skip to main content
Download PDF

Population-based analysis of ocular Chlamydia trachomatis in trachoma-endemic West African communities identifies genomic markers of disease severity

Genome Medicine volume 10, Article number: 15 (2018) Cite this article

Abstract

Background

Chlamydia trachomatis (Ct) is the most common infectious cause of blindness and bacterial sexually transmitted infection worldwide. Ct strain-specific differences in clinical trachoma suggest that genetic polymorphisms in Ct may contribute to the observed variability in severity of clinical disease.

Methods

Using Ct whole genome sequences obtained directly from conjunctival swabs, we studied Ct genomic diversity and associations between Ct genetic polymorphisms with ocular localization and disease severity in a treatment-naïve trachoma-endemic population in Guinea-Bissau, West Africa.

Results

All Ct sequences fall within the T2 ocular clade phylogenetically. This is consistent with the presence of the characteristic deletion in trpA resulting in a truncated non-functional protein and the ocular tyrosine repeat regions present in tarP associated with ocular tissue localization. We have identified 21 Ct non-synonymous single nucleotide polymorphisms (SNPs) associated with ocular localization, including SNPs within pmpD (odds ratio, OR = 4.07, p* = 0.001) and tarP (OR = 0.34, p* = 0.009). Eight synonymous SNPs associated with disease severity were found in yjfH (rlmB) (OR = 0.13, p* = 0.037), CTA0273 (OR = 0.12, p* = 0.027), trmD (OR = 0.12, p* = 0.032), CTA0744 (OR = 0.12, p* = 0.041), glgA (OR = 0.10, p* = 0.026), alaS (OR = 0.10, p* = 0.032), pmpE (OR = 0.08, p* = 0.001) and the intergenic region CTA0744–CTA0745 (OR = 0.13, p* = 0.043).

Conclusions

This study demonstrates the extent of genomic diversity within a naturally circulating population of ocular Ct and is the first to describe novel genomic associations with disease severity. These findings direct investigation of host-pathogen interactions that may be important in ocular Ct pathogenesis and disease transmission.

Background

The obligate intracellular bacterium Chlamydia trachomatis (Ct) is the leading infectious cause of blindness (trachoma) and the most common sexually transmitted bacterial infection [1, 2].

Ct strains are differentiated into biovars based on pathobiological characteristics and serovars based on serological reactivity for the major outer membrane protein (MOMP) encoded by ompA [3]. Serovars largely differentiate biological groups associated with trachoma (A–C), sexually transmitted disease (D–K) and lymphogranuloma venereum (LGV) (L1–L3). Despite diverse biological phenotypes, Ct strains share near complete genomic synteny and gene content [4], suggesting that minor genetic changes influence pathogen-host and tissue-specific infection characteristics [5,6,7,8]. All published African ocular Ct genomes are situated on the ocular branch within the T2 clade of non-LGV urogenital isolates [4]. Currently there are only 31 published ocular Ct genome sequences [4, 9,10,11,12].

The pathogenesis of chlamydial infection begins with epithelial inflammation and may progress to chronic immunofibrogenic processes leading to blindness and infertility, though many Ct infections do not result in sequelae [13, 14]. Strain-specific differences related to clinical presentation have been investigated in trachoma [8, 15, 16]. These studies examined a small number of ocular Ct isolates from the major trachoma serotypes and found a small subset of genes in addition to ompA that were associated with differences in in vitro growth rate, burst size, plaque morphology, interferon gamma –(IFNγ) sensitivity and, most importantly, intensity of infection and clinical disease severity in non-human primates (NHPs), suggesting that genetic polymorphisms in Ct may contribute to the observed variability in severity of trachoma in endemic communities [8].

The obligate intracellular development of Ct has presented significant technical barriers to basic research into chlamydial biology. Only recently has genetic manipulation of the chlamydial plasmid been possible, allowing in vitro transformation and modification studies, though this remains technically challenging, necessitating alternative approaches [17, 18].

Whole genome sequencing (WGS) has recently been used to identify regions of likely recombination in recent clinical isolates, demonstrating that WGS analysis may be an effective approach for the discovery of loci associated with clinical presentation [6]. Additionally, a number of putative virulence factors have been identified through WGS analysis and subsequent in vitro and animal studies [5, 19,20,21,22,23,24,25,26,27,28,29,30]. However, there are currently no published population-based studies of Ct using WGS with corresponding detailed clinical data, making it difficult to relate genetic changes to functional relevance and virulence factors in vivo.

There is an increasing pool of Ct genomic data, largely from archived samples following cell culture and more recently directly from clinical samples [31]. WGS data obtained directly from clinical samples can be preferable to using WGS data obtained from cell-cultured Ct, since repeated passage of Ct results in mutations that are not observed in vivo [32,33,34].

Ct bacterial load is associated with disease severity, particularly conjunctival inflammation, in active (infective) trachoma [35]. Conjunctival inflammation has previously been shown to be a marker of severe disease and plays an important role in the pathogenesis of scarring trachoma [36,37,38]. In this study we used principal component analysis (PCA) to reduce the dimensions of clinical grade of inflammation (defined using the P score from the follicles, papillary hypertrophy, conjunctival scarring (FPC) trachoma grading system [39]) and Ct bacterial load to a single metric to define an in vivo conjunctival phenotype in active (infective) trachoma. PCA is a recognized dimension reduction technique used to combine multiple correlated traits into their uncorrelated principal components (PCs) [40,41,42], allowing us to examine the relationship between Ct genotype and disease severity. These data from the trachoma-endemic region of the Bijagós Archipelago of Guinea-Bissau currently represent the largest collection of ocular Ct sequences from a single population and provide a unique opportunity to gain insight into ocular Ct pathogenesis in humans.

Methods

Survey, clinical examination and sample collection

Survey, clinical examination and sample collection methods have been described previously [43, 44]. Briefly, we conducted a cross-sectional population-based survey in trachoma-endemic communities on the Bijagós Archipelago of Guinea-Bissau. The upper tarsal conjunctivae of each consenting participant were examined, digital photographs were taken, a clinical trachoma grade was assigned and two sequential conjunctival swabs were obtained from the left upper tarsal conjunctiva of each individual using a standardized method [43]. DNA was extracted and Ct omcB (genomic) copies/swab quantified from the second conjunctival swab using droplet digital polymerase chain reaction (ddPCR) [44, 45].

We used the modified FPC grading system for trachoma [39]. The modified FPC system allows detailed scoring of the conjunctiva for the presence of follicles (F score), papillary hypertrophy (conjunctival inflammation) (P score) and conjunctival scarring (C score), assigning a grade of 0–3 for each parameter. A single validated grader conducted the examinations, and these were verified by an expert grader (masked to the field grades and ddPCR results) using the digital photographs. Grader concordance was measured using Cohen’s kappa, where a kappa > 0.9 was used as the threshold to indicate good agreement.

Conjunctival inflammation (P score) is known to have a strong association with Ct bacterial load in this and other populations [35, 46,47,48,49]. For this study we used PCA to combine the presence of inflammation (defined by the P score using the FPC trachoma grading system [39]) with Ct bacterial load (defined by tertile cut-offs illustrated in Additional file 1: Figure S1) [50]. The conjunctival disease phenotype is a dimension reduction of these two variables, defining what we observed in the conjunctiva at the time of sampling (Fig. 1). Dimension reduction using PCA to define complex disease phenotypes in genome-wide association studies (GWASs) is well recognized, as it allows multiple traits to be included to capture a more complex phenotype and accounts for correlation between traits. This approach therefore may reveal novel loci or pathways that would not be evident in a single-trait GWAS, where the full extent of genetic variation cannot be captured [40].

Fig. 1
figure 1

Composite in vivo conjunctival disease severity phenotype in ocular Chlamydia trachomatis infection. A composite in vivo phenotype was derived using principal component analysis (PCA) for dimension reduction of two phenotypic traits: a disease severity score (using the P score value) and C. trachomatis load (where C. trachomatis load was log transformed and cut-offs determined from the resulting density plot (see Additional file 1: Figure S1)). Each circle represents an individual infection (represented on the x-axis (Index), n = 81). Circle size reflects C. trachomatis load and circle colour reflects inflammatory P score (P0–P3) defined using the modified FPC (follicles, papillary hypertrophy, conjunctival scarring) grading system for trachoma [39]

Full size image

Preparation of chlamydial DNA from cell culture

For eight specimens, WGS data were obtained following Ct isolation in cell culture (from the first conjunctival swab) as a preliminary exploration of Ct genomic diversity in this population. Briefly, samples were isolated in McCoy cell cultures by removing 100 μl eluate from the original swab with direct inoculation onto a glass coverslip within a bijou containing Dulbecco’s modified Eagle’s medium (DMEM). The inocula were centrifuged onto cell cultures at 1800 rpm for 30 min. Following centrifugation the cell culture supernatant was removed and cycloheximide-containing DMEM was added to infected cells which were then incubated at 37 °C in 5% CO2 for 3 days. Viable Ct elementary bodies (EBs) were observed by phase contrast microscopy. Cells were harvested and further passaged every 3 days until all isolates reached a multiplicity of infection between 50 and 90% in 2xT25 flasks. Each isolate was prepared and the EBs purified as described previously [51]. DNA was extracted from the purified EBs using the Promega Wizard Genomic Purification kit according to the manufacturer’s protocol [52].

Pre-sequencing target enrichment

For the remaining specimens (n = 118), WGS data were obtained directly from clinical samples. DNA baits spanning the length of the Ct genome were compiled by SureDesign and synthesized by SureSelectXT (Agilent Technologies, UK). The total DNA extracted from clinical samples was quantified and carrier human genomic DNA added to obtain a total of 3 μg input for library preparation. DNA was sheared using a Covaris E210 acoustic focusing unit [31]. End-repair, non-templated addition of 3′-A adapter ligation, hybridization, enrichment PCR and all post- reaction clean-up steps were performed according to the SureSelectXT Illumina Paired-End Sequencing Library protocol (v1.4.1 Sept 2012). All recommended quality control measures were performed between steps.

Whole genome sequencing and sequence quality filtering

DNA was sequenced at the Wellcome Trust Sanger Institute using Illumina paired-end technology (Illumina GAII or HiSeq 2000). All 126 sequences passed standard FastQC quality control criteria [53]. Sequences were aligned to the most closely related reference genome, Chlamydia trachomatis A/HAR-13 (GenBank accession umber NC_007429.1 and plasmid GenBank accession number NC_007430.1), using the Burrows-Wheeler Aligner (BWA) [54]. SAMtools/BCFtools (SAMtools v1.3.1) [55] and the Genome Analysis Tool Kit (GATK) [56] were used to call SNPs. We used standard GATK SNP calling algorithms, where > 10× depth of coverage is routinely used as the threshold value [56, 57]. This has been shown to be adequate for SNP calling in this context [57,58,59].

Variants were selected as the intersection data set between those obtained using both SNP callers and SNPs were further quality-filtered. SNP alleles were called using an alternative coverage-based approach where a missing call was assigned to a site if the total coverage was less than 20× depth or where one of the four nucleotides accounted for at least 80% total coverage [60]. There was a clear relationship between the mean depth of coverage and the proportion of missing calls, based on which we retained sequences with greater than 10× mean depth of coverage over the whole genome (81 sequences retained).

Heterozygous calls were removed, and SNPs with a minor allele frequency (MAF) of less than 25% were removed. Samples with greater than 25% genome-wide missing data and 30% missing data per SNP were excluded from the analysis (n = 10, 71 sequences retained). All SNP positions with a MAF greater than 20% were identified using BCFtools v0.1.19 (https://samtools.github.io/bcftools/). Sequences were excluded from the final GWAS if more than 300 such positions were found using methods described by Hadfield et al. [61]. The quality assessment and filtering process is shown in Fig. 2. Details of the WGS data are provided in Additional file 2: Figure S2.

Fig. 2
figure 2

Whole genome sequencing (WGS) quality filtering processes and threshold criteria for inclusion in analyses. Ct DNA detected using droplet digital PCR [45]. WGS data were obtained using SureSelect target enrichment [31] (or chlamydial cell culture) and Illumina paired-end sequencing. FastQC [53] was used to assess basic WGS quality. SNP alleles were called against reference strain Ct A/HAR-13 using an alternative coverage-based approach where a missing call was assigned to a site if the total coverage was less than 20× depth or where one of the four nucleotides accounted for at least 80% total coverage [60]. There was a clear relationship between the mean depth of coverage and genome-wide proportion of missing calls; therefore, only sequences with greater than 10× mean depth of coverage over the whole genome were retained using the GATK Best Practices threshold [56, 57]. Heterozygous calls were removed and SNPs with a minor allele frequency (MAF) of less than 25% were removed. Samples with greater than 25% genome-wide missing data and 30% missing data per SNP were excluded from the analysis. WGS sequence quality is shown in detail in Additional file 12: Figure S12. *n = 157 including the 71 Bijagós sequences in addition to 48 Rombo District sequences and 38 reference sequences

Full size image

Phylogenetic reconstruction

Samples were mapped to the ocular reference strain Ct A/HAR-13 and SNPs were called as described above. Phylogenies were computed using RAxML v7.8.2 [62] from a variable sites alignment using a generalized time-reversible (GTR) + gamma model and are midpoint rooted. Recombination is known to occur in Ct [4, 6] and can be problematic in constructing phylogeny. We applied three compatibility-based recombination detection methods to detect regions of recombination using PhiPack [63]: the pairwise homoplasy index (Phi), the maximum χ2 and the neighbour similarity score (NSS) across the genome alignment. We also examined the confidence in the phylogenetic tree by computing RAxML site-based likelihood scores [62]. Phylogenetic trees were examined adjusting for recombination using the methods described above.

Additionally, sequence data for the tryptophan operon (CTA0182 and CTA0184–CTA0186), tarP (CTA0498), nine polymorphic membrane proteins (CTA0447–CTA0449, CTA0884, CTA0949–CTA0952 and CTA0954) and ompA (CTA0742) were extracted from the 81 ocular Ct sequences from Guinea-Bissau retained after quality control filtering described above, 48 ocular sequences originating from a study conducted in Kahe village, Rombo District, Tanzania [64] and 38 publicly available reference sequences. Phylogenies were constructed as described above.

Polymorphisms, insertions and deletions (indels) and truncations for the tryptophan operon were manually determined from aligned sequences using SeaView [65]. Tyrosine repeat regions and actin-binding domains in tarP were found using RADAR [66] and Pfam [67] respectively.

Pairwise diversity

A comparison was made between the two population-based Ct sequence data sets from the Bijagós (Guinea-Bissau) and Rombo (Tanzania) sequences whereby short read data from the 81 Bijagós sequences and 48 Rombo sequences were mapped against Ct A/HAR-13 using SAMtools. Within-population pairwise nucleotide diversity was calculated using the formula:

$$ \pi =2\times {\varSigma}_{i=1}^n\ {\varSigma}_{j=1}^{i-1}{x}_i{x}_j{\pi}_{ij} $$

where n is the number of sequences, x is the frequency of sequences i and j and π ij is the number of nucleotide differences per site between sequences i and j [68]. The frequency of sequences was considered uniform within the populations, and sites with missing calls were excluded on a per-sequence basis.

Genome-wide association analyses

To investigate the association between Ct polymorphisms with ocular localization and clinical disease severity, we used permutation-based logistic regression methods, which are powerful and well-recognized tools in GWAS, allowing for adjustment for population structure, age and gender in the model and accounting for multiple testing [69,70,71,72].

We used permutation analyses of 100,024 phenotypic re-samplings, where the distribution of the p value was approximated by simulating data sets through randomization under the null hypothesis of no association between phenotype and genotype. Genome-wide significance was determined as p* ≤ 0.05, where p* was defined as the fraction of re-sampled (simulated) data that returned p values that were less than or equal to the p values observed in the data [50]. All analyses were conducted using the R statistical package v3.0.2 (the R Foundation for Statistical Computing, https://www.r-project.org/) using MASS, GLM and lsr. All R script used for these analyses is contained within Additional file 3: Figure S3 and is released as a CC-BY open resource (CC-BY-SA 3.0).

Ocular localization

Tissue localization is defined as the localization (or presence) of a detectable Ct infection to either the conjunctival epithelium or the urogenital tract. Short read data from the 129 clinical ocular sequences from the pairwise diversity analysis and 38 publicly available reference sequences from ocular (n = 8), urogenital (n = 17) and rectal (n = 13) sites were mapped against Ct A/HAR-13 using SAMtools. Only polymorphic sites were retained, and SNPs were filtered as described above. The final analysis includes 1007 SNPs from 157 sequences, a phylogeny of which is contained within Additional file 4: Figure S4. A permutation-based generalized linear regression model was used to test the association between collection site (ocular or urogenital tissue localization) and polymorphic sites. For each SNP the standard error for the t statistic was estimated from the model and used to calculate the odds ratios (ORs) and 95% confidence intervals. A χ2 test was used to determine the association between ocular localization-associated SNPs and both gene expression stage and predicted localization of the encoded proteins. The developmental cycle expression stage for each transcript was based on data and groupings from Belland et al. [73]. Predicted localization of expressed proteins was defined using the consensus from three predictions using CELLO [74], PSORTb [75] and LocTree3 [76].

Clinical disease severity

A permutation-based ordinal logistic regression model was used to test the association between the disease severity score (using the in vivo conjunctival phenotype defined previously) and polymorphic sites. The final analysis includes 129 SNPs from 71 sequences derived as described in Fig. 2. For each SNP the standard error for the t statistic was estimated from the model and used to calculate the ORs and 95% confidence intervals. Individuals’ age and gender were included as a covariate to the regression analysis.

We investigated the effect of population structure on the results of the GWAS analysis using PCA [77]. The first three PCs captured the majority of structural variation, but including them in the model had no effect; therefore, they were not included in the final model.

We corrected for genomic inflation if the occurrence of a polymorphism in the population was more than 90% or if there was a MAF of 3%.

Results

Conjunctival swabs collected during a cross-sectional population-based trachoma survey on the Bijagós Archipelago yielded 220 ocular Ct infections detected by Ct plasmid-based ddPCR. Of the 220 Ct infections detected, 184 were quantifiable using Ct genome-based ddPCR.

We obtained WGS data from 126/220 samples using cell culture (n = 8) or direct sequencing from swabs with SureSelectXT target enrichment (n = 118), representing the largest cross-sectional collection of ocular Ct WGS. Eighty-one of these sequences were subsequently included in the phylogenetic and diversity analyses and 71 were retained in the final genome-wide association (tissue localization (derived from the anatomical site of sample collection) and disease severity) analyses. The quality filtering process is illustrated in Fig. 2 and detailed in Methods.

A total of 1034 unique SNP sites were identified within the 126 Bijagós Ct genomes relative to the reference strain Ct A/HAR-13. Following application of further threshold criteria based on MAF and genome-wide missing data thresholds, we retained only high-quality genomic data in the final association analyses (129 SNPs from 71 sequences). There were no significant differences between the 71 retained and the 55 excluded sequences with respect to demographic characteristics, bacterial load, disease severity scores or geographical location (Table 1). Clinical and demographic details of the survey participants in whom we did not identify Ct infection have been published previously [43]. Of the ten SNPs initially identified within the Ct plasmid sequences, none fulfilled the quality filtering criteria, and they were not retained for the genome-wide association analyses.

Table 1 Characteristics of ocular Chlamydia trachomatis sequences included in the disease severity association analysis
Full size table

Ocular C. trachomatis phylogeny and diversity

For the phylogeny and diversity analyses, 81 Bijagós Ct sequences were included on the basis of the quality filtering criteria described in detail in Fig. 2. SNP-based phylogenetic trees constructed using all 1034 SNPs for sequences above 10× coverage (n = 81), with 54 published Ct reference genomes, are shown in Fig. 3.

Fig. 3
figure 3

Maximum likelihood reconstruction of whole genome phylogeny of ocular Chlamydia trachomatis sequences from the Bijagós Archipelago (Guinea-Bissau). Maximum likelihood reconstruction of the whole genome phylogeny of 81 Ct sequences from the Bijagós Islands and 54 Ct reference strains. Bijagós Ct sequences (n = 81) were mapped to Ct A/HAR-13 using SAMtools [55]. SNPs were called as described by Harris et al. [4]. Phylogenies were computed with RAxML [62] from a variable sites alignment using a GTR + gamma model and are midpoint rooted. The scale bar indicates evolutionary distance. Bijagós Ct sequences in this study are coloured black, and reference strains are coloured by tissue localization (red = Ocular, green = Urogenital, blue = LGV). Branches are supported by > 90% of 1000 bootstrap replicates. Branches supported by 80–90% (orange) and < 80% (brown) bootstrap replicates are indicated

Full size image

The Bijagós sequences are situated within the T2 ocular monophyletic lineage with all other ocular Ct sequences [59] except those described by Andersson et al. [10]. However, our population-based collection of ocular Ct sequences has much greater diversity at whole genome resolution than previously demonstrated in African trachoma isolates [4, 8]. We used a pairwise diversity (π) metric to compare two populations of ocular Ct from regions with similar trachoma endemicity and studies with similar design, sample size and available epidemiological metadata. These data show much greater genomic diversity in the Bijagós ocular Ct sequences (π = 0.07167) compared to the Tanzanian (Rombo) ocular Ct sequences (π = 0.00047).

By ompA genotyping, 73 of the Bijagós sequences are genotype A and 8 are genotype B, supporting their classical ocular nature (Additional file 5: Figure S5). The high resolution of WGS data obtained directly from clinical samples captures diversity that may be useful in strain classification, particularly as we found some evidence of clustering at village level, although the very small number of sequences per village means that it is not possible to provide accurate estimates of clustering in this study (Fig. 4).

Fig. 4
figure 4

Maximum likelihood phylogenetic tree showing clustering of ocular Chlamydia trachomatis sequence types by village. RAxML maximum likelihood phylogenetic reconstruction including all ocular Ct sequences retained in the final disease severity association analysis after quality filtering (n = 71). Ocular Ct sequences are labelled by village (villages numbered and coloured), midpoint-rooted and mapped to reference Ct A/HAR-13. Branches are supported by > 90% of 1000 bootstrap replicates. Branches supported by 80–90% (orange) and < 80% (brown) bootstrap replicates are indicated

Full size image

Homoplasic SNPs and regions affected by recombination are shown in Additional file 6: Figure S6a. Removal of these regions of recombination identified using the pairwise homoplasy index had no effect on phylogenetic relationships. Additionally, a site-wise log likelihood plot demonstrated that there was no clear genomic region where there was significant lack of confidence in the tree construction due to recombination (Additional file 6: Figure S6b). Whether regions containing recombination were included or excluded, tree topology remained essentially identical, indicating that branching order is not affected by the removal of these regions.

Genome-wide analysis of C. trachomatis localization

Candidate genes thought to be involved in or indicative of ocular localization or preference were examined to further characterize this population of ocular Ct. Polymorphisms and truncations in the tryptophan operon have previously been implicated in the inability of ocular Ct to infect and survive in the genital tract [5]. All sequences contained mutations in trpA resulting in truncation. The majority (80/81) were truncated at the previously characterized deletion at position 533 [5]. Polymorphisms in trpB and trpR were less common (Additional file 7: Figure S7).

The variable domain structure of the translocated actin-recruiting phosphoprotein (tarP) has also been implicated in tropism [78]. Ocular strains possess more actin-binding domains (three or four) and fewer tyrosine repeat regions (between one and three). Urogenital strain tarP sequences have low copy numbers of both, and LGV strain sequences have additional tyrosine repeat regions. In this study, all sequences contain the expected three tyrosine repeat regions and three or four actin-binding domains (Additional file 7: Figure S7).

The nine virulence-associated polymorphic membrane proteins (Pmp) are variably related to tissue preference, with all encoding genes except pmpA, pmpD and pmpE clustering by tissue location [20]. In this population all phylogenies of the six tropism-clustering pmps show that all sequences cluster with other ocular sequences (Additional file 8: Figure S8).

Permutation-based re-sampling methods, commonly used in GWAS analyses, were used to account for multiple comparisons [69,70,71,72]. We tested 1007 SNPs in 157 Ct sequences (Fig. 2) for association with ocular localization (defined by anatomical site of sample collection), comparing 127 ocular, 17 urogenital and 13 LGV strains (Fig. 5a). One hundred and five SNPs were significantly associated with ocular localization (p* < 0.05), of which 21 were non-synonymous (details in Table 2a and Additional file 9: Figure S9). These were within a number of genes known to be polymorphic, genes previously identified as tropism-associated (CTA0156, CTA0498/tarP and CTA0743/pbpB) and virulence factors (CTA0498/tarP and CTA0884/pmpD). Four genes contained multiple non-synonymous SNPs (CTA_0733/karG, CTA_089/5sucD, CTA_0087 and CTA_0145/oppA_1), and ten genes contained multiple synonymous SNPs. Of the genes containing multiple synonymous SNPs, five contained more than three SNPs (CTA_0739/tsf, CTA_0733/karG, CTA_0156, CTA_0154 and CTA_0153). No predicted protein localization was over-represented in the ocular localization-related SNPs (p = 0.6174); however, early and very-late expressed genes were over-represented (p = 0.0197).

Fig. 5
figure 5

Single nucleotide polymorphisms on the Chlamydia trachomatis genome associated with (a) ocular localization and (b) disease severity at genome-wide significance. a Ocular localization-associated SNPs across the C. trachomatis genome. There were 1007 SNPs identified in coding and non-coding regions and included in permutation-based linear regression models in the Ct genome-wide association analysis. The threshold for genome-wide significance is indicated by the dashed line (p* < 0.05). The y-axis shows the –log10 p value. A –log10 p value of 1.3 is equivalent to a permuted p value of 0.05 (p* < 0.05). Synonymous (black) and non-synonymous SNPs (red) are indicated. Regions informative for ocular localization and genes of interest are labelled in blue. b Disease severity-associated SNPs across the Ct genome. From 129 SNPs identified in coding and non-coding regions, SNPs associated with the disease severity phenotype at genome-wide significance are identified using permutation-based ordinal logistic regression models adjusting for age in the Ct genome-wide association analysis. The threshold for genome-wide significance is indicated by the dashed line (p* < 0.05). The y-axis shows the –log10 p value. A log10 p value of 1.3 is equivalent to a permuted p value of 0.05 (p* < 0.05). Genes significantly associated with disease severity are labelled in blue

Full size image
Table 2 SNPs across the Chlamydia trachomatis genome identified using permutation-based genome-wide association analysis for (A) ocular localization (non-synonymous only) and (B) disease severity
Full size table

Markers of disease severity in ocular C. trachomatis infection

Using permutation-based re-sampling methods, eight SNPs were found to be significantly associated with disease severity (Fig. 5b). Seven of these are in coding regions (relative to Ct A/HAR-13). Five are present at nucleotide positions 465,330 (OR = 0.13, p* = 0.037), 32,779 (OR = 0.12, p* = 0.032), 875,804 (OR = 0.10, p* = 0.024), 939,488 (OR = 0.10, p* = 0.026) and 1,028,728 (OR = 0.08, p* = 0.013) (where p* is the permuted p value with a genome-wide threshold of 0.05) representing synonymous codon changes within the genes yjfH, trmD, alaS, glgA and pmpE respectively. Three further genome-wide significant synonymous SNPs were present at positions 827,184 (OR = 0.3, p* = 0.041) within the predicted coding sequence (CDS) CTA0744, 285,610 (OR = 0.12, p* = 0.027) within CTA0273 and 787,841 (OR = 0.13, p* = 0.043) in the intergenic region between loci CTA0744–CTA0745 (Table 2b and Additional file 10: Figure S10).

Discussion

This collection of clinical ocular Ct WGS from a single trachoma-endemic population to be characterized has enabled us to describe the population diversity of naturally occurring Ct in a treatment-naïve population. We used detailed clinical grading combined with microbial quantitation to perform a GWAS and investigated associations between Ct polymorphisms with ocular localization and disease severity in trachoma.

Unlike the recently published Australian Ct sequences [10], all Bijagós sequences clustered as expected within the T2 ocular clade derived from a urogenital ancestor [59, 61], each with loci typically associated with ocular tissue localization (trpA and tarP). Although the Bijagós sequences conform to the classical ocular genotype, the phylogenetic data show greater than expected diversity compared to historical reference strains of ocular Ct [4] and a population of clinical ocular Ct sequences obtained from cultured clinical conjunctival swab specimens collected from another African trachoma-endemic population [64] (Additional file 4: Figure S4). Our use of direct WGS from clinical samples reveals the natural diversity of a population-based collection of endemic treatment-naïve ocular Ct infections. This diversity may indicate genome-wide selection for advantageous mutations as demonstrated in other pathogens [79] or simply the naturally diverse circulation of endemic treatment-naïve ocular Ct.

The apparent village-level clustering provides new evidence that WGS has the necessary molecular resolution to fully investigate Ct transmission. Although the number of sequences from each village was very small, overall Ct genomic diversity supports our hypothesis of ongoing or recent transmission, since diversity requires mutation, recombination and gene flow. The data from this study demonstrate such mutation and indicate that WGS data may be useful in defining transmission networks and developing transmission maps, which have not been adequately defined using alternative Ct genotyping systems. Whole genome mapping has previously been shown to be a useful tool in the analysis of outbreaks and bacterial pathogen transmission [80, 81] and thus has multiple potential applications in epidemiological analysis and transmission studies. However, greater numbers of sequences per village are required to validate this finding.

Such diversity is likely to be representative of recombination present in Ct [82]. Genome-wide recombination was common and widespread within these sequences. Extensive recombination has been noted in previous studies and is thought to be a source of diversification with possible interstrain recombination [4, 82]. Recombination may represent fixation of recombination in regions that are under diversifying selection pressure [4].

Recently, a handful of bacterial GWASs have provided insight into the genetic basis of bacterial host preference, antibiotic resistance and virulence [83,84,85,86,87,88]. Until now, most inferences regarding disease-modifying virulence factors in chlamydial infection have been derived from a limited number of comparative genomic studies where only a few virulence factors were associated with disease severity. Chlamydial genomic association data have previously been used to highlight genes potentially involved in pathoadaptation [10, 89] and tissue localization [90].

In the current GWAS we found 21 genome-wide significant non-synonymous SNPs associated with ocular localization and eight genome-wide significant synonymous SNPs associated with disease severity.

Confidence that new SNPs identified in the ocular localization GWAS are candidate markers of pathoadaptation is supported by the observation that half of the SNPs identified have previously been described as polymorphic or recombinant within Ct and the ocular serovars [8, 91,92,93].

In support of the hypothesis that early events in infection and intracellular growth are crucial events in Ct survival and pathogenicity, we identified SNPs within genes that are expressed from the beginning of the chlamydial developmental cycle including CTA0156 (encoding early endosomal antigen 1 (EEA1) [73]), CTA0498 (encoding translocated actin-recruiting phosphoprotein (tarP) [94]) and CTA0884 (encoding polymorphic membrane protein D (PmpD) [95]), which have identified roles in entry to and initial interactions with host cells.

Two of the four genes containing multiple non-synonymous SNPs (karG and sucD) are involved in ATP metabolism and, more generally, chlamydial metabolism. Two of the genes with multiple synonymous mutations (ruvB and CTA_0284) are also involved in metabolism. Growth rates are known to vary significantly between biovars. The developmental cycle in ocular serovars is substantially longer than that in genital serovars [96]. These genes and the identified SNPs may therefore be important in the differential growth and development of Ct serovars. This is supported by the downregulation of sucD expression during in vitro persistence. Slower growth in ocular strains occurs primarily in the entry and early stages of differentiation, which may also indicate the role of previously described genes involved in entry into cells.

The eight disease severity-associated SNPs are within less well-characterized genes. Apart from pmpE, the remaining genes identified in this study have been shown to be relatively conserved [90]. This suggests that these SNPs may be important in ocular Ct pathogenesis, rather than in longer term chlamydial evolution. Three of these genes are putative Ct virulence factors, with functions in nutrient acquisition (glgA [24, 28, 97]), host-cell adhesion (pmpE [98]) and response to IFNγ-induced stress (trmD [73]). Homologues of alaS [99, 100] and CTA0273 [101, 102] are known virulence factors in related Gram-negative bacteria, suggesting that these genes are potentially important in Ct pathogenesis.

Transcriptome analysis of chlamydial growth in vitro has shown that there is highly upregulated gene expression of trmD (encoding a transfer RNA (tRNA) methyltransferase) associated with growth in the presence of IFNγ, thought to be important in the maintenance of chlamydial infection [73]. yjfH (renamed rlmB) is phylogenetically related to the TrmD family and encodes the protein RlmB, which is important for the synthesis and assembly of the components of the ribosome [103]. In Escherichia coli, Haemophilus influenzae and Mycoplasma genitalium, RlmB catalyses the methylation of guanosine 2251 in 23S ribosomal RNA (rRNA), which is of importance in peptidyl tRNA recognition but is not essential for bacterial growth [103, 104]. alaS encodes a tRNA ligase of the class II aminoacyl tRNA synthetase family involved in cytoplasmic protein biosynthesis. It is not known to have virulence associations in chlamydial infection, but has been described as a component of a virulence operon in Haemophilus ducreyi [99] and H. influenzae [100]. The CDS CTA0273 encodes a predicted inner membrane protein translocase component of the autotransporter YidC, an inner membrane insertase important in virulence in E. coli [101] and Streptococcus mutans [102]. Our study suggests that these loci may be important in disease severity and host-pathogen interactions in chlamydial infection. A summary of available literature for these key ocular localization and disease severity-associated SNPs is tabulated in Additional file 11: Figure S11. We cannot speculate further on the effect of these polymorphisms on expression. It is possible that the synonymous disease severity-associated SNPs are markers in linkage for disease-causing alleles that were not included in the final GWAS analysis. For both analyses, further mechanistic studies are required to establish causality and validity and to fully understand the nature of the associations presented.

Though we were intrinsically limited to those cases where infection was detectable and from which we were able to obtain Ct WGS data, our population-based treatment-naïve sample attempts to provide a representative picture of what is observed in ocular Ct infection. We acknowledge that there may be Ct genotypes that are cleared by the immune system such that we do not capture them in a cross-sectional study. We are limited to the small sample size in this study, but attempt to address the issues of statistical power and multiple testing by using a bi-dimensional conjunctival phenotype and permutation-based multivariable regression analysis. To date, many published microbial GWASs have sample sizes under 500 [105], including several key studies examining virulence [84] and drug resistance [85] in Staphylococcus aureus with sample sizes of 75 and 90 respectively.

Conclusions

The potential of bacterial GWASs has only recently been realized, and despite the limitations with sample size, their use to study Ct in this way is particularly important, since in vitro models are intrinsically difficult to develop, and it has not been possible to study urogenital Ct in the same way due to the lack of a clearly defined in vivo disease phenotype. The genomic markers identified in this study provide important direction for validation through in vitro functional studies and a unique opportunity to understand host-pathogen interactions likely to be important in Ct pathogenesis in humans. The greater than expected diversity within this population of naturally circulating ocular Ct and the clustering at village level demonstrate the potential utility of WGS in epidemiological and clinical studies. This will enable us to understand transmission in both ocular and urogenital Ct infection and will have significant public health implications in preventing and eliminating chlamydial disease in humans.

Abbreviations

ATP:

Adenosine triphosphate

Ct :

Chlamydia trachomatis

ddPCR:

Droplet digital PCR

DMEM:

Dulbecco’s modified Eagle’s medium

DNA:

Deoxyribonucleic acid

EB:

Elementary body

FPC:

Follicles, papillary hypertrophy, conjunctival scarring

GWAS:

Genome-wide association study

indels:

Insertions and deletions

LGV:

Lymphogranuloma venereum

MAF:

Minor allele frequency

MOMP:

Major outer membrane protein

NHP:

Non-human primate

NSS:

Neighbour similarity score

PC:

Principal component

PCA:

Principal component analysis

PCR:

Polymerase chain reaction

SNP:

Single nucleotide polymorphism

WGS:

Whole genome sequencing

References

  1. Hu VH, et al. Epidemiology and control of trachoma: systematic review. Tropical Med Int Health. 2010;15(6):673–91.

    Article  Google Scholar 

  2. World Health Organization. Sexually transmitted infections (STIs); 2016. http://www.who.int/mediacentre/factsheets/fs110/en/.

    Google Scholar 

  3. Rodriguez P, et al. Typing of Chlamydia trachomatis by restriction endonuclease analysis of the amplified major outer membrane protein gene. J Clin Microbiol. 1991;29:1132–6.

    CAS  PubMed  PubMed Central  Google Scholar 

  4. Harris SR, et al. Whole genome analysis of diverse Chlamydia trachomatis strains identifies phylogenetic relationships masked by current clinical typing. Nat Genet. 2012;44(4):413–s1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Caldwell HD, et al. Polymorphisms in Chlamydia trachomatis tryptophan synthase genes differentiates between genital and ocular isolates: implications in pathogenesis and infection tropism. J Clin Invest. 2003;111:1757–69.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Jeffrey BM, et al. Genome sequencing of recent clinical Chlamydia trachomatis strains identifies loci associated with tissue tropism and regions of apparent recombination. Infect Immun. 2010;78:2544–53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Nunes A, Borrego MJ, Gomes JP. Genomic features beyond C. trachomatis phenotypes: what do we think we know? Infect Genet Evol. 2013;16:392–400.

    Article  CAS  PubMed  Google Scholar 

  8. Kari L, et al. Pathogenic diversity among Chlamydia trachomatis ocular strains in non-human primates is affected by subtle genomic variations. J Infect Dis. 2008;197:449–56.

    Article  CAS  PubMed  Google Scholar 

  9. Butcher RMR, et al. Low prevalence of conjunctival infection with Chlamydia trachomatis in a treatment-naive trachoma-endemic region of the Solomon Islands. PLoS Negl Trop Dis. 2016;10(10):e0005051.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Andersson P, et al. Chlamydia trachomatis from Australian Aboriginal people with trachoma are polyphyletic composed of multiple distinctive lineages. Nat Commun. 2016;7:10688.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Feng L, et al. Survey, culture and genome analysis of ocular Chlamydia trachomatis in Tibetan boarding primary schools in Qinghai Province, China. Front Cell Infect Microbiol. 2017;6:207.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Borges V, et al. Complete genome sequence of Chlamydia trachomatis ocular serovar C strain TW-3. Genome Announc. 2014;2:e01204–13.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Darville T, Hiltke T. Pathogenesis of genital tract disease due to Chlamydia trachomatis. J Infect Dis. 2010;201(Supplement_2):S114–225.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. van Valkengoed IG, et al. Overestimation of complication rates in evaluations of Chlamydia trachomatis screening programmes—implications for cost-effectiveness analyses. Int J Epidemiol. 2004;33(2):416–25.

    Article  PubMed  Google Scholar 

  15. Bailey RL, et al. Molecular epidemiology of trachoma in a Gambian village. Br J Ophthalmol. 1994;78:813–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Andreasen AA, et al. Chlamydia trachomatis ompA variants in trachoma: what do they tell us? PLoS Negl Trop Dis. 2008;2(9):e306.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  17. Wang Y, et al. Development of a transformation system for Chlamydia trachomatis: restoration of glycogen biosynthesis by acquisition of a plasmid shuttle vector. PLoS Pathog. 2011;7(9):e1002258.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Wang Y, et al. Transformation of a plasmid-free genital tract isolate with a plasmid vector carrying a deletion in CDS6 revealed that this gene regulates inclusion phenotype. Pathogens Dis. 2013;67(2):100–3.

    Article  CAS  PubMed Central  Google Scholar 

  19. Longbottom D, et al. Molecular cloning and characterization of genes coding for the highly immunogenic cluster of 90 kilodalton envelope proteins from the Chlamydia psittaci subtype that causes abortion in sheep. Infect Immun. 1998;66:1317–24.

    CAS  PubMed  PubMed Central  Google Scholar 

  20. Gomes JP, et al. Polymorphisms in the nine polymorphic membrane proteins of Chlamydia trachomatis across all serovars: evidence for serovar Da recombination and correlation with tissue tropism. J Bacteriol. 2006;188:275–86.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Rockey DD, Heinzen RA, Hackstadt T. Cloning and characterization of a Chlamydia psittaci gene coding for a protein localized in the inclusion membrane of infected cells. Mol Microbiol. 1995;15:617–26.

    Article  CAS  PubMed  Google Scholar 

  22. Hefty PS, Stephens RS. Chlamydia trachomatis type III secretion system is encoded on ten operons preceded by a sigma 70-like promoter element. J Bacteriol. 2007;189:198–206.

    Article  CAS  PubMed  Google Scholar 

  23. Carlson JH, et al. In vivo and in vitro studies of Chlamydia trachomatis TrpR:DNA interaction. Mol Microbiol. 2006;59(6):1678–91.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. O’Connell CM, et al. Toll-like receptor 2 activation by Chlamydia trachomatis is plasmid dependent, and plasmid-responsive chromosomal loci are coordinately regulated in response to glucose limitation by C. trachomatis but not by C. muridarum. Infect Immun. 2011;79:1044–56.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  25. Hackstadt T, Scidmore-Carlson MA, Shaw EI, Fischer ER. Chlamydia trachomatis IncA protein is required for homotypic vesicle fusion. Cell Microbiol. 1999;1:119–30.

    Article  CAS  PubMed  Google Scholar 

  26. Nelson DE, et al. Inhibition of Chlamydiae by primary alcohols correlates with strain specific complement of plasticity zone phospholipase D genes. Infect Immun. 2006;74(1):73–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Carlson JH, Hughes S, Hogan D. Polymorphisms in the Chlamydia trachomatis cytotoxin locus associated with ocular and genital isolates. Infect Immun. 2004;72(12):7063–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Carlson JH, et al. The Chlamydia trachomatis plasmid is a transcriptional regulator of chromosomal genes and a virulence factor. Infect Immun. 2008;76:2273–83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Frazer LC, et al. Plasmid-cured Chlamydia caviae activates TLR2-dependent signaling and retains virulence in the guinea pig model of genital tract infection. PLoS One. 2012;7(1):e30747.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Song L, et al. Chlamydia trachomatis plasmid-encoded pgp4 is a transcriptional regulator of virulence-associated genes. Infect Immun. 2013;81(3):636.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Christiansen MT, et al. Whole genome enrichment and sequencing of Chlamydia trachomatis directly from clinical samples. BMC Infect Dis. 2014;14:591.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  32. Borges V, et al. Effect of long-term laboratory propagation on Chlamydia trachomatis genome dynamics. Infect Genet Evol. 2013;17:23–32.

    Article  CAS  PubMed  Google Scholar 

  33. Borges V, et al. Chlamydia trachomatis in vivo to in vitro transition reveals mechanisms of phase variation and downregulation of virulence factors. PLoS One. 2015;10(7):e0133420.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  34. Bonner C, et al. Chlamydia trachomatis virulence factor CT135 is stable in vivo but highly polymorphic in vitro. Pathog Dis. 2015;73(6):ftv043.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  35. Burton MJ, et al. Which members of a community need antibiotics to control trachoma? Conjunctival Chlamydia trachomatis load in Gambian villages. Invest Ophthalmol Vis Sci. 2003;44(10):4215–22.

    Article  PubMed  Google Scholar 

  36. Conway DJ, et al. Scarring trachoma is associated with polymorphisms in TNF-alpha gene promoter and with increased TNF-alpha levels in tear fluid. Infect Immun. 1997;65(3):1003–6.

    CAS  PubMed  PubMed Central  Google Scholar 

  37. West SK, et al. Progression of active trachoma to scarring in a cohort of Tanzanian children. Ophthalmic Epidemiol. 2001;8(2–3):137–44.

    Article  CAS  PubMed  Google Scholar 

  38. Burton MJ, et al. Pathogenesis of progressive scarring trachoma in Ethiopia and Tanzania: two cohort studies. PLoS Negl Trop Dis. 2015;9(5):e0003763.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Dawson CR, Jones BR, Tarizzo ML. Guide to trachoma control in programs for the prevention of blindness. Geneva: World Health Organization; 1981.

    Google Scholar 

  40. Reid JS, et al. A principal component meta-analysis on multiple anthropometric traits identifies novel loci for body shape. Nat Comms. 2016;7:13357.

    Article  CAS  Google Scholar 

  41. Yang J, et al. Conditional and joint multiple SNP analysis of GWAS summary statistics identifies additional variance influencing complex traits. Nat Genet. 2012;44:S1–3.

    Article  CAS  Google Scholar 

  42. Aschard H, et al. Maximising the power of principal component analysis of correlated phenotypes in genome wide association studies. Am J Hum Genet. 2014;94:662–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Last AR, et al. Risk factors for active trachoma and ocular Chlamydia trachomatis infection in treatment-naïve trachoma-hyperendemic communities of the Bijagós Archipelago, Guinea Bissau. PLoS Negl Trop Dis. 2014;8(6):e2900.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Last A, et al. Plasmid copy number and disease severity in naturally occurring ocular Chlamydia trachomatis infection. J Clin Microbiol. 2014;52(1):324–7.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  45. Roberts C, et al. Development and evaluation of a next generation digital PCR diagnostic assay for ocular Chlamydia trachomatis infections. J Clin Microbiol. 2013;51(7):2195–203.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Burton MJ, et al. Conjunctival chlamydial 16S ribosomal RNA expression in trachoma: is chlamydial metabolic activity required for disease to develop? Clin Infect Dis. 2006;42:463–70.

    Article  CAS  PubMed  Google Scholar 

  47. Faal N, et al. Conjunctival FOXP3 expression in trachoma: do regulatory T cells have a role in human ocular Chlamydia trachomatis infection? PLoS Med. 2006;3(8):e266.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  48. Derrick T, et al. Inverse relationship between micro RNA-155 and -184 expression with increasing conjunctival inflammation during ocular Chlamydia trachomatis infection. BMC Infect Dis. 2016;16:60.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  49. Last A, et al. Spatial clustering of high load ocular Chlamydia trachomatis infection in trachoma: a cross-sectional population-based study. Pathog Dis. 2017;75(5) https://doi.org/10.1093/femspd/ftx050.

  50. Suo C, et al. Analysis of multiple phenotypes in genome wide genetic mapping studies. BMC Bioinformatics. 2013;14:151.

    Article  PubMed  PubMed Central  Google Scholar 

  51. Skipp P, Robinson J, O’Connor CD, Clarke IN. Shotgun proteomic analysis of Chlamydia trachomatis. Proteomics. 2005;5(6):1558–73.

    Article  CAS  PubMed  Google Scholar 

  52. Seth-Smith HMB, et al. Co-evolution of genomes and plasmids within Chlamydia trachomatis and the emergence in Sweden of a new variant strain. BMC Genomics. 2009;10:239.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  53. Andrews S. FastQC: a quality control tool for high throughput sequence data; 2010. http://www.bioinformatics.babraham.ac.uk/projects/fastqc.

    Google Scholar 

  54. Langmead B, et al. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25.3.

    Google Scholar 

  55. Li H, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  56. McKenna A, et al. The Genome Analysis Tool Kit: a MapReduce framework for analysing next generation DNA sequencing data. Genome Res. 2010;20:1297–303.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Cheng AY, Teo YY, Ong RT. Assessing single nucleotide variant detection and genotype calling on whole genome sequenced individuals. Bioinformatics. 2014;30(12):1707–13.

    Article  CAS  PubMed  Google Scholar 

  58. Gudbjartsson DF, et al. Large-scale whole genome sequencing of the Icelandic population. Nat Genet. 2015;47(5):435–44.

    Article  CAS  PubMed  Google Scholar 

  59. Pickering H, et al. Genome-wide profiling of humoral immunity and pathogen genes under selection identifies immune evasion tactics of Chlamydia trachomatis during ocular infection. Sci Rep. 2017;7(1):9634.

    Article  PubMed  PubMed Central  Google Scholar 

  60. Coll F, et al. PolyTB: a genomic variation map for Mycobacterium tuberculosis. Tuberculosis (Edinb). 2014;94(3):346–54.

    Article  CAS  Google Scholar 

  61. Hadfield J, et al. Comprehensive global genome dynamics of Chlamydia trachomatis show ancient diversification followed by contemporary mixing and recent lineage expansion. Genome Res. 2017;27:1–10.

    Article  CAS  Google Scholar 

  62. Stamatakis A, Hoover P, Rougemont J. A rapid bootstrap algorithm for the RAxML Web servers. Syst Biol. 2008;57:758–71.

    Article  PubMed  Google Scholar 

  63. Bruen TC, Philippe H, Bryant D. A simple and robust statistical test for detecting the presence of recombination. Genetics. 2006;172:2665–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Solomon AW, et al. Strategies for control of trachoma: observational study with quantitative PCR. Lancet. 2003;362(9379):198–204.

    Article  PubMed  Google Scholar 

  65. Gouy M, Guindon S, Gascuel O. SeaView version 4: a multiplatform graphical user interface for sequence alignment and phylogenetic tree building. Mol Biol Evol. 2010;27(2):221–114.

    Article  CAS  PubMed  Google Scholar 

  66. Heger A, Holm L. Rapid automatic detection and alignment of repeats in protein sequences. Proteins. 2000;41(2):224–37.

    Article  CAS  PubMed  Google Scholar 

  67. Sonnhammer EL, Eddy SR, Durbin R. Pfam: a comprehensive database of protein domain families based on seed alignments. Proteins. 1997;28(3):405–20.

    Article  CAS  PubMed  Google Scholar 

  68. Nei M, Masatoshi N, Wen-Hsiung L. Mathematical model for studying genetic variation in terms of restriction endonucleases. Proc Natl Acad Sci. 1979;76(10):5269–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Che R, et al. An adaptive permutation approach for genome-wide association study. BioData Mining. 2014;7:9.

    Article  PubMed  PubMed Central  Google Scholar 

  70. Dudoit S, et al. Multiple hypothesis testing in microarray experiments. Stat Sci. 2003;18(1):71–103.

    Article  Google Scholar 

  71. Sham PC, Purcell SM. Statistical power and significance testing in large scale genetic studies. Nat Rev Genet. 2014;15(5):335–46.

    Article  CAS  PubMed  Google Scholar 

  72. Fairfax BP, et al. Genetics of gene expression in primary immune cells identifies cell type-specific master regulators and roles of HLA alleles. Nat Genet. 2012;44(5):502–10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Belland RJ, et al. Genomic transcriptional profiling of the developmental cycle of Chlamydia trachomatis. Proc Natl Acad Sci U S A. 2003;100(14):8478–83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Yu CS, Lin CJ, Hwang JK. Predicting subcellular localization of proteins for Gram-negative bacterial by support vector machines based on n-peptide compositions. Protein Sci. 2004;13(5):1402–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Yu NY, et al. PSORTb 3.0: improved protein subcellular localization prediction with refined localization subcategories and predictive capabilities for all prokaryotes. Bioinformatics. 2010;26(13):1608–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Goldberg T, et al. LocTree3 prediction of localization. Nucleic Acids Res. 2014;42(Web Server Issue):W350–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Thomas A, et al. Effect of linkage disequilibrium on the identification of functional variants. Genet Epidemiol. 2011;35(Suppl 1):S115–9.

    Article  PubMed  PubMed Central  Google Scholar 

  78. Lutter EI, et al. Phylogenetic analysis of Chlamydia trachomatis TARP and correlation with clinical phenotype. Infect Immun. 2010;78(9):3678–88.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Pepperell CS, et al. The role of selection in shaping diversity of natural Mycobacterium tuberculosis populations. PLoS Pathog. 2013;9(8):10.1371.

    Article  Google Scholar 

  80. Fey PD, et al. Assessment of whole genome mapping in a well-defined outbreak of Salmonella enterica serotype Saintpaul. J Clin Microbiol. 2012;50(9):3063–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Gilchrist CA, Turner SD, Riley MF, Petri WA Jr, Hewlett EL. Whole genome sequencing in outbreak analysis. Clin Microbiol Rev. 2015;28(3):541–63.

    Article  PubMed  PubMed Central  Google Scholar 

  82. Gomes JP, et al. Evolution of Chlamydia trachomatis diversity occurs by widespread interstrain recombination involving hotspots. Genome Res. 2007;17(1):50–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  83. Sheppard SK, et al. Genome-wide association study identifies vitamin B5 biosynthesis as a host specificity factor in Campylobacter. Proc Natl Acad Sci U S A. 2013;110:11923–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Laabei M, et al. Predicting the virulence of MRSA from its genome sequence. Genome Res. 2014;24(5):839–49.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  85. Alam MT, et al. Dissecting vancomycin intermediate resistance in Staphylococcus aureus using genome-wide association. Genome Biol Evol. 2014;6:1174–85.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  86. Chewapreecha C, et al. Comprehensive identification of single nucleotide polymorphisms associated with beta-lactam resistance within pneumococcal mosaic genes. PLoS Genet. 2014;10:e1004547.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  87. Maury M, et al. Uncovering Listeria monocytogenes hypervirulence by harnessing its biodiversity. Nat Genet. 2016;48:308–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  88. Earle SG, et al. Identifying lineage effects when controlling for population structure improves power in bacterial association studies. Nature Microbiol. 2016;1:16041.

    Article  CAS  Google Scholar 

  89. Borges V, Gomes JP. Deep comparative genomics among Chlamydia trachomatis lymphogranuloma venereum isolates highlights genes potentially involved in pathoadaptation. Infect Genet Evol. 2015;32:74–88.

    Article  CAS  PubMed  Google Scholar 

  90. Ferreira R, et al. In silico scrutiny of genes revealing phylogenetic congruence with clinical prevalence or tropism properties of Chlamydia trachomatis strains. G3. 2014;5(1):9–19.

    Article  PubMed  PubMed Central  Google Scholar 

  91. Carlson JH, Porcella SF, McClarty G. Comparative genomic analysis of Chlamydia trachomatis oculotropic and geniotropic strains. Infect Immun. 2005;73(10):6407–18.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  92. Brunelle BW, Nicholson TL, Stephens RS. Micro-array based genomic surveying of gene polymorphisms in Chlamydia trachomatis. Genome Biol. 2004;5(6):R42.

    Article  PubMed  PubMed Central  Google Scholar 

  93. Borges V, Nunes A, Ferreira R, Borrego MJ, Gomes JP. Directional evolution of Chlamydia trachomatis towards niche-specific adaptation. J Bacteriol. 2012;194(22):6143–53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  94. Clifton DR, et al. A Chlamydial type III translocated protein is tyrosine-phosphorylated at the site of entry and associated with the recruitment of actin. Proc Natl Acad Sci U S A. 2004;101(27):10166–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  95. Kari L, et al. Chlamydia trachomatis polymorphic membrane protein D is a virulence factor involved in early host-cell interactions. Infect Immun. 2014;82(7):2756–62.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  96. Miyairi I, et al. Different growth rates of Chlamydia trachomatis biovars reflect pathotype. J Infect Dis. 2006;194(3):350–7.

    Article  PubMed  Google Scholar 

  97. Da Cunha M, et al. Identification of type THREE secretion substrates of Chlamydia trachomatis using Yersinia enterocolitica as a heterologous system. BMC Microbiol. 2014;17(14):40.

    Article  CAS  Google Scholar 

  98. Becker E, Hegemann JH. All subtypes of the Pmp adhesion family are implicated in chlamydial virulence and show species-specific function. Microbiology. 2014;3(4):544–56.

    CAS  Google Scholar 

  99. Gangaiah D, et al. Carbon storage regulator A contributes to the virulence of Haemophilus ducreyi in humans by multiple mechanisms. Infect Immun. 2013;81(2):608–17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  100. Hogg JS, et al. Characterisation and modeling of the Haemophilus influenzae core and supragenomes based on the complete genomic sequences of Rd and 12 clinical non-typeable strains. Genome Biol. 2007;8:R103.1–8.

    Article  CAS  Google Scholar 

  101. Samuelson JC, et al. YidC mediates membrane protein insertion in bacteria. Nature. 2000;406(6796):637–41.

    Article  CAS  PubMed  Google Scholar 

  102. Palmer SR, et al. YidC1 and YidC2 are functionally distinct proteins involved in protein secretion, biofilm formation and cariogenicity of Streptococcus mutans. Microbiology. 2012;158(7):1702–12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  103. Michel G, et al. The structure of the RlmB 23S rRNA methyltransferase reveals a new methylation fold with a unique knot. Structure. 2002;10(10):1303–15.

    Article  CAS  PubMed  Google Scholar 

  104. Lovgren JM, Wikstrom PM. The rlmB gene is essential for formation of Gm2251 in 23S rRNA but not for ribosome maturation in Escherichia coli. J Bacteriol. 2001;183(23):6957–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  105. Power RA, Parkhill J, de Oliveira T. Microbial genome-wide association studies: lessons from human GWAS. Nat Rev Genet. 2017;18(1):41–50.

    Article  CAS  PubMed  Google Scholar 

  106. Le Negrate G, et al. ChaDub1 of Chlamydia trachomatis suppresses NF-KappaB activation and inhibits 1-KappaB-alpha ubiquitination and degradation. Cell Microbiol. 2008;10(9):1879–92.

    Article  CAS  PubMed  Google Scholar 

  107. Misaghi S, et al. Chlamydia trachomatis-derived deubiquinating enzymes in mammalian cells during infection. Mol Microbiol. 2006;61(1):142–50.

    Article  CAS  PubMed  Google Scholar 

  108. Li Z, et al. Induction of protective immunity against Chlamydia muridarum intravaginal infection with a chlamydial glycogen phosphorylase. PLoS One. 2012;7(3):e32997.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  109. Lu C, et al. Chlamydia trachomatis GlgA is secreted into host cell cytoplasm. PLoS One. 2013;8(7):e68764.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  110. Swanson KA, et al. Chlamydia trachomatis polymorphic membrane protein D is an oligomeric autotransporter with a higher-order structure. Infect Immun. 2009;77:508–16.

    Article  CAS  PubMed  Google Scholar 

  111. Werhl W, et al. From the inside out — processing of the chlamydial autotransporter Pmp D and its role in bacterial adhesion and activation of human host cells. Mol Micribiol. 2004;51(2):319–34.

    Article  CAS  Google Scholar 

  112. Nunes A, Gomes JP, Karunakaran KP, Brunham RC. Bioinformatic analysis of Chlamydia trachomatis PmpE, PmpF, PmpG and PmpH as potential vaccine antigens. PLoS One. 2015;10(7):e0131695.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  113. Tan C, Spitznagel JK, Shou HZ, et al. The polymorphic membrane protein gene family of the chlamydiaciae. In: Bavoil PM, Wyrick PB, editors. Chlamydia: genomics and pathogenesis. Norfolk: Horizon Bioscience; 2006. p. 195–218.

  114. Liu X, et al. Identification of Chlamydia trachomatis outer membrane complex proteins by differential proteomics. J Bacteriol. 2010;192(11):2852–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We extend thanks to colleagues at the Programa Nacional de Saúde de Visão in the Ministério de Saúde Publica in Bissau, the study participants and dedicated field research team in Guinea-Bissau and the Medical Research Council (MRC) Unit The Gambia for their support and collaboration in this work.

Funding

ARL was funded by the Wellcome Trust through a Clinical Research Training Fellowship (grant number 097330/Z/11/Z). MJH and ChR were funded by a Wellcome Trust Program Grant (grant number 079246/Z/06/Z). ChR was funded by the Wellcome Trust Institutional Strategic Support Fund (grant number 105609/Z/14/Z). Work undertaken by the Wellcome Trust Sanger Institute by NRT, HSS and JH was supported by the Wellcome Trust (grant number 098051) and the Wellcome Trust Institutional Strategic Support Fund (grant number 105609/Z/14/Z). TGC is funded by the MRC UK (grant numbers MR/K000551/1, MR/M01360X/1, MR/N010469/1). JP was funded by a Biotechnology and Biological Sciences Research Council (BBSRC) PhD studentship, and FC and HP were funded by Bloomsbury Colleges Research Fund PhD studentships.

Availability of data and materials

All sequence data are available from the European Bioinformatics Institute (EBI) short read archive. See Additional file 12: Figure S12 for details and accession numbers.

Author information

Authors and Affiliations

  1. Clinical Research Department, London School of Hygiene and Tropical Medicine, Keppel Street, London, UK

    A. R. Last, H. Pickering, C. h. Roberts, S. E. Burr, D. C. W. Mabey, R. L. Bailey & M. J. Holland

  2. Department of Pathogen Molecular Biology, London School of Hygiene and Tropical Medicine, Keppel Street, London, UK

    F. Coll, J. Phelan, T. G. Clark & N. R. Thomson

  3. Disease Control and Elimination Theme, Medical Research Council Unit The Gambia, Fajara, Gambia

    S. E. Burr

  4. Programa Nacional de Saúde de Visão, Ministério de Saúde Publica, Bissau, Guinea-Bissau

    E. Cassama & M. Nabicassa

  5. Pathogen Genomics, Wellcome Trust Sanger Institute, Wellcome Trust Genome Campus, Hinxton, UK

    H. M. B. Seth-Smith, J. Hadfield & N. R. Thomson

  6. Clinical Microbiology, Universitätsspital Basel, Basel, Switzerland

    H. M. B. Seth-Smith

  7. Applied Microbiology Research, Department of Biomedicine, University of Basel, Basel, Switzerland

    H. M. B. Seth-Smith

  8. Molecular Microbiology Group, University of Southampton Medical School, Southampton, UK

    L. T. Cutcliffe & I. N. Clarke

  9. Department of Infectious Diseases Epidemiology, London School of Hygiene and Tropical Medicine, Keppel Street, London, UK

    T. G. Clark

Authors
  1. A. R. Last
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. H. Pickering
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. C. h. Roberts
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. F. Coll
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. J. Phelan
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. S. E. Burr
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. E. Cassama
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. M. Nabicassa
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. H. M. B. Seth-Smith
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. J. Hadfield
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. L. T. Cutcliffe
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. I. N. Clarke
    View author publications

    You can also search for this author in PubMed Google Scholar

  13. D. C. W. Mabey
    View author publications

    You can also search for this author in PubMed Google Scholar

  14. R. L. Bailey
    View author publications

    You can also search for this author in PubMed Google Scholar

  15. T. G. Clark
    View author publications

    You can also search for this author in PubMed Google Scholar

  16. N. R. Thomson
    View author publications

    You can also search for this author in PubMed Google Scholar

  17. M. J. Holland
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

ARL, RLB, MJH, SEB and NRT designed the study. ARL, SEB, EC and MN conducted the field study. ARL, ChR and SEB conducted the molecular laboratory work. LTC and INC performed the chlamydial cell culture. HSS and JH designed and performed the whole genome sequencing and initial FastQC. ARL, ChR, HP, FC and TGC conducted the GWAS analysis. HP, JP, SH, JH and HSS supported the phylogenetic analysis. ARL, HP, MJH, DCWM, TGC and NRT wrote the paper. All authors have contributed to and reviewed the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to A. R. Last.

Ethics declarations

Ethics approval and consent to participate

This study was conducted in accordance with the declaration of Helsinki. Ethical approval was obtained from the Comitê Nacional de Ética e Saúde (Guinea-Bissau), the London School of Hygiene and Tropical Medicine Ethics Committee (UK) and The Gambia Government/MRC Joint Ethics Committee (The Gambia). Written informed consent to participate and publish anonymized patient data was obtained from all study participants or their guardians on their behalf if participants were children. A signature or thumbprint was considered an appropriate record of consent in this setting by the above ethical bodies. All communities received treatment for endemic trachoma in accordance with the World Health Organization and national policies following the survey.

Consent for publication

Written informed consent to publish anonymized patient data was obtained from all study participants as described above.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Figure S1. Histogram and density plot showing log-transformed C. trachomatis load (omcB copies/swab) data. (PDF 111 kb)

Additional file 2:

Figure S2. Detailed summary of whole genome sequence (WGS) data quality control of Bijagós Chlamydia trachomatis sequences. (PDF 76 kb)

Additional file 3:

Figure S3. R Script used for (A) tissue localization and (B) disease severity Chlamydia trachomatis GWAS. (PDF 180 kb)

Additional file 4:

Figure S4. Maximum likelihood reconstruction of whole genome phylogeny of Chlamydia trachomatis sequences examined in the tissue localization analysis. (PDF 357 kb)

Additional file 5:

Figure S5. Maximum likelihood reconstruction of the ompA (CTA0742) phylogeny. (PDF 450 kb)

Additional file 6:

Figure S6. Recombination present across Bijagós Chlamydia trachomatis genome sequences using the pairwise homoplasy index (Phi) and the site-wise log likelihood support for the best-scoring maximum likelihood tree. (PDF 153 kb)

Additional file 7:

Figure S7. Tyrosine repeat regions and actin-binding domains in tarP (CTA0948) and polymorphisms in the trp operon (CTA0182–CTA0186) (trpR, trpB and trpA) within Bijagós (Bissau-Guinean) ocular Chlamydia trachomatis sequences. (PDF 49 kb)

Additional file 8:

Figure S8. Maximum likelihood reconstruction of phylogeny by polymorphic membrane protein (Pmp) genes A–I. (PDF 1738 kb)

Additional file 9:

Figure S9. Ocular localization-associated SNPs (p value < 0.1). (PDF 150 kb)

Additional file 10:

Figure S10. SNPs across the Chlamydia trachomatis genome associated with disease severity using permutation-based genome-wide association analysis. (PDF 158 kb)

Additional file 11:

Figure S11. Summary of published studies supporting the key ocular localization and disease severity-associated SNPs [106,107,108,109,110,111,112,113,114]. (PDF 105 kb)

Additional file 12:

Figure S12. European Nucleotide Archive (ENA) (European Bioinformatics Institute (EBI)) accession numbers relating to C. trachomatis sequence data analysed in this study. (PDF 75 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Last, A.R., Pickering, H., Roberts, C.h. et al. Population-based analysis of ocular Chlamydia trachomatis in trachoma-endemic West African communities identifies genomic markers of disease severity. Genome Med 10, 15 (2018). https://doi.org/10.1186/s13073-018-0521-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13073-018-0521-x

Keywords

Genome Medicine

ISSN: 1756-994X

Contact us