Distribution of β-Lactamase Genes in Clinical Isolates from California Central Valley Hospital Deviates from the United States Nationwide Trends

The evolution and dissemination of antibiotic resistance genes throughout the world are clearly affected by the selection and migration of resistant bacteria. However, the relative contributions of selection and migration at a local scale have not been fully explored. We sought to identify which of these factors has the strongest effect through comparisons of antibiotic resistance gene abundance between a distinct location and its surroundings over an extended period of six years. In this work, we used two repositories of extended spectrum β-lactamase (ESBL)-producing isolates collected since 2013 from patients at Dignity Health Mercy Medical Center (DHMMC) in Merced, California, USA, and a nationwide database compiled from clinical isolate genomes reported by the National Center for Biotechnology Information (NCBI) since 2013. We analyzed the stability of average resistance gene frequencies over the years since collection of these clinical isolates began for each repository. We then compared the frequencies of resistance genes in the DHMMC collection with the averages of the nationwide frequencies. We found DHMMC gene frequencies are stable over time and differ significantly from nationwide frequencies throughout the period of time we examined. Our results suggest that local selective pressures are a more important influence on the population structure of resistance genes in bacterial populations than migration. This, in turn, indicates the potential for antibiotic resistance to be controlled at a regional level, making it easier to limit the spread through local stewardship.


Introduction
β-lactam antibiotics have been in use since the discovery of penicillin in the 1940s, and they continue to be the most widely used antibiotics due to their high effectiveness, ease of delivery, and low toxicity [1,2]. The longstanding use of β-lactam antibiotics has led to the emergence of resistant strains in clinical care settings [3]. The continuous selection and evolution of β-lactamase genes by β-lactam antibiotic use has led to the diversification of successful β-lactamase genes: bla TEM , bla SHV , bla CTX-M , and bla OXA [4]. β-lactamase genes produce extended spectrum β-lactamase (ESBL) enzymes that work by hydrolyzing βlactam antibiotics, rendering them ineffective. bla TEM and bla SHV were the first β-lactamase enzymes identified in 1963 and 1972, respectively, and were implicit in outbreaks in the 1990s [5][6][7]. Today, bla SHV composes 10% of ESBLs, and bla TEM has become somewhat less common in the U.S. [8]. bla CTX-M was first identified in 1989, and was identified with increasing frequency throughout the 1990s [9]. By the 2000s, the frequency of bla CTX-M enzymes surpassed those of bla TEM and bla SHV . Although first discovered in 1976, bla OXA enzymes have been increasing in prevalence due to the frequent association of bla OXA-1 with bla CTX-M-15 [10,11]. Today, bla CTX-M enzymes are the most identified ESBLs, and have displaced bla TEM and bla SHV in many individual hospitals [6,9,[12][13][14]. However, this trend is not uniform across publications originating from different surveillance locations [11,15,16]. Bajpai et al. (2017) found bla TEM to be the most abundant ESBL enzyme in a single hospital, although other reports detail different ESBL gene frequencies [17]. In the United States, few recent nationwide surveillance studies have specifically examined the frequencies of specific ESBL genes. One recent survey of 26 hospitals identified bla TEM as the most abundant ESBL enzyme in clinical isolates (47%), followed by bla CTX-M (36%), bla SHV (35%), and bla OXA (20%) [18].
Regional variance in the frequencies of ESBLs enables the assessment of which factors are contributing the most to ESBL frequencies. Due to the strong selection that bacteria experience from antibiotics and the rapid migration of bacteria that occurs in human populations, selection and migration were the two factors we chose to investigate. To understand the relative contributions of selection and migration, it was important to obtain and compare updated ESBL gene frequencies. We chose to compare the frequencies of ESBLs in a local repository of ESBL positive isolates collected from a single hospital, with average frequencies nationwide across the U.S. obtained from ESBL-positive clinical isolates whose genomic sequences have been deposited in the NCBI Isolates Browser Database.
When comparing genetic variations over two populations, there are four possible outcomes depending on the stability of gene frequencies within a site and comparisons of those frequencies between sites. First, gene frequencies that are stable over time within a population and non-uniform across populations indicate low migration between bacterial populations and that selection for resistance within a given population is strong and constant. However, if the gene frequencies are unstable over time within a population and non-uniform across populations, this would indicate alternating local selective pressures and rapid migration as "immigrants would increase the mutation supply rate" [19] and would compete with "better-adapted residents maintaining the population away from the local fitness optimum" [20]. Stability over time and uniformity between populations suggest rapid continuous migration between populations and strong consistent selection resulting in a highly resistant and optimized strain [20]. Finally, unstable (alternating) frequencies over time and uniformity between populations indicate strong alternating selective pressures in large areas (or populations). Moreover, this signal also indicates rapid migration because variation between populations averages out as immigration leads to a decrease in genetic differentiation between populations [21]. We compared ESBL gene frequencies from Dignity Health Mercy Medical Center (DHMMC) and the rest of the U.S. over a period of six years as follows.

Regional Gene Frequencies
We performed a molecular surveillance study of common β-lactamases among isolates. At DHMMC, the most common ESBL gene we identified was bla CTX-M , followed by bla OXA , bla TEM , and bla SHV (Figure 1a). Their yearly frequencies are provided in Table 1. Mathematical analysis of those frequencies over time revealed no significant differences over months or agricultural seasons. However, there were some significant differences (p-value < 0.05) in yearly frequencies (Table 1). bla SHV and bla TEM frequencies were stable over time. bla CTX-M frequencies increased after the first year in 2014 and remained stable over time (Table S3). bla OXA frequencies significantly decreased in 2016 from previous years but returned to stable in 2017 (Table S4).  We then measured the frequencies at which these genes co-occurred in each isolate population. Statistical analysis suggests a genetic linkage (or correlation) between resistance genes in isolates from DHMMC ( Table 2). Two statistical methods (Pearson's chi-square test and the phi coefficient) revealed a significant positive correlation between bla TEM and bla SHV (p-value < 0.05) and between bla CTX-M and bla OXA (p-value < 0.05). There are significant negative correlations between bla TEM and bla CTX-M (p-value < 0.05) and bla TEM and bla OXA (p-value < 0.05).
When stratified by species, the resistance genes in E. coli and K. pneumoniae isolates from the DHMMC clinical isolates show unique correlations ( Table 2). In E. coli isolates, bla CTX-M and bla OXA are positively correlated with one another (p-value < 0.05) and negatively correlated with bla TEM (p-value < 0.05). However, in K. pneumoniae, all the resistance genes are positively correlated with each other (p-value < 0.05) but significance is lost for bla TEM and bla OXA after the false discovery rate (FDR)-controlling procedure.

U.S. Database Gene Frequencies
We conducted an analogous surveillance study using a nationwide database of ESBL clinical isolates from the NIH Pathogen Detection Isolates Browser (Figure 1b). Nationwide, the resistance gene frequencies were different from the DHMMC repository. The most common ESBL gene was bla TEM , followed by bla SHV , bla CTX-M , and bla OXA . Those frequencies are provided in Table 3. Analysis of these gene frequencies over time also showed no significant differences in the frequencies of common ESBL genes in the U.S. over a period of months. However, there were significant differences (p-value < 0.05) in the yearly frequencies for bla SHV , bla TEM , and bla CTX-M (Table 3). bla SHV significantly increased in 2015, followed by a significant decrease in 2016, but returned to stable in 2017 (Table S5). bla TEM frequencies significantly decreased in 2017 and 2018 from previous years (Table S6). There was a non-significant decrease in bla CTX-M frequency in 2015 from previous years. The 2015 bla CTX-M frequency was significantly different than that of 2018 due to an increase that year (Table S7). Table 2. Linkage analysis summary for DHMMC isolates. The p-value for a chi-square test for linkage, the phi coefficient, and the associated p-value are presented for each resistance marker pair comparison. The number of isolates for each species is given in parenthesis. An asterisk (*) indicates a statistically significant comparison after the FDR-controlling procedure (q* = 0.025) for both the chi-square test and the phi coefficient. In the nationwide database, the co-occurrence of resistance genes differs from that of our local samples from DHMMC (Table 4). We observed a negative association between bla TEM and the other resistance markers (p-value < 0.05), and a positive association between bla CTX-M and bla OXA (p-value < 0.05). The negative association of bla TEM with both bla CTX-M and bla OXA was almost double that found at DHMMC for all species (Tables 2 and 4). When stratified by species, E. coli upholds these trends, but in K. pneumoniae, bla SHV is negatively correlated with the other resistance genes (p-value < 0.05) and bla CTX-M is positively correlated bla OXA (p-value < 0.05). Table 4. Linkage analysis summary for the Nationwide Database Isolates. The p-value for a chisquare test for linkage, the phi coefficient, and the associated p-value are presented for each resistance marker pair comparison. The number of isolates for each species is given in parenthesis. An asterisk (*) indicates a statistically significant comparison after the FDR-controlling procedure (q* = 0.025) for both the chi-square test and the phi coefficient.

Comparison of DHMMC and U.S. Populations
We then performed a formal statistical comparison between the frequencies of resistance genes in the DHMMC repository and the U.S. database to determine significant differences (Table 5). In the nationwide database, bla SHV and bla TEM occur more frequently (all p-values < 0.05) than at DHMMC. At DHMMC, bla CTX-M and bla OXA occur at higher frequencies than they do nationwide (all p-values < 0.05). A further breakdown of the gene frequencies by species and FDR controlling revealed no significant frequency difference in bla SHV within species between datasets. bla TEM occurs in 57-68% more E. coli isolates in the U.S. database than E. coli isolates at DHMMC (p-values < 0.05). There is no significant difference in bla TEM frequency in K. pneumoniae isolates. bla CTX-M occurs much more frequently in E. coli and K. pneumoniae isolates from DHMMC than in the U.S. database (all p-values < 0.05). bla CTX-M occurs in 47-58% more E. coli isolates and 33-52% more K. pneumoniae isolates at DHMMC than isolates from the nationwide U.S. database. This trend is similar but less drastic for bla OXA , which occurs in 26-36% more E. coli isolates and 26-44% more K. pneumoniae isolates than isolates nationwide. Table 5. Percent frequency differences between resistance markers at DHMMC and the Nationwide Database. The frequency of a resistance marker from DHMMC is denoted F M , and the frequency of a resistance marker from the National Database is denoted by F N . An asterisk (*) indicates a statistically significant comparison after the FDR-controlling procedure (q* = 0.025). The last column provides the 95% confidence interval for the percent difference for a particular resistance marker between the two datasets.

Discussion
We answered our initial question about the relative importance of selection and migration in small and large regions. The stability of resistance genes over time in a distinct community that differ from the nationwide frequencies strongly suggests that local selective pressures have a larger impact on frequencies than migration. DHMMC is unique with regards to the presence of bla CTX-M and bla TEM genes. At DHMMC, bla CTX-M occurs more frequently than in the nationwide database, while bla TEM occurs less frequently at DHMMC. The negative correlation at DHMMC between bla TEM and bla CTX-M and between bla TEM and bla OXA implies incompatibilities between bla TEM and at least one of the other genes. Since bla CTX-M and bla OXA are commonly linked with each other, it is not surprising that bla TEM is negatively associated with both of them, and genetic incompatibilities may exist for only one of those pairings. As those genetic compatibilities were not observed throughout the U.S., it is likely that they are the product of local selective pressures. This negative relationship results mainly from E. coli isolates; this relationship is not observed in K. pneumoniae isolates. This result indicates that either the genetic background of K. pneumoniae eliminates the incompatibilities of these genes, or that the antibiotic exposures of these pathogens is different from E. coli. Additionally, at DHMMC, there appears to be strong selection for bla CTX-M , which may be displacing bla TEM likely due to antagonism between these genes [22].
In terms of antimicrobial stewardship, our results suggest resistance may be modulated at a regional level, facilitating the implementation of effective strategies to limit and control selection of the antibiotic resistance genes. However, we also found evidence that resistance must be modulated differently for separate species, which may be difficult to conduct because of the environmental presence of antibiotics and the use of empiric therapies without species identification.
High bla SHV frequencies are unique to our nationwide dataset. However, when we considered E. coli and K. pneumoniae isolates separately, DHMMC and nationwide frequen-cies of bla SHV are similar. The greater number of K. pneumoniae isolates in the nationwide database likely accounts for the overall higher nationwide frequencies of bla SHV (Simpson's paradox [23]). A high bla TEM frequency is also unique to our nationwide dataset. There are 57-68% more E. coli isolates with bla TEM nationwide than isolates from DHMMC. E. coli isolates nationwide have a negative correlation between bla TEM and the other resistance genes, meaning these isolates are likely to have bla TEM and no other resistance genes. In K. pneumoniae isolates nationwide, there is a negative correlation between bla SHV and the other resistance genes, meaning these isolates are likely to have bla SHV and no other resistance genes, which is the opposite of the K. pneumoniae isolates from DHMMC. ES-BLs were initially derived from bla SHV and bla TEM , explaining the relatively high bla TEM frequencies nationwide, as bla TEM has proceeded to fixation [24]. bla SHV and bla TEM have been responsible for most ESBL infections since at least the 1980s, so it is reasonable for them to be widely distributed in a large-scale dataset [25].
We observed stable gene frequencies over time within the DHMMC population which differ from nationwide frequencies. This implies there is low migration or low survival of immigrants between populations, and that the selective forces within the DHMMC population are strong and constant. Overall, there is strong evidence that local selective pressures have a much stronger effect on the frequencies of ESBLs in local populations. This suggests that communities and specific regions have the potential to effectively manage ESBL frequencies through intentional antibiotic stewardship practices.

Hospital Isolates
Clinical patient isolates (n = 872) were collected from patients at DHMMC in Merced, CA, USA, from 2013 to 2018. These isolates were flagged as ESBL using Vitek 2 (bioMérieux, Inc. Hazelwood, MO, USA). These patient samples were collected from urine, blood, sputum, and wounds.

Molecular Methods
Genomic DNA was isolated using the boil preparation method by adding a single colony to 100 µL sterile water and boiling at 100 • C for 15 min. From this 100 µL solution, 1 µL was used in the PCR reaction for the respective genes with the primers listed in Table 6. Multiplex PCR was used to determine the presence of bla CTX-M , bla TEM , and bla OXA . Detection of bla SHV was run in a separate reaction. Each PCR reaction consisted of 1 µL of template DNA, 10 µM of each primer, and Taq 2X master mix (NEB) at a final concentration of 1X, and the reactions were run under the following conditions: initial denaturation at 94 • C for 10 min, 30 cycles of 94 • C for 40 s, 60 • C for 40 s, 72 • C for 1 min, and a final elongation at 72 • C for 7 min [26]. PCR amplicons were run out on 2% agarose gel at 100 V for 30 min and visualized using a ChemiDoc™ Touch Imaging System. (Figures S1 and S2).

U.S. Database
We obtained clinical isolate genomes from the NCBI RefSeq database [27], using the NIH Isolate Browser [28] to identify clinical isolates of E. coli and K. pneumoniae from the United States from 2013 to 2018. Using the Comprehensive Antibiotic Resistance Database (CARD) [29], we identified isolate genomes containing ESBL genes to compile an ESBL clinical database (n = 1060) using the BLAST+ program. In combination with a 98% identity cut-off to positively identify the frequency of bla TEM , bla OXA , bla CTX-M , and bla SHV , we applied an additional base pair match cutoff for each gene to limit partial gene matches. For bla TEM , we required a base pair (bp) match at or above 753 bp; for bla OXA , we required 831 bp; for bla CTX-M, we required 876 bp; and for bla SHV , we required 861 bp. The metadata for nationwide clinical isolate genomes were downloaded from the NIH Isolate Browser and included date, species, and location. The list of genome assemblies used to perform this analysis can be found in the Supplementary Materials.

Statistical Analysis
We used one-way analysis of variance (ANOVA) to compare the means of the resistance gene frequencies to identify significant differences between frequencies across months. We compared the same months across years, different months within the same year, different months across all years, and bins of 2, 3, 4, and 6 months. We tested for significant differences between the means of the resistance gene frequencies across years at DHMMC and the nationwide database using a Z-test [30]. We tested for significant differences in the proportions of a resistance marker between isolates from DHMMC and the nationwide database using a Z-test. Pairwise linkage among the resistance alleles in each of the two clinical isolate populations, DHMMC and the nationwide database, was assessed using a chi-square test [31]. The phi coefficient (PC) was used as a chi-square measure of directional deviation from the null relationship of independent assortment [32]. The PC has the desired property of accounting for sample size (often >500 in this work) and because it has a known sampling distribution, it allows us to compute significance and to form confidence intervals. Controlling for multiple statistical tests was conducted via the FDR-controlling procedure [33], a Bonferroni-type multiple testing procedure, with a false discovery control level of q* = 0.025. Only those results that remained significant after using FDR are reported as significant. All analyses were performed using the Statistics and Machine Learning Toolbox of MATLAB R2020a [34].
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/antibiotics10050498/s1. Table S1: Pairwise yearly frequency comparison for blaSHV in the DHMMC database, Table S2: Pairwise yearly frequency comparison for blaTEM in the DHMMC database, Table S3: Pairwise yearly frequency comparison for blaCTX-M in the DHMMC database, Table S4: Pairwise yearly frequency comparison for blaOXA in the DHMMC database, Table S5: Pairwise yearly frequency comparison for blaSHV in the Nationwide database, Table S6: Pairwise yearly frequency comparison for blaTEM in the Nationwide database, Table S7: Pairwise yearly frequency comparison for blaCTX-M in the Nationwide database, Table S8: Pairwise yearly frequency comparison for blaOXA in the Nationwide database, Figure S1: 2% agarose gel of blaCTX-M, blaTEM, and blaOXA PCR products, Figure S2: 2% agarose gel of blaSHV PCR products,  Data Availability Statement: All data used for this study are available through Supplemental Materials and through the NCBI Isolates Browser database.