Information on populations and ethics for samples collected previously are included in the publications by Dean et al. (2009) and Batteiger et al. (2014). Russian women were enrolled in the original studies following verbal informed consent after approval by the Local Institutional Review Board at DO Ott Institute of Obstetrics and Gynecology and the Russian Academy of Medical Sciences (RAMS), St. Petersburg, Russia (Shipitsyna et al., 2007). The studies were additionally approved by the Department of Clinical Investigations and Intellectual Property, St. Petersburg Medical Academy of Postgraduate Studies (North-Western State Medical University named after I.I. Mechnikov since 2011), under the Federal Agency of Public Health and Social Development of Roszdrav, in accordance with the Declaration of Helsinki. For the 60 Dutch samples, each was obtained as previously described (Dean et al., 2009) according to the Declaration of Helsinki. The eight Boston samples were sent to Dr. Dean as de-identified samples with no trace to patient name and were not considered human subjects according to IRB at Children’s Hospital Oakland Research Institute and NIH guidelines.

For the Russian samples, endocervical swabs were obtained consecutively from January 2006 to January 2008 in two university clinics in St. Petersburg, Russia, as described elsewhere (Smelov et al., 2009). To minimize the potential for low response from the enrolled women, samples were collected without obtaining personal information. After removal of any mucopus with a cotton swab, a Dacron swab was inserted into the endocervix, rotated and placed in an empty 5 mL vials. Specimens were kept at 4°C (39°F) for up to 4 days before they were shipped to the laboratory where they were stored at 4°C (39°F). Within 1–3 days (Shipitsyna et al., 2007) all samples were tested for the presence of Ct by commercial NAAT assays (Shalepo et al., 2006; Smelov et al., 2009). Additionally, either a conventional PCR (Lytech, Moscow, Russia) or a real-time PCR (Central Research Institute of Epidemiology, Moscow, Russia) were used. The results were confirmed in Amsterdam, Netherlands by the commercial real-time PCR (TaqMan, Applied Biosystems, United States) (Morré et al., 1999; Smelov et al., 2009) and CE-IVD certificated Presto CT-NG Assays (Goffin Molecular Technologies, Beek, Netherlands) (de Waaij et al., 2015).

A total of 323 reference and clinical samples were analyzed that included 58 endocervical samples from St. Petersburg, 60 additional endocervical samples from Amsterdam and eight endocervical samples from Boston, Massachusetts; 265 were already in the MLST database¹ including the 20 reference strains: A/SA-1, A/HAR13, B/TW5/OT, Ba/Apache2, C/TW3/OT, D/UW3/Cx, Da/TW-448, E/Bour, F/IC-Cal3, G/UW57/Cx, H/UW4/Cx, I/UW12/Ur, Ia/UW202, J/UW36/Cx, Ja/UW92, K/UW36/Cx, L₁/440, L₂/434, L_2a/UW396, L₃/404.

Genomic DNA was purified from clinical isolates and ompA genotyped as described previously (Dean et al., 2009; Batteiger et al., 2014). Only the samples from St. Petersburg were provided as swabs; DNA was extracted and purified for these samples using our protocol as described in Joseph et al. (2014). MLST analysis examined seven housekeeping genes: glyA, mdhC, pdhA, yhbG, pykF, lysS, and leuS, with primers as described (Dean et al., 2009) (Supplementary Table 1). All seven housekeeping genes were amplified and sequenced as described (Batteiger et al., 2014). A consensus sequence was created from the forward and reverse sequence. The genes for each of the St. Petersburg, Amsterdam and Boston samples were each concatenated and queried against the 265 samples in the MLST database in addition to including these samples in the database. Allelic numbers and STs were assigned based on this query as described previously (Batteiger et al., 2014).

Using the concatenated sequences, dataset strain clustering and Single Nucleotide Polymorphism (SNP) analyses were performed as described (Batteiger et al., 2014). Briefly, this included visualizing clusters of related STs and non-related STs using eBURST² (Feil et al., 2004). Founder STs were identified by the highest number of single locus variants (SLV) branching from that particular ST (i.e., the clonal ancestor that diversifies into other STs). Clonal complexes generated by eBURST were defined as a group of STs separated by one SLV.

Phylogenetic trees were created by Maximum likelihood using the Symmetric+GI model, which provided the best fit for the data, in the R package phangorn (Schliep, 2010) to analyze the nucleotide sequence variation between the seven MLST loci for each ST. Tree nodes were verified with 1,000 bootstrap replicates. Alternative evolutionary pathways, such as horizontal gene transfer, were analyzed with SplitsTree³ using the splits decomposition method as described (Dean et al., 2009; Tamura et al., 2013). In addition, the sequence for each of the seven MLST loci for a sample were compared across the dataset and to the ompA genotype of the same sample to determine evidence for putative recombination.

Fisher’s exact test was performed in R⁴ to test for significant region-specific ompA and ST clustering; a p-value of <0.05 was considered significant. Confidence intervals were determined based on the method of Clopper-Pearson (Borkowf, 2006). Simpson’s Diversity Index, D, was calculated for the MLST data as described (Simpson, 1949; Hunter and Gaston, 1988). A D-value of ≥0.95 was considered ideal for molecular typing techniques (van Belkum et al., 2007). The ompA genotypes were excluded from the analysis. The Benjamin-Hochberg FDR method (Benjamini and Hochberg, 1995) was used to correct p-values for multiple comparisons.

Samples were classified as putative recombinants when the sequences of the seven gene sequences that comprise the ST or any of the seven genes were non-concordant with each other or with the ompA genotype of the same sample.

PROC FREQ in SAS was used to identify SNPs associated with disease phenotype and Haplotype as described previously (Dean et al., 2009). Levene’s test evaluated the variance across the dataset of the 323 samples. The Classification Index was used to determine significance of each SNP with a disease phenotype where a p-value of <0.05 was considered significant.

DnaSP v5.10 (Librado and Rozas, 2009) was used to calculate nucleotide (nt) and haplotype (hd) to determine the genetic diversity and differentiation for regional subgroups on the concatenated sequences of the seven MLST genes. DnaSP considers the frequency of variants (STs) present in a population and also genetic distances that separate these variants from each other. Genetic population differentiation between regional subgroups was assessed using the pairwise fixation index (Fst) in Arlequin v3.5 (Excoffier and Lischer, 2010) with significance testing by permutation.

The characteristics of the alleles for each gene locus are shown in Supplementary Table 2. The number of alleles varied by gene locus, ranging from seven to 18, as did the number of polymorphic sites. We determined the allele frequencies by geographic region for the 78 alleles (Table 1). Thirty-two (41%) alleles were observed once. The highest number of unique alleles for a geographic region was in Western Europe at 16 alleles but the highest frequency was 43.3% for Russia, which was not statistically significant.

For the 323 samples, 84 unique STs were identified (Supplementary Table 3). STs novel to the dataset were numbered consecutively in order of identification. Of those 84, 57 (67.9%) were singletons, with a relatively even distribution by geographic region (excluding Asia and Africa where the sample sizes were very small) with a higher percentage of singletons in Europe that was not significantly different (Table 2; p = 0.08). Table 2 also shows that the percentage of unique STs per region was highly similar. However, the ST/sample ratio was significantly greater for Russian than for European and American samples (p < 0.009). There were also significant differences in the distribution of STs (Table 3). Dutch females (n = 79) were significantly more likely to be infected with ST23 (p = 0.032; 95% CI: 1–23) and ST34 (p = 0.02; 95% CI: 1.1–25) compared to Russian females (n = 58) and with ST19 (p = 0.001; 95% CI: 1.7–34.7) compared to American females (n = 108). Supplementary Table 4 shows the distribution of STs by geographic region.

There were 26 ompA genotypes observed in the dataset resulting in a 3.23 lower resolution than for STs. Excluding the STDs samples from South Africa, all samples available from Asia and Africa were from trachoma patients; 1 (7.1%) of 14 was a urogenital Da ompA genotype (ST 37) in Asia while 1 (5.9%) of 17 was a urogenital E genotype (ST 39) in Africa.

ompA genotype E (n = 72) was the most prevalent and associated with 18 STs with an ST to sample ratio of 0.25 (Supplementary Table 3). The distribution for the remaining ompA genotypes in descending order was G (38; 16 STs; ratio: 0.42); D (33; 16 STs; ratio: 0.48); F (30; 5 STs; ratio: 0.17); Ia (17; 4 STs; ratio: 0.25); K (15; 7 STs; ratio: 0.47); J (13; 9 STs; ratio: 0.69); B, H and I at 12 samples each (7, 6, and 3 STs, respectively; ratios: 0.58, 0.50, and.25, respectively); C (11: 5 STs; ratio: 0.45); A (9; 4 STs; ratio: 0.44); D2 and Ja at 4 samples each (1 and 3 STs, respectively; ratios: 0.25 and 0.75, respectively); Ba, D1, E6, F4 and Ia4 at 2 samples each (1 ST for each; ratio: 0.50 for each); and Da (1; 1 ST; ratio: 0.10).

The ompA distribution of urogenital strains varied across Europe, the Americas, and Russia. In Europe, ompA genotype D was significantly more prevalent than in other geographic regions (p = 0.046) and more frequent than the globally prevalent genotype E. In Russia, ompA genotype G was significantly more prevalent than in all other regions (p = 0.001). In the Americas, ompA genotype Ia was significantly more prevalent than the other regions (p = 0.001). In comparing female STD cohorts, Russian women were significantly more likely to be infected with E and G (p = 0.001 and 0.026, respectively) while Dutch women were significantly more likely to be infected with D and I (p = 0.026 and 0.002, respectively).

Simpson’s Discriminatory Index was 0.98 (95% CI 0.97–0.99) for MLST and 0.70 (95% CI 0.67–0.73) for ompA genotyping.

Based on DnaSP, nucleotide diversity was lower in Russia (Pi = 0.00216) compared to Netherlands (Pi = 0.00321) and North America (Pi = 0.00294) (Table 4). However, Netherlands and North American datasets contained isolates from men and LGV samples. The phylogenetically disparate LGV samples appeared to contribute substantially to nucleotide diversity; nucleotide diversity dropped from 0.00250 to 0.00196 in Dutch women when the four LGV isolates were removed. When comparing non-LGV isolates from women in these regions, Russian women had the highest nucleotide diversity (Pi = 0.00216), followed by North American women (Pi = 0.00210), and Dutch women (Pi = 0.00196).

Assessing population differentiation between regional subgroups by Fst revealed significant differences between most regional subgroup pairwise comparisons (Table 5). African and Asian subgroups exhibited high Fst values in all comparisons, indicating those regions were distinct from others in the dataset. This is not surprising given the small sample size and the fact that these samples are from two distinct trachoma populations. Though Western Europe was significantly different from all subgroups, the Russia and Americas subgroups were not significantly different by Fst. However, this could be attributed to the higher LGV representation in the Western European samples (15.5%) in comparison to the American samples (7.3%), and Russian samples (0%). Concordantly, when LGV samples were removed from the Western European dataset, Western Europe was no longer significantly different from Russia (Table 5). When evaluating only women without LGV from the Netherlands, Russia, and North America, the Russian women were not significantly different from the American and Netherlands women (Supplementary Table 8).

The phylogenetic relationships were initially evaluated by eBURST, which revealed clonal clusters (CC) similar to what we reported previously (Figure 1) (Dean et al., 2009; Batteiger et al., 2014) but with the addition of an LGV cluster. These included CC-A encompassing trachoma STs, CC-B with non-invasive, non-prevalent urogenital STs, CC-C with non-invasive prevalent STs and CC-D that included LGV STs. The predicted founders were ST19, ST23, ST34, and ST39 (Supplementary Figure 1 and Table 5). The 57 singleton STs are denoted by a gray circle alone or within a colored circle, representing a specific geographic region, but were not associated with any specific region.

The tree (Figure 2) revealed ST branches similar to the eBURST clusters; both had branches or clusters for the disease phenotypic groups of LGV, non-invasive prevalent urogenital, and non-invasive non-prevalent urogenital STs. However, the trachoma STs formed a subgroup of the non-invasive non-prevalent urogenital branch. In addition, within disease groups, STs branched from central nodes by geographic region. These nodes contained, in general, large numbers of STs from diverse locations. For example, the founders ST19, ST23, ST34, and ST39 located at nodes in the tree contain STs from Europe, Russia, and the Americas. The amino acid tree showed a similar phylogeny (Supplementary Figure 1).

The Splitstree decomposition tree revealed evidence for a network structure consistent with homologous recombination (Figure 3). This was demonstrated by the interconnecting networks specifically among the ST founders ST19, ST23, ST34, and ST39 and other STs on the network in addition to the canonical evolutionary pathway shown in the tree (Figure 2). Supplementary Figure 2 shows the amino acid tree.

The Splitstree data are consistent with findings in the MLST and ompA genotyping data. Samples were classified as recombinant when the sequences of the seven genes that determine the ST were non-concordant with each other or with the seven gene sequences of the Ct genotype associated with the ompA genotype (Table 6), denoted in bold in Supplementary Table 3. There was no evidence of any recombination among the seven ST genes for any sample, although this would conceivably be possible.

A total of 106 (32.8%) samples were considered putative recombinants. Excluding Asia and Africa where the sample sizes were small, Europe had a significantly higher number of recombinants than Russia and the Americas (p < 0.04) (Table 2), where, of the 109 samples, 49 were recombinant with 43 (87.8%) from the Netherlands. While ompA genotypes were generally consistent across samples of the same ST, many cases of recombinant strains were observed. For example, ST19 (n = 40) was primarily associated with ompA genotype G (37.5%) but eight different ompA genotypes (B, D, E, G, H, I, J, K) were also associated with this ST. The majority of these samples were from St. Petersburg and Amsterdam (Supplementary Table 3). ST23 (n = 32) was also associated with eight ompA genotypes (B, D, G, H, I, Ia, J, K) but the most frequent was Ia (47%). In contrast, the most geographically prevalent ST was ST39 (n = 45) where 95% were associated with ompA genotypes E. There were also 24 singletons that were recombinants (Table 6).

Table 7 shows the ompA genotype, ST and allelic SNPs, if present, associated with each of the eight Boston samples added to the dataset. Under the column denoted ST sequence homology are the Ct ompA genotype sequences to which the sequences of the seven MLST genes are identical. For example, sample J/259b has four SNPs that are different from the sequences of the seven genes for reference and clinical J strains in the database. These SNPs are identical to the sequences of strains G and K, which suggests that this sample is a recombinant between J and G or K strains. Another example is sample D/256b; the ompA genotype is D but the MLST sequences of the seven genes are identical to reference and clinical F strains in the database. Therefore, it is a recombinant between D and F. Five (62.5%) of the eight Boston samples were putative recombinants. Supplementary Table 5 shows similar data for the 60 new samples from Amsterdam. For the Amsterdam samples, 43 (54.4%) of 79 samples were recombinants. Likewise, Supplementary Table 6 represents the data from St. Petersburg where 20 (34.5%) of 58 samples were recombinants.

Table 8 shows the SNPs that are specific for each Haplotype and disease phenotype group. These groups include invasive urogenital disease caused by LGV strains (Haplotype 1), non-invasive urogenital disease group that includes the most prevalent worldwide strains including D, Da, E, and F (Haplotype 2), and trachoma caused by ocular A, B, Ba, and C strains (Haplotype 3). The non-invasive group includes clinical Ja strains that are recombinants of E and F strains. There are 14 SNPs that are specific for the LGV disease group, and each SNP independently identifies this group and Haplotype. The same is true for the non-invasive disease group. However, the two SNPs for the trachoma group together are specific for this group and Haplotype 3. Based on the Classification Index, the SNPs were not uniformly distributed.