WGS to predict antibiotic MICs for Neisseria gonorrhoeae

Background Tracking the spread of antimicrobial-resistant Neisseria gonorrhoeae is a major priority for national surveillance programmes. Objectives We investigate whether WGS and simultaneous analysis of multiple resistance determinants can be used to predict antimicrobial susceptibilities to the level of MICs in N. gonorrhoeae. Methods WGS was used to identify previously reported potential resistance determinants in 681 N. gonorrhoeae isolates, from England, the USA and Canada, with phenotypes for cefixime, penicillin, azithromycin, ciprofloxacin and tetracycline determined as part of national surveillance programmes. Multivariate linear regression models were used to identify genetic predictors of MIC. Model performance was assessed using leave-one-out cross-validation. Results Overall 1785/3380 (53%) MIC values were predicted to the nearest doubling dilution and 3147 (93%) within ±1 doubling dilution and 3314 (98%) within ±2 doubling dilutions. MIC prediction performance was similar across the five antimicrobials tested. Prediction models included the majority of previously reported resistance determinants. Applying EUCAST breakpoints to MIC predictions, the overall very major error (VME; phenotypically resistant, WGS-prediction susceptible) rate was 21/1577 (1.3%, 95% CI 0.8%–2.0%) and the major error (ME; phenotypically susceptible, WGS-prediction resistant) rate was 20/1186 (1.7%, 1.0%–2.6%). VME rates met regulatory thresholds for all antimicrobials except cefixime and ME rates for all antimicrobials except tetracycline. Country of testing was a strongly significant predictor of MIC for all five antimicrobials. Conclusions We demonstrate a WGS-based MIC prediction approach that allows reliable MIC prediction for five gonorrhoea antimicrobials. Our approach should allow reasonably precise prediction of MICs for a range of bacterial species.


Introduction
Antimicrobial-resistant Neisseria gonorrhoeae is a risk to public health, particularly as emerging resistance to available treatment such as azithromycin and ceftriaxone leaves few treatment options. 1 Treatment and control strategies depend on reliable monitoring of trends in N. gonorrhoeae antimicrobial resistance, which are a major focus of national surveillance programmes, such as GRASP (England and Wales), 2 Euro-GASP (Europe), 3 GISP (USA) 4 and the Canadian national programme. 5 Several studies have demonstrated the potential of WGS to predict antimicrobial susceptibilities across a range of pathogens, including N. gonorrhoeae, 6,7 Staphylococcus aureus, 8 Enterobacteriaceae 9 and Mycobacterium tuberculosis. 10 N. gonorrhoeae antimicrobial resistance mechanisms are well described, 11 allowing their identification using WGS from national surveillance collections, 7 cohorts with cefixime 12 and azithromycin resistance, 6 and reference collections. 13 However, the reported ability of WGS to predict antimicrobial susceptibilities in N. gonorrhoeae has been variable, e.g. single-resistance determinants accurately predict cefixime, but not azithromycin susceptibilities. 7 Most WGS studies predict antimicrobial susceptibility categorically as 'susceptible' or 'resistant'. While useful for managing individual patients, this provides less surveillance information than phenotypic MIC measurements, which allow trends over time, and the development and spread of resistance to be monitored with more precision. This is particularly important where prevalent strains have antibiotic MICs near clinical breakpoints.
WGS has been successfully used to predict b-lactam MICs in Streptococcus pneumoniae 14,15 and azithromycin MICs in N. gonorrhoeae. 6 Therefore, using sequence collections from England, the USA and Canada, we investigated whether WGS and simultaneous analysis of multiple resistance determinants using multivariate regression can be used to predict antimicrobial susceptibilities to the level of MIC for multiple drugs in N. gonorrhoeae.

Sample collections and antimicrobial susceptibility testing
Three collections, from Brighton, England, 16 the USA 12 and Canada 6 with antimicrobial susceptibilities determined as part of national surveillance programmes, and WGS from previous studies, were studied. The Brighton samples were unselected consecutive clinical isolates from 3 months of each calendar year, the US samples were enriched for cefixime resistance and the Canadian samples for azithromycin resistance. MICs were determined by agar dilution, for Brighton using the GRASP method, 17 and in the USA 18 and Canada using CLSI methods. 6 Quality control was performed using the 2008 WHO gonorrhoea reference strain panel 19 and ATCC 49226, with MICs and sequences obtained for these strains in Canada used in the final analysis. The distribution of MICs of cefixime, penicillin, azithromycin, ciprofloxacin and tetracycline was sufficient to determine predictors of resistance. There were insufficient numbers of samples with ceftriaxone resistance, and therefore insufficient numbers with important resistance determinants, to include it in the analysis. EUCAST breakpoints were used to categorize samples as susceptible and resistant: azithromycin, 0.25 and .0.5; cefixime, 0.125 and .0.125; ciprofloxacin, 0.03 and .0.06; penicillin, 0.06 and .1; and tetracycline, 0.5 and .1 mg/L. 20

Genetic determinants of antimicrobial susceptibility
Previously reported genetic determinants conferring reduced antimicrobial susceptibility to the different agents are shown in Table 1. Using WGS, SNPs and resultant amino acid substitutions were determined following mapping to the NCCP11945 reference genome (NC_011035.1) and quality filtering as in De Silva et al., 16 with the exception of the penA and penB genes and mtrR promoter variants, which were more variable, and were therefore identified from Velvet 21 de novo assemblies using BLAST and subsequently aligned using MUSCLE. 22 BLAST searches of de novo assemblies were used to determine the presence/absence of accessory genes. To identify rRNA variants, sequence reads were mapped against a single copy of the NCCP11945 23S RNA gene using BWA mem 23 with default settings. Base counts were determined using SAMtools, 24 enabling estimation of the proportion of gene copies with relevant mutations.

Statistical methods
Multivariate linear regression models were used to identify genetic predictors of MIC. Measured MIC values were converted to a log 2 scale, such that a one-unit change is equivalent to a single doubling dilution, e.g. log 2 (MIC) " 3 represents MIC " 8 mg/L and log 2 (MIC) " 4 represents MIC " 16 mg/L. Log 2 (MIC) values were used as the (approximately normally distributed) outcome in multivariate linear regression models, and the Table 1 genetic determinants as potential predictors, such that the model coefficients associated with each genetic determinant can be used for MIC prediction. The models are additive on the log 2 scale, i.e. the log 2 (MIC) value is predicted by adding the coefficients for each genetic determinant to the model constant term [which is equivalent to the WT log 2 (MIC)]. The predicted MIC is then 2 to the power log 2 (MIC). Predicted MICs are presented rounded to the nearest doubling dilution, rounding log 2 (MIC), prior to conversion back to the absolute scale.
Where phenotypic MIC values were observed to be below or above the quantification limits, e.g. ,0.06 or .16 mg/L, the actual MIC was assumed to be the adjacent doubling dilution, e.g. 0.03 and 32 mg/L, respectively, for the purposes of model fitting. For genetic determinants with multiple non-WT alleles, where there were ,10 observations for a given allele, the variable was collapsed to a binary WT/non-WT variable, identifying amino acid substitutions by the mutation site only. Otherwise, determinants were included as categorical variables with each non-WT allele specified separately. The only exception made to this was for a single isolate in the dataset with penA A501P, 25 as this, combined with penA XXXIV, confers high-level cefixime resistance. A separate coefficient for penA A501P is presented in the results, but not used for model predictions, as only a single sample was available.
penA alleles were determined using a maximum likelihood phylogenetic tree of previously published penA alleles and those available in GenBank. The closest matching published penA allele was identified for each sample using BLAST. Where ,10 samples matched a given allele, this allele was merged with the nearest neighbour on the tree, and the process repeated until all alleles considered had 10 samples ( Figure S1, available as Supplementary data at JAC Online). Previously reported SNPs in penA were also considered as potential additional predictors, as individual SNPs may affect resistance, against a background of otherwise identical/similar penA alleles. 25 Univariate regression coefficients were calculated for all antimicrobialrelevant/specific genetic determinants listed in Table 1. Given the large number of potential genetic determinants, to reduce model over-fitting, the final model was chosen using backwards selection, starting from the model including all genetic determinants, and then removing one determinant at a time, until the lowest possible Akaike information criterion (AIC) was obtained. AIC was used rather than more stringent P value thresholds (such as P . 0.05) because a priori determinants were considered plausibly resistance associated. Alleles were considered as categorical variables, i.e. all alleles were included if the factor as a whole improved the AIC. The rRNA A2059G mutation was only present at 0 or 4 copies and was therefore treated as a categorical variable. The rRNA C2611T mutation was present at 0, 1, 2, 3 and 4 copies, and was modelled as a continuous variable, after ensuring model fit using the AIC was not improved by treating it as a categorical variable. Excluded determinants were added back one at a time to the final model to check this did not improve the AIC; any that did were retained in the final model. Interactions between genetic determinants were then tested, e.g. to detect a saturation effect, where further determinants in the same pathway do not increase the MIC beyond a threshold, or conversely for synergy between determinants. Pairwise interactions between all model terms were tested, initially one at a time; all interactions with a Wald P value ,0.01 were added to the model, and backwards selection repeated, requiring multivariate Wald P values ,0.01 to include an interaction in the final model. This pre-specified approach to interactions is deliberately more conservative than the AIC, which is used to identify main effects, to account for multiple potential interaction pairs tested and minimize overfitting. Errors and model MIC predictions were estimated using leave-oneout cross-validation, where the final model (including interaction terms) is fitted using all samples except one, which is then used to predict the outcome for the excluded sample, repeated over all samples in the dataset. Analyses were performed using STATA version 14.1 (StataCorp, College Station, TX, USA). Eyre et al.  Using WGS, all genetic determinants present in the 2008 WHO reference strain panel 13,19 were recovered with 100% concordance. The univariate coefficients for each genetic determinant and the final multivariate model for each antibiotic are given in Table S1. The numbers of samples included in the final multivariate model with completely determined genotypes were 670 for cefixime, 672 for penicillin (2 ambiguous penA alleles; relevant penA SNPs not determined for 9 and 7 samples, respectively, due to incomplete de novo assemblies), 676 for ciprofloxacin (5 samples without parC variants determined) and 681 for azithromycin and tetracycline.
For cefixime, penA allele, and mtrR45, penA501/512, mtr120 and penB mutations were independent predictors of MIC, with penA XXXIV associated with the greatest MIC increase, 8.3-fold (multivariate regression coefficient 3.06, 95% CI 2.19-3.98, P , 0.001; MIC fold increase calculated by 2 to the power of the regression coefficient, i.e. 2 3.06 ). The single sample with A501P, showed a predicted MIC increased by a further 26.7-fold. penA XV had a negative coefficient (-0.44), indicating this lineage may represent WT MIC, rather than the closely related M32091 ( Figure S1) chosen as the baseline allele. penB121 mutations also had a negative coefficient, but as this finding was based on three samples it may represent model over-fitting. All potential determinants of penicillin MIC were included in the final model with the exception of penB121; the bla TEM b-lactamase increased MICs the most, 23.4-fold (coefficient 4.55, 3.80-5.31, P , 0.0001). rRNA mutations were the strongest determinants of azithromycin MIC, along with erm genes (although these were only found in three isolates), and mtrR mutations and promoter disruption. The rRNA A2059G allele was present in either zero or four copies, four copies increased MICs 375-fold (coefficient 8.55, 7.69-9.40, P , 0.001), and C2611T rRNA mutations present in between zero and four copies, increased MICs 2.9-fold (coefficient 1.52, 1.25-1.79, P , 0.001) per copy. gyrA and parC mutations predicted ciprofloxacin MIC, with gyrA95 mutations having the greatest effect, 16.8-fold MIC increase (coefficient 4.07, 2.86-5.29, P , 0.001). Finally, tetracycline MICs were most impacted by the presence of tetM, resulting in a 126-fold MIC increase (coefficient 6.98, 6.36-7.60, P , 0.001), but also by penB and rpsJ mutations and mtrR mutations and promoter disruption.
In addition, the study country [a composite of the testing laboratory and method: CLSI versus GRASP, the isolate origin, and the sampling frame (US samples enriched for cefixime resistance, and Canadian samples for azithromycin resistance)] was independently predictive of MIC for all agents, and not just those enriched for in the sampling frame. Compared with England, cefixime, penicillin and azithromycin MICs were all higher for the USA and Canada (all P , 0.001), ciprofloxacin MICs were lower for the USA (P , 0.001) and tetracycline MICs were lower for the USA and Canada (P , 0.001). For example, cefixime MICs were 2.6-fold (coefficient 1.36, 1.18-1.54, P , 0.001) and 1.5-fold (coefficient 0.60, 0.42-0.78, P , 0.001) higher for the USA and Canada, respectively, compared with England.
Several interaction terms had negative coefficients (Table S1), i.e. the combination of two genetic determinants produced a smaller increase in MIC than the multiple (sum on the log scale) of their effects when found alone. For example, for tetracycline with rpsJ 57 and tetM, the addition of rpsJ 57 to tetM does not further increase MIC [6.3-fold increase in MIC with rpsJ 57 alone (coefficient 2.66), 126-fold with tetM alone (coefficient 6.98) and 124-fold with both (coefficients 2.66!6.98 -2.69)]. Other examples include tetM with mtrR 45 and mtrR promoter deletion, for ciprofloxacin parC87 and 91 mutations. Positive interaction coefficients indicate potential synergy, e.g. for tetracycline, the combination of penB120 and mtrR promoter deletion increases MIC an additional 1.5-fold (coefficient 0.55) compared with the multiple of their effects alone.
If epistasis, beyond simple pairwise interactions, involving multiple known determinants were to explain the majority of the variation not captured by the models, fixed combinations of known variants would be expected to produce fixed MICs. Therefore, for each antimicrobial, and for each combination of genetic determinants with 10 samples, we plotted observed MICs ( Figures S2-S6), demonstrating similar variability within these fixed combinations to that seen in the model residual plots (Figure 2).

Eyre et al.
We also investigated if analysis of single determinants could better explain which isolates were resistant for cefixime. In keeping with previous reports 12 and other data, 7 penA XXXIV was a sensitive predictor of resistance, present in 106/110 resistant strains [sensitivity 96%, VME rate 3.6% (95% CI 1.0%-9.0%)]. However, taken alone it lacks specificity when all included cohorts, including those not selected on the basis of their cefixime susceptibility, are considered, with 83/560 susceptible strains (42 Brighton, 34 Canada, 7 USA) also containing penA XXXIV [specificity 85%, ME rate 14.8% (12.0%-18.0%)]. Additionally, other resistance Match I Match R Minor ME VME Figure 2. Accuracy of predicted MICs. The difference between the measured and predicted MICs is shown on a log 2 scale, such that a difference of 1 represents a difference of 1 doubling dilution. A VME is where the predicted susceptibility is susceptible when the measured susceptibility is resistant; an ME is where the predicted susceptibility is resistant, but the measured susceptibility is susceptible; and minor errors denote where the phenotype is susceptible and the prediction intermediate, or the phenotype is intermediate and the prediction is resistant, or vice versa in both cases. S, susceptible; I, intermediate; R, resistant; ME, major error; VME, very major error.
Eyre et al.  20 A VME is where the predicted susceptibility is susceptible when the measured susceptibility is resistant; an ME is where the predicted susceptibility is resistant, but the measured susceptibility is susceptible; and minor errors denote where the phenotype is susceptible and the prediction intermediate, or the phenotype is intermediate and the prediction is resistant, or vice versa in both cases. The bottom panel presents the same analysis as the middle panel, but allowing for the true phenotype to be + 1 doubling dilution from the observed phenotype. penA genotypes could not be determined for two samples; these samples were excluded from the cefixime and penicillin models. A further nine and seven samples were excluded from the cefixime and penicillin models, respectively, as the penA SNPs could not be determined due to incomplete de novo assemblies. parC variants could not be determined for five samples; these samples were excluded from the ciprofloxacin model.

Discussion
We demonstrate a WGS-based MIC prediction approach and show that it allows reliable MIC prediction for five gonorrhoea antimicrobials. We predicted MIC within one doubling dilution for 93% of samples, and within two for 98%. This is comparable to routine phenotypic performance: acceptable variation in MIC measurement for type strains in national surveillance programmes is typically +1 doubling dilution and sometimes greater. 27 Combining MIC predictions and clinical breakpoints resulted in VME and ME rates within regulatory targets, 26 with the exception of VME rates for cefixime and ME rates for tetracycline. The majority of cefixime errors arose from single doubling dilution discrepancies in MIC prediction near the breakpoint. Using an alternative approach, 7,12 based on the presence/absence of penA XXXIV as a determinant of cefixime resistance, VMEs approached acceptable levels [3.6% (95% CI 1.0-9.0%)], but at the expense of increased MEs (14.8%). Almost identical VME/ME rates (3.6%/15.2%) were obtained using a breakpoint one doubling dilution lower than the EUCAST value with the predicted MICs. Some of the very marked phenotype/prediction differences may also be the result of laboratory labelling errors, but we were unable to re-phenotype and resequence all such isolates in this study.
We detected MIC variation based on the country of testing, despite phenotypes being determined using gold-standard methods in national reference laboratories. This variation might be partly explained by the sampling frames used, with US isolates enriched for cefixime resistance and Canadian isolates enriched for azithromycin resistance. These isolates might be more likely to contain otherwise unexplained resistance to cefixime and azithromycin, respectively, which is attributed to the country of testing. However, we observed variation across all antimicrobials investigated, and with by-country in MICs following different directions by antimicrobial MICs: higher MICs for cefixime, penicillin and azithromycin in North America, and lower MICs for ciprofloxacin and tetracycline. The differences could also be due to differences between the CLSI testing methods used in North America (coefficients for the USA and Canada tended to be similar), and the GRASP methods used in England and Wales. We chose to model the study country rather than the testing method to be transparent with respect to any possible biases introduced by the country-specific sampling frames. However, there were also differences between the USA and Canada despite the same testing method being used, which may be due to differences in unmeasured strain background in these countries, but may also represent inter-laboratory measurement variation. Determining the exact contribution of each mechanism to differences between countries is not possible with the current study; however, the differences highlight a potential benefit from greater standardization of phenotyping to facilitate evaluation of trends globally.
The genetic determinants with the greatest impact on MIC are broadly in keeping with previously published observations. 6,11 It should be noted that the approach taken does not necessarily define the mechanistically most important resistant determinants, but simply those that are statistically most informative to predict MIC. Where multiple mechanisms typically coexist, this collinearity may lead to only one mechanism remaining in the final model, and its associated coefficient represents the overall effect of these multiple mechanisms.
Our data demonstrate that several different mechanisms can together contribute to determine MICs: a strength of the multivariate regression approach used is that these mechanisms can be accounted for simultaneously. We noted examples where the presence of the same mechanism can result in a MIC just above or below the breakpoint, either due to the presence/absence of other mechanisms as observed for azithromycin in this study, or potentially due to testing methods, e.g. for cefixime. While breakpoints are selected across all strains, based on the likelihood of clinical treatment success, it remains to be seen if some mechanisms are more important than others in determining efficacy of treatment. Wider deployment of WGS may facilitate determination of the extent to which mechanism or absolute MIC is more important in determining treatment success.
The pseudo-R 2 values of 0.80-0.96 obtained suggest that our models explain most of the MIC variation. The residual, unexplained, variation probably relates to a combination of inherent phenotypic variability (which would not be addressed by a better model) and unidentified resistance determinants (or their interactions). Repeat testing of the same collection/subcollection of isolates would allow the former to be measured, and should ideally form part of future studies, to allow the latter to be estimated.
A limitation of the current modelling framework and sampling frames is that the available power allowed only determination of the effect of relatively common resistance determinants. For example, we did not have sufficient numbers of resistant isolates to formulate a similar model for ceftriaxone resistance. To avoid model over-fitting we selected the predictors in our final models based on the AIC; as such, it is possible that exclusion from the model of less common variants may be due to lack of power rather than lack of a true effect. An alternative approach would have been to retain all potential predictors identified a priori in the final model, albeit with an increased risk of model over-fitting.
Another related limitation is that the models are only as good as the resistant determinant catalogue on which they are based. Uncommon, but known, variants may have been excluded, as were likely multiple unknown resistance determinants. These unknown determinants have been referred to as 'factor X' in transformation experiments, where penA, mtrR and penB determinants increase resistance, but do confer donor levels of resistance to penicillin or extended-spectrum cephalosporins. 28 Future WGS has the potential to identify novel resistance determinants through genome-wide association studies. 29 The dataset used also restricts the resistance determinants for which coefficients can be determined. Several known resistance determinants were not Eyre et al. present in our WGS dataset, and we were unable to determine their impact on MIC. One approach that might bridge this shortcoming would be to fit the model in a Bayesian framework where prior information from strains not previously subject to WGS, but with a known mechanism and MIC increment, e.g. from laboratory experiments with isogenic mutants, informs the final model as well as paired phenotypic and WGS data. There is also a need to derive models on more diverse datasets, e.g. many of the cefixime-resistant strains in this study are from a single clone dispersed throughout North America and England, 12,16 and other mechanisms of cefixime resistance may be common elsewhere. Similarly, our approach of assigning each sequence to the most closely related penA allele grouping for cefixime and penicillin MIC prediction is only likely to be appropriate where similar alleles have been previously identified, and therefore sufficiently diverse derivation datasets are required to ensure the probability of encountering a very divergent novel allele is low.
The current study does not include an independent validation dataset, as we chose to include all available samples to maximize power. However, a further independent validation should be undertaken before applying our method to patient care.
WGS may become part of routine diagnostic microbiology workflows in the next 5-10 years. 30 The approach taken in this study could be adopted, with sequencing following culture simultaneously providing information about antimicrobial susceptibilities and strain relatedness. Alternatively, if reliable inferences can be made from sequencing clinical samples directly this might offer further advantages, such as allowing culture-negative samples to be analysed and reducing turnaround times. However, further work will be required to develop methods to undertake WGSbased MIC prediction from clinical samples, including accounting for coexisting human and potentially similar commensal Neisseria species DNA.
In summary, WGS allows MIC prediction for a range of antimicrobials for N. gonorrhoeae. The approach taken here should allow reasonably precise prediction of MICs from genetic data for a range of bacterial species, at least to levels of variation considered acceptable for routine diagnostics. If WGS becomes a widely used diagnostic tool, large amounts of surveillance data may become available from routine clinical activity. WGS may therefore provide reproducible and readily exchangeable data on the spread of antimicrobial resistance, alongside providing informative data on the strains that are carrying it and how resistance is being transmitted.

Data deposition
NCBI short-read archive accession numbers for sequences used in this publication and associated phenotypes and MIC predictions can be found in Table S2.