NIH Public Access Author Manuscript Infect Genet Evol. Author manuscript; available in PMC 2008 March 1.
NIH-PA Author Manuscript
Published in final edited form as: Infect Genet Evol. 2007 March ; 7(2): 271–278.
Temporal Trends in Gonococcal Population Genetics in a High Prevalence Urban Community Marcos Pérez-Losada, PhD1, Keith A. Crandall, PhD1,2, Jonathan Zenilman, MD3, and Raphael P Viscidi, MD4 1Department of Integrative Biology, Brigham Young University, Provo, UT 2Department of Microbiology and Molecular Biology, Brigham Young University, Provo, UT 3Department of Medicine, Division of Infectious Diseases, Johns Hopkins University School of Medicine, Baltimore MD 4Department of Pediatrics, Johns Hopkins University School of Medicine, Baltimore MD
NIH-PA Author Manuscript
Abstract Molecular evolutionary studies can provide insights into the spread of infectious diseases and inform infection control measures. We performed a population genetic analysis of gonococcal isolates obtained over a 15-year interval in Baltimore MD, where gonorrhea is highly prevalent. Categorical analysis of genetic differentiation revealed temporal structuring of the gonococcal population. The use of a new method to determine the historical demography of N. gonorrhoeae from sequence data showed a strong correlation with trends in the number of reported cases of N. gonorrhoeae and, in addition, may reflect the influence of social and demographic factors and the impact of antimicrobial resistance on the molecular epidemiology of gonorrhea in Baltimore over the past two decades. The strong correlation between the population genetic inferences over the last twenty years and the demographic data collected over the same time period demonstrates the utility of these approaches for the accurate inference of complex population dynamics using multilocus sequence data. The real time application of population genetic analysis can provide sentinel data on gonococcal prevalence, antibiotic resistance patterns and changing epidemiology of gonococcal infections.
Keywords
NIH-PA Author Manuscript
N. gonorrhoeae; gonorrhea; population genetics; MLST; Bayesian skyline plot model and antimicrobial resistance
INTRODUCTION Gonorrhea is a disease with an estimated worldwide annual incidence of 62 million cases (Anonymous, 1995). The Centers for Disease Control and Prevention estimate that approximately 700,000 cases of gonorrhea (twice the reported number) occurred in the United States in 2004 (Division of STD Prevention, 2005). The etiologic agent of gonorrhea is Neisseria gonorrhoeae, a fastidious Gram negative diplococcus, whose only natural host is humans. Spread of the organism occurs directly from person-to-person usually by sexual
Corresponding author: Raphael P. Viscidi, M.D. The Johns Hopkins Hospital, Blalock Building Room 1105, 600 North Wolfe Street, Baltimore MD 21287, 410-614-1494 FAX 410-955-3723,
[email protected]. Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Pérez-Losada et al.
Page 2
NIH-PA Author Manuscript
contact (for review see (Handsfield and Sparling, 1995). In the absence of a vaccine, prevention and control of infection within communities depends on surveillance for prompt identification and rapid treatment of infected individuals in order to interrupt chains of transmissions. A better understanding of the spread of N. gonorrhoeae in human populations can have direct application to public health control efforts.
NIH-PA Author Manuscript
Molecular genetic approaches can provide a powerful tool to assess spread of infectious diseases and identify epidemiological patterns that can inform infection control efforts. If genes of an infecting agent evolve in a predictable manner during transmission from one human to another, then the reconstructed phylogeny will tell how the organism spread. Phylogenetic analysis has been widely and successfully applied to studies of infectious diseases epidemiology (Letiner, 2002), and among other applications has provided insights into past population dynamics (Rambaut et al., 2004) and been used to predict future evolutionary trends (Bush et al., 1999). However, the performance of phylogenetic methods can be adversely affected by recombination, as it violates the underlying assumption of a bifurcating genealogical history (Posada and Crandall, 2002). Since recombination is a dominant evolutionary force that structures genetic variation of N. gonorrhoeae populations (O'Rourke and Stevens, 1993; O'Rourke and Spratt, 1994; Posada et al., 2000; Perez-Losada et al., 2005), phylogenetic analysis is of limited use in molecular epidemiological studies of this organism. Here we apply molecular evolutionary analysis to obtain insights into the spread of the gonococcus through studies of bacterial population structure. Gonorrhea does not spread in a uniform and random fashion in human populations; rather its transmission is strongly influenced by sources of heterogeneity, which may be behavioral, spatial or temporal in nature (Rothenberg, 1983; Rothenberg and Potterat, 1988; Potterat et al., 1985; Aral and Holmes, 1990). Heterogeneity in transmission is expected to give rise to population subdivision and where there is population subdivision there is almost inevitably some genetic differentiation among subpopulations.
NIH-PA Author Manuscript
Baltimore, Maryland has historically had high gonorrhea rates (Becker et al., 1998). For over two decades, we have collected and stored gonococcal strains obtained from clients attending Baltimore City STD Clinics. To study the evolution of N. gonorrhoeae in this high prevalence urban community, we analyzed strains from a time window spanning 15 years. We examined three different classes of molecular markers: porB, a gene under positive selection (Posada et al., 2000; Smith et al., 1995), 7 housekeeping genes which are presumably neutral evolving genetic markers (Viscidi and Demma, 2003; Perez-Losada et al., 2005), and gyrA and parC, the target loci for fluoroquinolone resistance (Belland et al., 1994). Fluoroquinolone resistant N. gonorrhoeae strains were first detected in the last year of our survey period. To test for genetic differentiation, we first estimated population genetic parameters of overall population structure, including allele frequencies, nucleotide diversity, mutation and recombination. Then we tested for population substructure using a categorical analysis of allele frequencies, and also estimated past population dynamics using a recently developed method that estimates effective population size through time from heterochronous sequence data.
MATERIALS AND METHODS Isolate Collection and Genes Gonococcal strains were obtained from two clinics, the Eastern clinic and the Druid clinic, that are run by the Baltimore City Health Department. Approximately 40% of reported gonorrhea cases in the city of Baltimore are diagnosed at the two clinics, with an approximately equal number from each clinic. The total number of cases of N. gonorrhoeae reported by the Baltimore City Health Department in the survey years of 1991, 1996, 2001 and 2005 were 13,007, 6495, 5914 and 3,489, respectively. Thirty-nine isolates of N. gonorrhoeae were obtained from October-November 1991, 76 isolates from June-September 1996, 51 isolates Infect Genet Evol. Author manuscript; available in PMC 2008 March 1.
Pérez-Losada et al.
Page 3
NIH-PA Author Manuscript
from April-November 2001, and 32 isolates from October-November 2005. All isolates were obtained from heterosexual men. Thirty-eight isolates, collected exclusively in 1996, came from clients attending the Eastern Baltimore City STD Clinic and the remaining 160 isolates were from clients attending the Druid STD Clinic (West Baltimore). The study was approved by the Institutional Review Board of the Johns Hopkins University School of Medicine.
NIH-PA Author Manuscript
We sequenced partial regions of 7 of the core housekeeping genes described in Viscidi and Demma (Viscidi and Demma, 2003). The genes were fumC, gdh, glnA, gnd, pilA, pyrD, and serC. These genes were selected based on our previous survey of polymorphism in N. gonorrhoeae housekeeping genes and with the intent to include markers distributed throughout the N. gonorrhoeae FA1090 reference genome. We also sequenced partial regions of gyrA and parC, which encompassed the fluoroquinolone resistance-determining regions, and a partial region of the porB, which spanned hypervariable loops 3 to 6. Genomic DNA was prepared from N. gonorrhoeae cells by using a Clontech nucleospin tissue kit (BD Biosciences, Clontech, Palo Alto, Calif.). PCR reactions were performed in a 50-μl volume containing ∼10 ng DNA, 0.5 μm forward and reverse primer (Table 1), 1X Expand PCR buffer (Boehringer Mannheim Biochemicals, Indianapolis Ind.), 0.2 mM dNTPs, and 1.75 U Expand High Fidelity enzyme (Boehringer Mannheim Biochemicals). Reaction conditions were 94°C for 1 min, and then 94°C for 40 s, 55°C for 40 s, and 72°C for 2 min for 30 cycles and a final 10 min at 72° C. The nucleotide sequences were determined on each DNA strand in a 15-μl volume containing 2 μl DNA, 5 μl BetterBuffer (Gel Company, San Francisco, Calif.), 0.167 μM sequencing primer (Table 1) and 1 μl BigDye Ready Reaction mix (PE Biosystems, Foster City, Calif.). Reaction products were separated and detected with an ABI Prism 3700 automated DNA sequencer from PE Biosystems. Trace data were edited and nucleotide sequences were assembled with the SeqMan software program (DNASTAR Inc., Madison, Wis.). Genetic Analysis A total of 1,856 unique DNA sequences were used for the genetic analyses; 1,162 housekeeping gene sequences and 136 porB sequences were previously reported (Perez-Losada et al., 2005), and 224 housekeeping gene sequences, 306 fluoroquinolone resistance gene sequences and 28 porB sequences were generated for this study. Unique sequences were deposited in GenBank under accession numbers DQ644023 through DQ644038 for housekeeping and fluoroquinolone resistance genes, and DQ645752 through DQ645782 for porB. Sequences were translated into amino acids using the universal reading frame in MacClade 4.05 and then aligned in ClustalX using the default options. Gaps were introduced by Clustal within the porB sequences and the alignments were slightly refined by eye. The 7 housekeeping genes and the 2 fluoroquinolone resistance genes were separately concatenated into two sequences for each isolate.
NIH-PA Author Manuscript
Preliminary analyses showed no significant genetic differences between strains from the two clinics, therefore, all strains were combined for subsequent analyses. Models of nucleotide and codon substitution were assessed using the maximum likelihood approach (Huelsenbeck and Crandall, 1997) with best-fit models selected by the computer program Modeltest v3.06 (Posada and Crandall, 1998) using Akaike information criterion (Posada and Buckley, 2004). The substitution model GTR+⌈+I was selected as the best-fit model of molecular evolution. Nucleotide sequence variation (Na = number of alleles, Ps = number of polymorphic sites and π = nucleotide diversity per site) was estimated separately for the concatenated housekeeping genes, the concatenated fluoroquinolone resistance genes and the porB gene. Genetic diversity (θ), recombination rates (r), and the ratio of the per site chance of recombination to the persite chance of mutation (c/μ) were estimated for each set of genes under the best–fit model of nucleotide substitution using the maximum likelihood coalescent approach implemented in
Infect Genet Evol. Author manuscript; available in PMC 2008 March 1.
Pérez-Losada et al.
Page 4
LAMARC ver2.0.2 (Kuhner et al., 2004). An estimate of the diversity generated per recombination (C) was obtained by multiplying r and θ (for review see (Posada et al., 2002)).
NIH-PA Author Manuscript
To detect genetic subdivision by year of collection we used a χ2 test of sequence absolute frequencies (Hudson et al., 1992) as implemented in Achaz et al. (2004). Simulations by these latter authors showed that this test is robust to variation in samples sizes, asymmetry of sample sizes, and mutation rates, and to presence/absence of recombination, and that the power to detect temporal structure is high. We performed 10,000 permutations using the algorithm of Roff and Bentzen (Roff and Bentzen, 1989) implemented in the program CHIPERM (David Posada, available upon request).
NIH-PA Author Manuscript
The past population dynamics of N. gonorrhoeae in Baltimore were inferred using the Bayesian skyline plot model (Drummond et al., 2005) as implemented in BEASTv1.3 . This coalescentbased demographic model uses standard Markov chain Monte Carlo (MCMC) sampling procedures to estimate a posterior distribution of effective population size through time directly from heterochronous sequence data under a best-fit substitution model. The hyperparameter m was set to ¼ of the sequences in each data set. Four independent MCMC analyses 2 × 107 steps long were performed sampling every 1,000th generation, with the burn-in set at 2 × 106 generations. Log normal relaxed clock evolutionary models were used (Drummond et al., 2006). All the Bayesian MCMC outputs generated by BEAST were analyzed in Tracer v1.3 (Rambaut and Drummond, 2003). BEAST assumes no recombination, an assumption clearly violated by N. gonorrhoeae. Hence before we performed our BEAST analysis we tested for the presence of recombinant sequences and gene regions in all our data sets by using RDP3 (Martin et al., 2005). All the sequences and regions detected by at least two of the six methods included in this program as potential recombinants were excluded in the subsequent BEAST analysis. This criterion is more conservative than the default option suggested by the authors (at least three methods). After running RDP3, 8 to 14% of the sequences were excluded from each data set and ∼0.5 kilo bases were deleted from the MLST concatenated data set. All of the remaining sequences were used in the analysis including identical haplotypes from more than one subject. N. gonorrhoeae Surveillance Data: Cases of N. gonorrhoeae reported by the Baltimore City Health Department for the years 1970 to 2005 were obtained from the Division of STD/HIV Prevention, Centers for Disease Control and Prevention.
RESULTS
NIH-PA Author Manuscript
As we have observed in previous analyses of N. gonorrhoeae sequences (Posada et al., 2000; Perez-Losada et al., 2005), neighbor joining trees of housekeeping genes and porB, constructed using the best-fit substitution model, did not show clear temporal phylogenetic structuring (data not shown). Similarly negative results for temporal structuring were obtained using network approaches, despite their ability to account for moderate levels of recombination (data not shown). To test for temporal differences in population structure, we estimated population genetic parameters, including polymorphism (Na and Ps), nucleotide diversity (θ), genetic diversity (π), and recombination (r and C), for the housekeeping genes, the fluoroquinolone resistance genes and porB (Table 2). As expected for slowly evolving genes, the housekeeping genes showed a low level of polymorphism and the value of the recombination parameter (r) was greater than one indicating that the per-site chance of recombination is greater than the per-site chance of mutation, consistent with the proposed nonclonal population structure of N. gonorrhoeae (O'Rourke and Stevens, 1993). In terms of temporal trends, values for genetic diversity and diversity due to recombination were similar between years. The π and C values for the fluoroquinolone resistance genes were low and did not show differences between years, but the LAMARC estimates for these two genes varied considerably between runs and thus
Infect Genet Evol. Author manuscript; available in PMC 2008 March 1.
Pérez-Losada et al.
Page 5
NIH-PA Author Manuscript
are not very reliable. The porB gene showed a higher level of polymorphism and higher rates of genetic diversity and recombination than the housekeeping genes and the ratio of the persite chance of recombination to the per-site chance of mutation indicates that mutation plays a greater role than recombination in the evolution of porB. Although C values were similar between years, the 1996 sequences showed a greater level of genetic diversity than the 1991, 2001 and 2005 sequences. We next tested for evidence of temporal substructure using a categorical measure of genetic differences among subpopulations. The χ2 test of sequence absolute frequencies was significant (P
