Associations of Anaplasma phagocytophilum Bacteria Variants in Ixodes scapularis Ticks and Humans, New York, USA

Anaplasmosis, caused by the tickborne bacterium Anaplasma phagocytophilum, is an emerging public health threat in the United States. In the northeastern United States, the blacklegged tick (Ixodes scapularis) transmits the human pathogenic genetic variant of A. phagocytophilum (Ap-ha) and a nonpathogenic variant (Ap-V1). New York has recently experienced a rapid and geographically focused increase in cases of anaplasmosis. We analyzed A. phagocytophilum–infected I. scapularis ticks collected across New York during 2008–2020 to differentiate between variants and calculate an entomological risk index (ERI) for each. Ap-ha ERI varied between regions and increased in all regions during the final years of the study. Space-time scan analyses detected expanding clusters of Ap-ha located within documented anaplasmosis hotspots. Ap-ha ERI was more positively correlated with anaplasmosis incidence than non-genotyped A. phagocytophilum ERI. Our findings help elucidate the relationship between the spatial ecology of A. phagocytophilum variants and anaplasmosis.

ing 2010-2018; the median was 454 (range 220-1,112) cases/year. Anaplasmosis incidence increased statewide nearly 4-fold over a period of a decade, from 2.0 cases/100,000 persons in 2010 to 7.6 cases/100,000 persons during 2018; that increase was not spatially homogenous (12,13). Specifically, the largest increase in anaplasmosis incidence occurred in the Capital District region of New York, where incidence increased from 3.0/100,000 persons in 2008 to 5.3/100,000 persons in 2018 (13). Focal increases in incidence of anaplasmosis may result from a change in the abundance and spatial extent of Ap-ha-infected I. scapularis ticks, potentially related to changes in the deer-tick-rodent cycle (10); those changes may reflect the relative abundance of Ap-ha competent reservoir hosts, increased contact between host-seeking ticks (unfed ticks of any lifestage actively seeking a host bloodmeal) and small mammal reservoirs of Ap-ha, or a greater overlap between human residential or recreational areas and habitats conducive for the enzootic amplification of Ap-ha. Thus, further examination of the relationship between Ap-ha and Ap-V1 may broaden the understanding of anaplasmosis etiology and A. phagocytophilum ecology, refining risk assessment for this emerging disease and enabling targeted prevention efforts to reduce anaplasmosis incidence.
We analyzed A. phagocytophilum-infected hostseeking I. scapularis tick specimens to elucidate the spatial differences in entomological risk for anaplasmosis in New York. We used a TaqMan single-nucleotide polymorphism (SNP) genotyping assay (ThermoFisher Scientific, https://www.thermofisher.com) to differentiate between the Ap-ha and Ap-V1 variants as previously described (14). Next, we calculated a measure of human-infection risk as a function of I. scapularis tick density and A. phagotycophilum genotype prevalence, known as the entomological risk index (ERI) for both Ap-ha and Ap-V1.We then tested for spatiotemporal differences in counts of Ap-ha-and Ap-V1-infected I. scapularis ticks by using statistical modeling and tested for correlations between Ap-ha ERI and anaplasmosis incidence. We also used a scan statistic to search for spatiotemporal clusters of Ap-ha and Ap-V1 in New York tick populations for 2008-2020. We then compared spatiotemporal clusters of Apha and Ap-V1 in I. scapularis ticks to documented regions of increased anaplasmosis incidence in New York over the same timeframe (13,15). The results help to illuminate the relationship between the spatial ecology of each variant and the outcome of human disease.

Active Tick Sampling
We collected host-seeking ticks primarily from publicly accessible lands across New York during 2008-2020 using standardized drag-cloth or flag surveys of vegetation and forest leaf litter, as previously described (16). Generally, we visited 1 site within each county twice annually; we visited additional sites as weather and resources permitted. Collection sites had suitable tick habitat and potential risk for human exposure to ticks. We typically sampled >1,000 m 2 of suitable habitat per site during each collection event. In some instances, we did not find suitable habitat for nymphal and adult I. scapularis tick sampling at the same collection site; as a result, we sampled some sites for 1 I. scapularis lifestage, resulting in separate nymphal and adult I. scapularis tick sampling sites ( Figure 1). We stored ticks in 70%-100% ethanol at 4°C until they were sorted by developmental stage and identified to species using dichotomous keys (17,18), placed into sterile microcentrifuge tubes containing 100% ethanol, and stored at −20°C until DNA extraction (16,19).

Pathogen Detection and Ap-ha/Ap-V1 Differentiation
Individual I. scapularis nymphs and adults underwent total genomic DNA extraction as previously described (16,19). Using a quadplex real-time PCR (20), we tested for the following pathogens and target genes: A. phagocytophilum (msp 2), Babesia microti (18S rDNA), Borrelia burgdorferi (16S rDNA), and Borrelia miyamotoi (16S rDNA). Samples testing positive for A. phagocytophilum by quadplex PCR were further tested using a custom TaqMan SNP genotyping PCR to differentiate between the Ap-ha and Ap-V1 variants of A. phagocytophilum as described previously (14), with the following modifications: each 25 uL reaction contained 0.625 uL of 80× Custom TaqMan SNP Genotyping Assay primer/probe mix, 12.5 uL TaqMan Universal Master Mix (ThermoFisher), 1.875 uL nuclease-free water, and 10 uL of gDNA template (or nuclease-free water for negative controls). We performed SNP assays and variant assignment and generated allelic discrimination plots using Applied Biosystems 7500 Real-time PCR System version 2.0.5 (ThermoFisher).

ERI Calculation
To estimate risk for human exposure to Ap-ha and Ap-V1 variants of A. phagocytophilum, we calculated ERI (20) using the equation for I. scapularis nymphs and adults. A high ERI value denotes an increased risk for human exposure to a particular pathogen. Of note, Ap-V1 does not cause disease in humans, so the term entomologic risk index is used only to provide a metric for comparison between Ap-ha and Ap-V1 variants. For sites with multiple sampling events within the same season, we averaged ERIs for all events.

Bernoulli Space-Time Scan Statistic
We compared the spatiotemporal distribution of Apha-infected and Ap-V1-infected I. scapularis adults and nymphs by using the Bernoulli space-time scan statistic (21), implemented in SaTScan version 9.6 (https:// www.satscan.org) (22,23). SaTScan searched for statistically significant clusters of Ap-ha or Ap-V1; we considered a cluster to be a location at which the relative risk (RR) of the presence of a variant is >1.0 inside a given cluster compared to outside. We selected maximum spatial and temporal cluster sizes a priori for our analysis. We set maximum spatial cluster size as 10% of the collected ticks to allow for new clusters to form as time progressed and to show the movement of each variant. We selected maximum temporal cluster size as 90% of the study period to allow for the identification of established populations of either A. phagocytophilum variant.

Statistical Analysis
We tested the Spearman rank correlation between mean Ap-ha ERI in I. scapularis nymphs and adults and anaplasmosis incidence at the postal (ZIP) code tabulation area level gathered from the New York State Department of Health (NYSDOH) Communicable Disease Electronic Surveillance system as previously described (13). We assessed the correlation for each year, 2010-2018. We selected the Spearman rank test because of the underlying count data used to generate anaplasmosis incidence and Ap-ha ERI. We corrected the 18 correlation tests for multiple testing using the Bonferroni-Holm adjustment (24). We compared results of the Spearman rank tests with the results from a previous analysis using non-genotype-specific A. phagocytophilum ERI (13).
We tested for spatiotemporal interaction in the number of Ap-ha-and Ap-V1-infected I. scapularis nymphs and adults by year and across latitude and longitude categories using a generalized linear mixed model (GLMM) extension of zero-inflated negative We assessed interaction between year and latitude/ longitude categories using the likelihood ratio test (LRT). We used the natural log of the total number of I. scapularis ticks of the target developmental stage (nymphs during sampling events in late May, June, July, and August; adult ticks during April, early May, October, November, and December) as an offset. We conducted data cleaning, generation of summary statistics, and data visualization using R version 4.0.3 (http://www.rstudio.com) and the dplyr (https:// CRAN.R-project.org/package=dplyr), sf, ggplot2, and tmap R packages (27)(28)(29). We used the glmmTMB package in R (30) for modeling.

Active Tick Sampling and Pathogen Detection
We recorded active tick sampling and A. phagocytophilum genotyping results for I. scapularis tick specimens (Table 1). We categorized sampling sites according to the presence of A. phagocytophilum genetic variants (Figure 1  nymphs were determined to be positive. A. phagocytophilum genotyping deterimined that I. scapularis adults had a higher prevalence of Ap-ha (5.4%) than Ap-V1 (1.7%), whereas I. scapularis nymphs had a higher prevalence of Ap-V1 (2.6%) than Ap-ha (1.7%). We observed co-infection of Ap-ha and Ap-V1 in I. scapularis adults (0.04%).

ERI
The Hudson Valley and Capital District regions in the eastern portion of the state generally exhibited higher Ap-ha ERI than the Central and Western regions of New York (Figures 2, 3). Ap-V1 ERI of I. scapularis nymphs was generally higher in the Western and Hudson Valley regions than in the Capital District and Central NY regions, but the levels were highly variable and exhibited no obvious temporal trend. Overall, ERI for Ap-ha increased in the later years of the study period for all 4 regions.

Retrospective Bernoulli Space-Time Cluster Analysis
We detected increased RR of Ap-ha or Ap-V1 in 9 clusters of adult (

Regression Models
Anaplasmosis incidence was significantly correlated with Ap-ha ERI in I. scapularis adults for all 9 years analyzed, whereas 6 of the 9 years analyzed were correlated for I. scapularis nymphs (Table 3). Statistically significant correlation coefficients in I. scapularis adults were 0.36-0.75, an increase in the minimum and maximum correlation coefficients compared with non-variant-specific PCR results (13). Statistically significant correlation coefficients in I. scapularis nymphs were 0.19-0.68, a decrease in the minimum coefficient and an increase in the maximum coefficient compared with correlations calculated using the non-variant-specific PCR results (13). ZINB regression models of Ap-ha and Ap-V1 in I. scapularis nymphs failed to converge, likely because of an insufficient number of observations within the year, latitude, and longitude covariate combinations. We detected notable spatiotemporal interaction only in Ap-ha in I. scapularis adults ( Table 4). The final model for Ap-ha-infected I. scapularis adults indicated high interaction between both year and latitude and year and longitude in the negative binomial portion of the model; however, the LRT indicated the model with year  The results from the ZINB regression model support our variant cluster detection; Ap-ha expanded northward at an increasing rate over time, and some  westward expansion is evident. Of note, these directions are relative only to the spatial extent of New York; our analyses only examined New York data. Including data from Massachusetts, Connecticut, and Vermont could indicate a westward expansion of entomological risk from Ap-ha. If the true spread of Ap-ha is occurring radially from neighboring states east of New York, the phenomenon would appear as the northward and westward increase relative to the borders of New York, as we observed in our results. The geographic dynamics of Ap-ha and Ap-V1 are also likely linked to the deer-tick-rodent cycle and the varied reservoir competency of key vertebrate hosts. The observed geographic range expansion of Ap-ha, whereas that of Ap-V1 remained stable, indicates that the variants may not have a directly inverse relationship. Furthermore, co-infection of Ap-ha and Ap-V1 within individual ticks was rarely observed in our study and others (14), suggesting some competitive interaction between pathogen variants within the vector. The varied A. phagocytophilum genotype prevalence across I. scapularis tick life stages points to a developmental stage-specific association; the exact ecologic mechanism is unknown. One possibility is that Ap-V1 acquired during a larval tick bloodmeal may be later outcompeted by a subsequent infection of Ap-ha acquired during a nymphal bloodmeal. This phenomenon may be plausible because infection with either variant is maintained within the tick vector from one developmental stage to the next; higher prevalence of Ap-V1 in the nymphal stage did not yield a higher rate of Ap-V1 infection in adult ticks of the same cohort (31). Other possibilities include that I. scapularis larvae may be more likely to feed on deer than mice in certain geographic regions, that smallmammal populations in certain regions may not harbor Ap-ha, or that another small mammal may serve as an alternate reservoir for Ap-V1 in nature. The difference in Ap-ha prevalence between I. scapularis adults and nymphs could cause higher anaplasmosis incidence during the early spring and autumn months, when I. scapularis adults are most active, but investigating this possibility was beyond the scope of this study.
It is likely that competition between Ap-V1 and Ap-ha occurs primarily between particular Ap-V1 and Ap-ha clusters (Figures 4, 5). The region between these clusters should be a continued area of focus for epidemiologic research. Competition of variants at the local scale will likely result in spatial changes in the incidence of anaplasmosis; the  Our study showed tests for correlation using only Ap-ha in ERI calculations moderately increased the correlation between ERI and reported anaplasmosis incidence; this finding suggests that surveillance testing to detect A. phagocytophilum in host-seeking ticks must be variant-specific to yield the most accurate assessment of anaplasmosis risk. Our study had several limitations, including spatiotemporal variability in tick sampling and the limited spatial extent of our data. Collection sites are more numerous and in closer proximity to one another in the Capital District region than other regions ( Figure 1). Heterogenous sampling effort likely resulted in larger clusters of Ap-V1 than Ap-ha (Table 2; Figures 4, 5), because the a priori parameters of the cluster analysis forced clusters to include the same maximum number of ticks, regardless of the distance between sites. Therefore, Ap-V1 clusters in the Western region required a larger radius to include the same number of ticks than Ap-ha clusters in the Capital District region, potentially limiting spatial resolution in western New York. Lower sample sizes at varying latitudes and longitudes also likely reduced statistical power of the ZINB models. In addition, the directionality of the emergence of new clusters and spatiotemporal interaction in the ZINB models are limited by the lack of data outside of New York, possibly biasing results.
Given the changes in spatial distribution of the Ap-ha variant of A. phagocytophilum, we suggest medical providers in newly emergent areas familiarize themselves with the signs and symptoms of anaplasmosis to streamline prompt and accurate diagnosis and treatment to ensure the best patient outcomes. Tickborne disease prevention education campaigns should target populations along the leading edge of Ap-ha advancement in New York and elsewhere. Continued differentiation and monitoring of the Ap-ha and Ap-V1 variants of A. phagocytophilum to document rate and directionality of spread is critical; further studies will elucidate the ecologic factors driving the expansion of Ap-ha and the resulting increase in anaplasmosis. The results of our study and others can be used to educate medical practitioners and to guide public health policy and disease prevention efforts in New York.