Prevalence,Distribution,and Development of an Ecological Niche Model of Dermacentor variabilis Ticks Positive for Rickettsia montanensis

Abstract

Rickettsia montanensis has long been considered a nonpathogenic member of the spotted fever group rickettsiae. However, the infection potential of R. montanensis is being revisited in light of its recent association with a case of human infection in the United States and the possibility that additional cases may have been misdiagnosed as Rocky Mountain spotted fever. To this end, DNA was extracted from American dog ticks (Dermacentor variabilis) removed from Department of Defense (DoD) personnel and their dependents at DoD medical treatment facilities (MTFs) during 2002–2012 (n = 4792). These 4792 samples were analyzed for the presence of R. montanensis (n = 136; 2.84%) and all vector DNA was confirmed to be of D. variabilis origin using a novel Dermacentor genus-specific quantitative real-time polymerase chain reaction procedure, Derm, and a novel Dermacentor species multilocus sequence typing assay. To assess the risk of R. montanensis infection, the positive and negative samples were geographically mapped utilizing MTF site locations. Tick localities were imported into a geographical information systems (GIS) program, ArcGIS, for mapping and analysis. The ecological niche modeling (ENM) program, Maxent, was used to estimate the probability of tick presence in eastern United States using locations of both R. montanensis-positive and -negative ticks, climate, and elevation data. The ENM for R. montanensis-positive D. variabilis estimated high probabilities of the positive ticks occurring in two main areas, including the northern Midwest and mid-Atlantic portions of the northeastern regions of United States, whereas the R. montanensis-negative D. variabilis tick model showed a wider estimated range. The results suggest that R. montanensis-positive and -negative D. variabilis have different ranges where humans may be at risk and are influenced by similar and different factors.

Introduction

Spotted fever group rickettsiae (SFGR) are obligate intracellular bacteria, many of which are the etiological agents of tick-borne diseases (Hackstadt 1996). The SFGR currently include 26 known species with varying degrees of human virulence ranging from the very pathogenic Rickettsia rickettsii (Rocky Mountain spotted fever [RMSF]) to the nonpathogenic Rickettsia peacockii (Wood and Artsob 2012, Parola et al. 2013). Reports of tick-borne rickettsial diseases show an increasing trend in the number of reported cases as RMSF has increased from 1.9 cases per million in 2000 to 8.4 cases per million in 2008, with fewer in 2010 (6 cases per million) (http://cdc.gov/rmsf/stats/). In 2010, the CDC created a new reporting category, spotted fever group rickettsiosis that still includes RMSF, but encompasses reportings of non-R. rickettsii rickettsiosis cases (http://cdc.gov/rmsf/stats/).

Until recently, Rickettsia montanensis had long been considered a nonpathogenic member of the SFGR with a geographical distribution that colocalizes with R. rickettsii due to its shared association with Dermacentor variabilis (American dog tick), which ranges across eastern United States and California (Fig. 1). However, in light of its recent implication in a case of human infection (presenting as an afebrile rash with headache and nausea) in the United States (McQuiston et al. 2012) as well as the possibility that additional rickettsial cases may have been misdiagnosed as RMSF, the possibility of humans as accidental hosts for R. montanensis-infected ticks should be considered a plausible possibility. In addition, a recent PCR study utilizing blood samples taken from dogs that were seropositive for SFGR found that the causative agent was R. montanensis, instead of the more pathogenic R. rickettsii (Barrett et al. 2014).

FIG. 1.

CDC estimated distribution of Dermacentor variabilis in the United States (www.cdc.gov/ticks/geographic_distribution.html).

The high frequency at which D. variabilis parasitize humans and animals warrants investigation into the distribution of R. montanensis-positive D. variabilis as well as the risk of humans serving as accidental hosts. The distributions of tick-borne diseases are driven by the vector that transmits the disease and by vector–host movements, which for D. variabilis is complicated because several species of small animals serve as primary hosts for this tick (Azad and Beard 1998, Stein et al. 2008). Beyond the presence of the tick vector, the epidemiology of rickettsial diseases also is dependent upon the prevalence of the rickettsiae in the tick population and the interaction of the ticks with humans (Hackstadt 1996), which indicates the need to conduct surveillance of R. montanensis-positive D. variabilis.

Although R. montanensis was first reported in 1965 (Lackman et al. 1965), few studies have addressed R. montanensis-positive D. variabilis. In the 1980s and 1990s, R. montanensis was found in locations such as Connecticut (Anderson et al. 1986); Ohio (Pretzman et al. 1990); New Jersey, Oklahoma, Wisconsin, and Virginia (Stromdahl et al. 2001); and Maine (Feng et al. 1980). More recently, it was identified in Maryland (Ammerman et al. 2004), California (Dergousoff et al. 2009), North Carolina (Smith et al. 2010), Canada (Wood and Artsob 2012), Georgia (McQuiston et al. 2012), Oklahoma (Barrett et al. 2014), and Virginia (Nadolny et al. 2014). The prevalence of R. montanensis has not been well characterized and estimates range from 1.5–2.0% in Virginia in 2010/2011 to 19–19.4% in North Carolina in 2006 (Feng et al. 1980, Dergousoff et al. 2009, Smith et al. 2010, Nadolny et al. 2014). The prevalence of R. montanensis in D. variabilis removed from humans in the Department of Defense (DoD) Human Tick Test Kit Program (HTTKP) of the Army Public Health Center (APHC) ranged from 2.9% to 4.6% with an average of 3.2% across the states west of the Rocky Mountains and in California (Stromdahl et al. 2011). Since R. montanensis was considered nonpathogenic and the identification of Rickettsia in some of the earliest prevalence and distribution studies was only to genus level, the prevalence and distribution of R. montanensis-positive D. variabilis have likely been underestimated (Ammerman et al. 2004).

In studying tick-borne pathogens, variables such as pathogen expansion, often-shifting environments, and the links between disease emergence and changes in climate and habitats make identifying a simple relationship between humans and tick-borne pathogens a matter of increasing complexity (Wimberly et al. 2008b). The availability of current mapping platforms and species habitat modeling software, enhanced vector distribution reports, and the production of niche models allow for better predictions of infected vectors in locations not previously assessed. Typically, range maps are generalized over broad areas and do not account for local changes in environmental conditions such as elevation and temperature (Wimberly et al. 2008a). Recently, ecological niche modeling (ENM) has become widely used to estimate species ranges for vectors, pathogens, and/or pathogen-infected vectors (Peterson 2006, Colacicco-Mayhugh et al. 2010, Haverkost et al. 2010). ENMs examine the environmental conditions at the known locations of a species of interest and subsequently identify areas of similar environmental conditions to estimate the probability of the presence of selected species (Phillips et al. 2006).

In our study, D. variabilis localities obtained from the HTTKP were mapped and used to create ENMs estimating the distribution of R. montanensis-positive and -negative D. variabilis ticks in the United States. The results from our study increase the understanding of R. montanensis prevalence, distribution, and potential risk to humans. This was achieved by examining the distribution within the endemic range of D. variabilis with particular emphasis on those R. montanensis-positive D. variabilis associated with humans and predictive modeling (ENM). In presenting our data, R. montanensis-positive D. variabilis ticks may be referred to as positive ticks and R. montanensis-negative D. variabilis ticks may be referred to as negative ticks.

Materials and Methods

Dermacentor variabilis collection, nucleic acid preparation, and rickettsial detection and identification

R. montanensis-positive D. variabilis tick samples were obtained from the HTTKP from the years 2002–2012 (n = 28,689 for all tick species), where 4792 (16.7%) were morphologically identified as D. variabilis and prepared for rickettsial DNA assessment, as previously described (Stromdahl et al. 2001, 2011). Individual tick DNA preparations were produced and screened for pathogens using standard PCR, restriction fragment length polymorphism, and Rickettsia genus-specific qPCR as previously described (Stromdahl et al. 2001, 2011, Jiang et al. 2005, 2012). Of these, 163 were further analyzed by the Rmont qPCR assay as previously described (Smith et al. 2010).

Dermacentor genus-specific molecular identification

To confirm the identification of the R. montanensis-positive ticks as Dermacentor, a novel genus-specific qPCR assay, Derm, was utilized using primers and probes designed to target the ITS2 region DermITS2F2 (forward primer): 5′-TACCACGCAACTGCTCGT; DermITS2R2 (reverse primer): 5′-TCGACGCCCTACCACGT, and DermITS2P2 (probe): 5′-FAM-CCTCGTCTATAAGTAGCTCCGTGTTGG-BHQ1. The primers and probe oligonucleotide sequences were selected based on Dermacentor ITS2 sequences from D. andersoni, D. variabilis, and D. marginatus in GenBank (AY365361, AY365355, AY365358, AY365357, AY365360, AY365359, AY365356, AY365362, S83088, FN296275, FN296269, FN296273, and FN296278). The sequences were compared for regions of similarity, and primers and probe were designed using Beacon Design Software (Premier Biosoft, Palo Alto, CA). The qPCR assay was optimized on a SmartCycler II (Cepheid, Sunnyvale, CA) thermocycler system. The final reaction was performed in a volume of 25 μL, consisting of 12.5 μL or 0.75 U of UDG Supermix, 0.6 μM of forward and reverse primer, 0.5 μM of probe, 5 mM of MgCl₂, and 1 μL of template. The thermocycler reaction parameters consisted of a hold at 50°C for 2 min, a hold at 95°C for 2 min, and 45 cycles at 95°C for 15 s and 61°C for 30 s. Fluorescence was recorded during the annealing step and data were analyzed with SmartCycler software (Cepheid, Sunnyvale, CA).

To confirm D. variabilis speciation, a fragment of 230bp of the 16s rRNA gene was amplified from randomly selected 16 (∼10%) R. montanensis-positive D. variabilis samples using Dermacentor-specific standard PCR primers (Derm16SF [forward primer]: 5′-CAGGGAAAAACTGTTTTAAACT TAA and Uni16SR [reverse primer]: 5′-CGCTGTTATCCC TAGAGTATTT), followed by sequencing. The primers were designed by comparing Dermacentor 16S sequences from D. variabilis, D. marginatus, D. andersoni, Dermacentor occidentalis, Dermacentor albipictus, Dermacentor imitans, Dermacentor halli, Dermacentor hunteri, and Dermacentor nitans found in GenBank (AY010238, AY010235, AY010241, AY010237, AF001248, AF001233, AF001247, AF001246, AF001249, AF001232, AF001231, AF001252, L34299, L34300, AY676458, U14144, and Z97879) and utilizing Beacon Design Software (Premier Biosoft). The PCR used 0.3 μM of forward and reverse primers and was run under conditions consisting of a hold at 94°C for 3 min; 40 cycles at 94°C for 30 s, 51°C for 30 s, and 68°C for 1 min; and a hold at 72°C for 7 min. The amplicons were run on a 1.5% agarose gel and purified using the QIAquick PCR Purification Kit (Qiagen, Valenica, CA). Sequencing reactions, purification, and sequence assembly were performed as previously described (Luce-Fedrow et al. 2015).

Mapping and modeling

Locations provided by the HTTKP (ticks were accompanied by submission forms that requested where the tick bite occurred) for the R. montanensis-positive D. variabilis ticks were available for 135 of the 136 ticks (location information was not available for one sample) using either the medical treatment facility (MTF) or the reported tick collection site. Text locations were converted to latitude/longitude coordinates using Google Maps or DistancesFrom (also Google Maps based) (https://maps.google.com/ and www.distancesfrom.com/). In some instances, only the closest locality such as the county could be determined for a tick location. Using the geographic information system (GIS) program, ArcGIS (Environmental Systems Research Institute, Redlands, CA), the coordinates for the confirmed positive R. montanensis tick localities were converted into a map layer and plotted on a base map of the United States (GADM Database, www.gadm.org).

To obtain negative tick coordinates, the records of 3746 D. variabilis ticks having tested negative for R. montanensis as part of the HTTKP sampling protocol were downloaded from TickMap (www.tickmap.org). The localities of R. montanensis-negative D. variabilis ticks included localities in eastern United States, California, Washington, Colorado (n = 1), and parts of Canada. To increase the accuracy of the model, a measure of uncertainty was selected for the latitude–longitude coordinate that was less than or equal to 21,000 meters to eliminate the most inaccurately located points and mirror the highest measure of uncertainty for the positive samples (range = 2500–21,000 meters, mean = 10,533.82 meters). Duplicate presence records and records outside the eastern part of United States where too few records were found were eliminated from the study, leaving a total of 394 unique locations for modeling the negative ticks.

One kilometer resolution bioclimatic data and global elevation data were obtained from WorldClim for modeling (www.worldclim.org) (Table 1). The 20 environmental data layers (1 altitude, 12 temperature, and 7 precipitation) were reformatted in ArcGIS for use in the Maxent 3.3.3 ENM program (www.cs.princeton.edu/∼schapire/maxent/). Maxent, which utilizes a maximum entropy algorithm, was used to analyze the environmental variables and tick localities to estimate the probability of the presence over the eastern half of United States. Maxent uses a jackknife test to calculate the importance of each of the environmental variables by measuring the training gain of the variables used individually in models as well as in combination with other variables. Each variable is considered important to the model when the training gain for the variable is high when independently incorporated in the model and lowered when the variable is not incorporated in the model (Masuoka et al. 2010). The model was considered accurate when the area under the curve (AUC) was greater than 0.8 and when the minimum training presence p value was less than 0.05, indicating that the model is better than random (Phillips and Dudík 2008).

Table 1.

Ecological Niche Modeling Variables and Variable Statistics

		Rickettsia montanensis-positive Dermacentor variabilis ticks (n = 135)			R. montanenisis-negative D. variabilis ticks ^a (n = 3595)
Variable name	Variable description	Minimum	Maximum	Mean	Minimum	Maximum	Mean	Variable units	Student's t-test ^b
Alt	Elevation above sea level	5	427	135.56296	−3	1107	175.96495	Meters	7.4897E-04
Bio1	Annual mean temperature	48	153	113.17778	48	197	100.85229	°C^*10	1.6336E-06
Bio2	Mean diurnal range	92	139	117.14815	87	145	119.22726	°C^*10	1.2150E-02
Bio3	Isothermality	25	38	31.348148	25	42	30.570236	---	2.7142E-03
Bio4	Temperature seasonality	7459	12082	9007.7481	6166	12101	9495.0618	---	1.0120E-05
Bio5	Max. temperature of warmest month	265	324	297.96296	250	364	293.79332	°C^*10	9.8594E-04
Bio6	Min. temperature of coldest month	−201	0	−72.83704	−201	56	−95.16634	°C^*10	2.6055E-06
Bio7	Temperature annual range	301	477	370.8	265	477	388.95967	°C^*10	7.5242E-06
Bio8	Mean temperature of wettest quarter	16	247	198.93333	−1	269	196.95549	°C^*10	6.5542E-01
Bio9	Mean temperature of driest quarter	−118	202	26.303704	−118	268	−5.311266	°C^*10	7.1863E-05
Bio10	Mean temperature of warmest quarter	194	248	225.74815	176	279	218.60056	°C^*10	3.4985E-06
Bio11	Mean temperature of coldest quarter	−118	56	−7.081481	−118	116	−27.03255	°C^*10	1.7266E-06
Bio12	Annual precipitation	695	1356	1047.5333	604	1995	992.38081	Millimeters	7.0432E-05
Bio13	Precipitation of wettest month	95	139	112.48148	93	203	110.39833	Millimeters	1.8533E-02
Bio14	Precipitation of driest month	17	85	63.17037	15	142	54.519054	Millimeters	2.2283E-05
Bio15	Precipitation seasonality	9	53	19.718519	6	53	24.933797	---	1.0971E-04
Bio16	Precipitation of wettest quarter	271	373	317.06667	237	546	315.29708	Millimeters	3.6539E-01
Bio17	Precipitation of driest quarter	58	283	207.88889	54	455	181.03255	Millimeters	3.9882E-05
Bio18	Precipitation of warmest quarter	249	347	304.27407	183	504	304.27316	Millimeters	9.9953E-01
Bio19	Precipitation of coldest quarter	58	343	218.31852	54	520	187.65063	Millimeters	3.1845E-05

Ticks obtained from TickMap with localities with ≤21,000 meters of uncertainty.

Two-Sided, type 3 Student's t-test with (Bonferroni correction/adjustment for multiple variables) ([p ≤ α/n] = [p ≤ 0.05/20] = 0.0025).

--- Bio3, Bio4, and Bio15 are unitless bioclimatic variables (www.worldclim.org).

To compare the positive and negative tick models, the raster output probability maps from the Maxent model were imported into ArcGIS and processed with the Raster Calculator program. Using a threshold of 0.17 (based on the 10 percentile training presence of each model), images containing predicted presence and predicted absence of the vector were produced. These binary images of positive and negative D. variabilis ticks were used in Raster Calculator to produce a final raster image representing predicted areas of R. montanensis-negative D. variabilis ticks, -positive D. variabilis ticks, and -positive and -negative D. variabilis ticks.

Statistical analysis

To understand the possible reasons for differences in the estimated range of the positive versus negative ticks, the environmental conditions for each tick location were examined using statistical analysis. ArcGIS was used to extract environmental variable values at the collection sites of positive and negative ticks. Using Excel (Microsoft, Redmond, WA), a two-tailed heteroscedastic t-test with a Bonferroni adjustment/correction for multiple variables (p ≤ 0.05/20) was applied to the two data sets to determine if there was a significant difference between the positive and negative tick datasets.

Results

R. montanensis was confirmed in 136 of 4792 (2.84%) D. variabilis ticks (16.7%, 4792 of 28,689 total ticks) that were submitted to MTFs as part of the HTTKP. The DNA preparations were confirmed as R. montanensis positive by the Rmont qPCR assay, showing an average cycle threshold (Ct) value of 31.71 (range: 20.18–41.08). These positive DNA preparations were confirmed to have been isolated from Dermacentor ticks by the novel genus-specific qPCR assay, Derm, showing an average Ct value of 18.66 (range: 12.14–28.53). Furthermore, 16 randomly chosen DNA preparations from a total of 155 (136 R. montanensis-positive and 19 R. montanensis-negative D. variabilis) DNA sample preparations identified as Dermacentor by the Derm qPCR assay were amplified by the novel species-specific PCR and sequencing assay, Derm16SF/Uni16SF, yielding a 230 bp amplicon that was sequenced with Dermacentor-targeted primers to confirm 100% species identification as D. variabilis.

The locations of ticks obtained through the HTTKP were ascertained from the accompanying questionnaire form used when submitting a tick. Eighty-two percent of submitters reported that the tick bite was acquired either on a military installation or at a descriptive location (e.g., city, town, and creek). The remaining 12% did not specify a location for where the tick bite occurred. This study used the location of the MTF when a tick was submitted without a specific location reported. Using the detail provided on the questionnaire, the latitude and longitude were determined by the installation/MTF (74%), city (18.5%), county (3.0%), or descriptive location (e.g., creek, school, or land fill) (4.4%).

The distribution of R. montanensis-positive D. variabilis ticks was found to occur across the eastern half of United States, with the highest prevalence in Virginia (n = 31), Maryland (n = 23), New Jersey (n = 18), Wisconsin (n = 15), Kentucky (n = 12), and Pennsylvania (n = 10) (Fig. 2A). Virginia, followed by Maryland, contained the majority of the positive ticks. The distribution area of the positive ticks (Fig. 2A) was smaller compared with the negative ticks (Fig. 2B). Negative ticks were found to occur at least once in California, Washington, and Colorado, but primarily were found in states east of the Mississippi River, excluding Michigan and Vermont.

FIG. 2.

Distribution of D. variabilis ticks collected as part of the Human Tick Test Kit Program (HTTKP). Each point represents a location where one or more tick(s) were collected. (A) The distribution of Rickettsia montanensis-positive D. variabilis. (B) The distribution of R. montanensis-negative D. variabilis.

In developing the ENMs, a total of 62 unique site locations or presence records (each geographical coordinate was only represented once) were used in the Maxent ENM to construct the R. montanensis-positive D. variabilis model (Fig. 3A). Fifty points were randomly designated as training points to build the model and 12 points were used to test the model. The statistical evaluation indicated a high accuracy as shown by the AUC of .952 for training data and the AUC of .932 for test data and p value of 4.619E-5. The Maxent jackknife test indicated that the variable that produced the highest training gain in a model using one environmental variable was the mean temperature of the coldest quarter. The variable that decreases the gain when it is not included in the model was altitude. Altitude therefore contains information that the other variables do not. In addition to mean temperature of the coldest quarter, annual mean temperature, precipitation seasonality, and minimum temperature of coldest month were influential variables in this model.

FIG. 3.

Maxent ecological niche model of D. variabilis ticks collected as part of the HTTKP. Each point represents a spatially unique tick presence record that was incorporated in the ecological niche model; duplicate presence records were deleted. (A) Estimated probability of R. montanensis-positive D. variabilis. (B) Estimated probability of R. montanensis-negative D. variabilis.

The R. montanensis-negative D. variabilis model utilized 316 unique presence records for training and 78 for testing (Fig. 3B). The AUC was 0.949 and the p value was 1.505E-21. The most useful variable by itself was temperature seasonality, which accounts for the differences in temperature between the summer and winter seasons (Andreasson and Schmitz 2000), and the variable with the most information not in the other variables was altitude. In addition, other influential variables were mean temperature of coldest quarter, mean temperature of driest quarter, and precipitation seasonality, which reflects the changes of seasonal receipt of precipitation (Pryor and Schoof 2008).

The influential variables of the positive and negative tick models were compared. Both models indicate altitude, mean temperature of the coldest quarter, and precipitation seasonality to be the predominantly influential variables. Furthermore, the t-test p values associated with these variables are below the threshold value (p ≤ 0.0025), indicating that these conditions associated with these variables differ between the models (Table 1). The p values for the other influential variables (important to one, but not both models) listed above for each model were also below the same threshold value (Table 1). Finally, the variables that showed low significance levels by the t-test were not influential variables in either model (Table 1). Overall, the variables that statistically appear to differ between the models, but only influenced one or the other model, were temperature-based variables: annual mean temperature, temperature seasonality, minimum temperature of coldest month, and mean temperature of driest quarter.

An image representing predicted areas of R. montanensis-negative ticks, -positive ticks, and -positive and -negative ticks was created. According to the probability map produced by Maxent, the areas with higher probabilities of positive ticks were located in Delaware, Maryland, New Jersey, Virginia, and the Kentucky–Tennessee location; the areas with the lowest probability were Alabama, Georgia, Florida, and Michigan. The negative tick model shows three areas of increased probability and the positive tick models show two areas of increased probability. More specifically, the upper and lower areas of the easternmost regions found in the negative tick model are absent in the positive tick model. The positive tick model shows a low probability of positive ticks to be found in Oklahoma, Kansas, southern North Carolina, and South Carolina. Furthermore, it appears that the probability is not as great in Connecticut, Maine, and Rhode Island. To further compare the probabilities indicated by the two models, we were able to use Raster Calculator to create a prediction that incorporated both models and reinforced the probabilities each individual model reported (Fig. 4) by demonstrating that the positive only area was smaller in extent compared with the areas that may be inhabited by both positive and negative ticks.

FIG. 4.

Overlay of the results of two tick models: R. montanensis-positive and -negative D. variabilis. The purple regions on the map indicate the overlap of the models for R. montanensis-positive and -negative ticks; the yellow and red regions indicate where the two models differ: yellow represents areas predicted to only have R. montanensis-negative ticks and red represents areas predicted to only have R. montanensis-positive ticks.

Discussion

The recent discovery that R. montanensis may actually be pathogenic to humans prompted a reexamination of the prevalence and geographic distribution of R. montanensis in D. variabilis, along with a determination of the risk to humans based on models that indicate the probability of positive D. variabilis ticks in the eastern half of United States. The study investigated ticks found on military personnel and their dependents, particularly in the North East (including the mid-Atlantic), South, and much of the Midwest, as classified by the US Census Bureau (www.census.gov/geo/reference/gtc/gtc_census_divreg.html). The distribution maps developed visualized the spatial distribution of both R. montanensis-positive and -negative D. variabilis. The spatial data in conjunction with environmental variables were used to construct ENMs used to predict locations where positive ticks could be found in comparison with those locations where negative ticks could be found.

Over the 10-year tick collection period, the 2.84% prevalence of R. montanensis found in D. variabilis ticks parasitizing DoD personnel and their dependents is similar to other reportings, including field-collected ticks. Reports of other tick species being field collected include Amblyomma americanum, Amblyomma maculatum, and Ixodes scapularis (Stromdahl et al. 2011, Stromdahl and Hickling 2012). In relation to other studies, prevalence reports of R. montanensis detected in D. variabilis ticks collected include 1.5–2.0% by flagging in Virginia (Nadolny et al. 2014), 3.3% by various collection methods (dragging, flagging, CO₂-baited traps, and human collections) in Kentucky and Tennessee (military installations) (Pagac et al. 2014), 3.8% by flagging in Maryland (Ammerman et al. 2004), and 19.4% by flagging in North Carolina (Smith et al. 2010). In ticks collected in Tennessee from humans, animals, and tick drags, reports indicated 0.3% R. montanensis infections in A. americanum ticks and 10% in D. variabilis (Moncayo et al. 2010). Similarly, R. montanensis was found in D. andersoni in Montana and in D. variabilis in Ohio, New York (Long Island), and Maryland (Azad and Beard 1998). Interestingly, serological studies of dogs demonstrated that 15% of dogs in North Carolina (Breitschwerdt et al. 1987) and 50% in Oklahoma (Barrett et al. 2014) tested positive for R. montanensis. Although D. andersoni and A. americanum have shown to be infected with R. montanensis, D. variabilis is the primary vector and R. montanensis has not been found in A. maculatum, I. scapularis, D. albipictus, Ixodes affinis, Haemaphysalis leporispalustris, or Ixodes brunneus (Azad and Beard 1998, Socolovschi et al. 2009, Smith et al. 2010, Stromdahl et al. 2011, Barrett et al. 2014, Nadolny et al. 2014). In relation to the other studies, our prevalence rate falls within similar ranges. However, our collection methods used were strictly from parasitized humans and did not include environmental collections.

The localities of the tick samples were mapped and used to create ENMs for R. montanensis-positive and R. montanensis-negative D. variabilis. The ENMs provide insight into the ecological and environmental factors that a given species depends upon for survival to estimate the range of the species, as previously described (Peterson 2006). The spatial representation and resulting ENMs predicting the risk of encountering R. montanensis-positive D. variabilis indicated regions where ticks may be found. More specifically, the models included potential regions that were not part of the assessed dataset, indicating that positive ticks could inhabit more regions than what were observed by our collection distribution and could pose a risk to humans in those areas.

When modeling with Maxent, presence-only data are incorporated, and in many instances, negative or absence data are not available (Phillips et al. 2006). Utilizing our tick data that were negative for R. montanensis does not indicate an absolute absence of R. montanensis in an area; instead, it demonstrates that it has just not been yet identified there. Therefore, it was advantageous to have access to a negative dataset for modeling, where the models could be discerned individually or in comparison. The separate modeling of the R. montanensis-negative data and -positive data allowed for comparison of the pathogen and vector, particularly in regard to predicted areas and environmental requirements. The positive tick model showed clustered areas, whereas the negative tick model was more dispersed. In addition, the predicted extent/range was smaller for the positive tick model than the negative tick model. When examining the influential variables of these models, the variables that include cold temperatures, seasonal temperature changes, and altitude were important and were common to both models, although by different magnitudes.

The Maxent software produced habitat species models of the probability of presence for positive and negative D. variabilis ticks. When considering the influential variables utilized in these models, it is not surprising that temperature was shown as an important environmental condition. Previous studies have shown that ticks are sensitive to temperature extremes (Atkinson 2012) and that covariates related to Dermacentor behavior include temperature (maximum, minimum, and mean) variables (Monello and Gompper 2007). Furthermore, it has been reported that host-seeking D. variabilis prefers temperatures above 10°C (Harlan and Foster 1986) and that this behavior decreases on cooler, less clear days (Sonenshine et al. 1966). Therefore, it stands to reason that our study's prediction as to where R. montanensis-positive D. variabilis ticks are is so highly based on temperature. Altitude can affect temperature as well, but the impact of altitude in these models could be biased as MTFs are generally not located in higher altitude/mountainous regions. The temperature and time of year when the HTTKP ticks were reported were not included in our study.

As for vector causes, nature suggests that there is competition between closely related species that may limit or prevent the coexistence of one or more species, such as rickettsial organisms (Telford 2009). This is demonstrated in D. variabilis tick populations where only one SFGR was detected per tick based on collections from diverse geographical locations (Socolovschi et al. 2009). These phenomena could be the result of transovarial interference, which could potentially interfere or block a secondary infection, resulting in D. variabilis ticks only maintaining a single rickettsial species (Macaluso et al. 2002, Socolovschi et al. 2009, Telford 2009). These factors could, in part, potentially explain the smaller region shown for R. montanensis compared with that of the D. variabilis having tested negative for R. montanensis.

Published spatial studies on non-R. rickettsii SFGR, such as R. montanensis and Rickettsia parkeri, appear to be limited due to the lack of recognition of their potential importance. This is largely due to the determination that species detection by nonmolecular serological methods has not shown species specificity, indicating that many less pathogenic rickettsiae could be the cause of spotted fever illnesses (Apperson et al. 2008, Stromdahl et al. 2011, McQuiston et al. 2012, 2014, Wood and Artsob 2012, Barrett et al. 2014). Distribution maps have been generated for D. variabilis and other tick species to investigate tick behavior and/or rickettsial species (or closely related species), but at the individual state level (Clark et al. 1998, Smith et al. 2010, Stromdahl et al. 2011, Gaines et al. 2014, Nadolny et al. 2014). As for ENMs, R. montanensis models have not been published to date, although statistical and mathematical models have been produced. Kollars et al. published an ANOVA model for D. variabilis and RMSF, Harlan et al. published a model for D. variabilis questing behavior based on multivariate analysis of micrometeorologic and microenvironmental parameters, and Gaff et al. created mathematical models to investigate tick control and disease dynamics primarily in regard to Ehrlichia chaffeensis (rickettsial relative) (Harlan and Foster 1990, Kollars Jr. 1996, Gaff and Gross 2007). In addition, global habitat models and predictive risk models for Anaplasma phagocytophilum and E. chaffeensis were created using environmental and ecological variables based on logistic regression (Manangan et al. 2007, Wimberly et al. 2008a, 2008b). In 2012, Atkinson published preliminary results where a lack of spatial concordance between D. variabilis and RMSF was found (Atkinson et al. 2012). Furthermore, Atkinson's study used Maxent-based habitat suitability and disease incidence rates to evaluate spatial concordance (Atkinson et al. 2012).

Further evaluation of the model's factors determined to be influential in relation to the predicted range is warranted and could help elucidate why the range of R. montanensis is more restricted than the range of D. variabilis. Factors that influence R. montanensis-positive D. variabilis ticks that need to be ascertained and investigated include sampling/collection, environmental and topographical impact, vector characteristics (coinfection or infection suppression by endosymbionts), and host density. The data obtained by analyzing factors will be useful to help assess why certain regions are more (or less) favorable for R. montanensis-positive D. variabilis ticks. Additional studies using traditional tick collection methods (e.g., tick dragging, flagging, and CO₂ trapping) with accurate records of collection locations (GPS coordinates) over a wide area of diverse habitats will be helpful in confirming our findings. Creating a larger dataset that incorporates greater habitat and geographical diversity is important in further understanding R. montanensis and its relationship to D. variabilis migration and host movement (Sonenshine et al. 1966). Therefore, in conjunction with expanding the sampling of habitat diversity, particularly beyond primarily urbanized MTFs, collections from states that are underrepresented in the current study should be included in future datasets. These states or U.S. regions include southeastern United States (southern halves of Louisiana, Mississippi, Alabama, Georgia, and Florida), California, and Washington State, where D. variabilis has been reported as part of the HTTKP or CDC (www.cdc.gov/ticks/geographic_distribution.html). In particular, including D. variabilis from Oregon could be insightful as the northern and southern states showed evidence of this vector. Therefore, sampling in these locations would prove beneficial in the R. montanensis-positive D. variabilis distribution, allowing the incorporation of new variables (land cover could be descriptive in a non-MTF-only based study) and improving the predictive power of the ENM.

The prevalence, distribution, and ENM results of R. montanensis-positive ticks provide further insight into the evolving field of rickettsiology and into the possibility of misdiagnosed SFGR, the evaluation of fevers of unknown origin, and the better understanding of emerging infectious diseases (Paddock et al. 2004, Paddock 2005, Whitman et al. 2007, Stromdahl et al. 2011, Wood and Artsob 2012, Parola et al. 2013, Barrett et al. 2014, McQuiston et al. 2014). As underestimated tick-borne febrile and afebrile diseases are better understood, the ability to treat and prevent the related illnesses also improves, as has been the case for Tidewater Spotted Fever caused by R. parkeri found in the Gulf Coast Tick (Whitman et al. 2007). Furthermore, the findings from this and related studies could be adapted to other non-tick-borne rickettsioses, such as flea-borne spotted fever caused by Rickettsia felis (Reif and Macaluso 2009, Luce-Fedrow et al. 2014), human monocytic ehrlichiosis caused by E. chaffeensis (Paddock and Childs 2003), and rickettsialpox caused by Rickettsia akari (Paddock et al. 2004, Whitman et al. 2007).

It has been documented that determining and predicting the geographical distribution of any pathogen can be complicated by the dynamic interaction between bacteria, their vectors, and their hosts (Manangan et al. 2007). Although the information and results obtained from our study are MTF based, the potential distribution of R. montanensis provides a preliminary view of where R. montanensis may be found and potentially impact human health. Future studies need to incorporate areas, both geographical and environmental, which are not represented within this study. These sorts of studies will help to create a larger data set of unique locations and additional environmental variables that can be compared with the currently developed models. Our results will help to create a more detailed understanding of R. montanensis prevalence and distribution among people parasitized with D. variabilis in the eastern half of United States. These data will be useful in assisting with future R. montanensis surveillance and prevention programs. We acknowledge that the results of this study may be limited due to the tick sampling method and lack of spatial accuracy of the collection samples. However, and most importantly, the data provide the first study concerning the details of the distribution and potential risk areas associated with R. montanensis-positive D. variabilis found parasitizing humans.

Footnotes

Acknowledgments

The authors thank Dr. Chad Porter for his statistical counseling. The views expressed in this article are those of the authors and do not necessarily represent the official policy or position of the Department of the Navy, Department of Defense, the U.S. Government, the Uniformed Services University of the Health Sciences, or the Henry M. Jackson Foundation. This project was funded by the Department of Defense Global Emerging Infections System (GEIS), work unit 0000188M.0931.001.A0074. A.L.R. and E.Y.S. are employees of the U.S. Government and their work was prepared as part of their official duties. Title 17 U.S.C. §105 provides that “Copyright protection under this title is not available for any work of the United States Government.” Title 17 U.S.C. §101 defines a U.S. Government work as a work prepared by a military service member or employee of the U.S. Government as part of that person's official duties.

Author Disclosure Statement

No competing financial interests exist.

References

Ammerman

, Swanson

, Anderson

, Schwartz

, et al. Spotted-fever group Rickettsia in Dermacentor variabilis, Maryland. Emerg Infect Dis, 2004; 10:1478–1481.

Anderson

, Magnarelli

, Philip

, Burgdorfer

. Rickettsia rickettsii and Rickettsia montana from Ixodid ticks in Connecticut. Am J Trop Med Hyg, 1986; 35:187–191.

Andreasson

, Schmitz

. Temperature seasonality in the early middle Eocene North Atlantic region: evidence from stable isotope profiles of marine gastropod shells. Geologic Soc Am Bull, 2000; 112:628–640.

Apperson

, Engber

, Nicholson

, Mead

, et al. Tick-borne diseases in North Carolina: is “Rickettsia amblyommii” a possible cause of rickettsiosis reported as Rocky Mountain spotted fever?. Vector Borne Zoonotic Dis, 2008; 8:597–606.

Atkinson

, Sarkar

, Avina

, Schuermann

, et al. Modelling spatial concordance between Rocky Mountain spotted fever disease incidence and habitat probability of its vector Dermacentor variabilis (American dog tick). Geospat Health, 2012; 7:91–100.

Azad

, Beard

. Rickettsial pathogens and their arthropod vectors. Emerg Infect Dis, 1998; 4:179–186.

Barrett

, Little

, Shaw

. “Rickettsia amblyommii” and R. montanensis infection in dogs following natural exposure to ticks. Vector Borne Zoonotic Dis, 2014; 14:20–25.

Breitschwerdt

, Moncol

, Corbett

, MacCormack

, et al. Antibodies to spotted fever-group rickettsiae in dogs in North Carolina. Am J Vet Res, 1987; 48:1436–1440.

Clark

, Oliver

Jr. , McKechnie

, Williams

DC.

Distribution, abundance, and seasonal activities of ticks collected from rodents and vegetation in South Carolina. J Vector Ecol, 1998; 23:89–105.

10.

Colacicco-Mayhugh

, Masuoka

, Grieco

. Ecological niche model of Phlebotomus alexandri and P. papatasi (Diptera: Psychodidae) in the Middle East. Int J Health Geogr, 2010; 9:2.

11.

Dergousoff

, Gajadhar

AJA

, Chilton

. Prevalence of Rickettsia species in Canadian populations of Dermacentor andersoni and D. variabilis . Appl Environ Microbiol, 2009; 75:1786–1789.

12.

Feng

, Murray

, Burgdorfer

, Spielman

, et al. Spotted fever group rickettsiae in Dermacentor variabilis from Cape Cod, Massachusetts. Am J Trop Med Hyg, 1980; 29:691–694.

13.

Gaff

, Gross

. Modeling tick-borne disease: a metapopulation model. Bull Math Biol, 2007; 69:265–288.

14.

Gaines

, Operario

, Stroup

, Stromdahl

, et al. Ehrlichia and spotted fever group Rickettsiae surveillance in Amblyomma americanum in Virginia through use of a novel six-plex real-time PCR assay. Vector Borne Zoonotic Dis, 2014; 14:307–316.

15.

Hackstadt

The biology of rickettsiae. Infect Agents Dis, 1996; 5:127–143.

16.

Harlan

, Foster

. Temperature effects on host-seeking by larval American dog ticks, Dermacentor variabilis (Say). The Ohio J Sci, 1986; 86:210–212.

17.

Harlan

, Foster

. Micrometeorologic factors affecting field host-seeking activity of adult Dermacentor variabilis (Acari: Ixodidae). J Med Entomol, 1990; 27:471–479.

18.

Haverkost

, Gardner

, Peterson

. Predicting the distribution of a parasite using the ecological niche model, GARP. Rev Mex Biodivers, 2010; 81:895–902.

19.

Jiang

, Blair

, Olson

, Stromdahl

, et al. Development of a duplex quantitative real-time PCR assay for the detection of tick-borne spotted fever group rickettsia and Rickettsia rickettsii . Int Rev Armed Forces Med Sev, 2005; 78:174.

20.

Jiang

, Stromdahl

, Richards

. Detection of Rickettsia parkeri and Candidatus Rickettsia andeanae in Amblyomma maculatum Gulf Coast ticks collected from humans in the United States. Vector Borne Zoonotic Dis, 2012; 12:175–182.

21.

Kollars Jr

TM.

Interspecific differences between small mammals as hosts of immature Dermacentor variabilis (Acari: Ixodidae) and a model for detection of high risk areas of Rocky Mountain spotted fever. J Parasitol, 1996; 82:707–710.

22.

Lackman

, Bell

, Stoenner

, Pickens

. The rocky mountain spotted fever group of Rickettsias. Health Lab Sci, 1965; 2:135–141.

23.

Luce-Fedrow

, Macaluso

, Richards

. Growth of Rickettsia felis in Drosophila melanogaster S2 cells. Vector Borne Zoonotic Dis, 2014; 14:101–110.

24.

Luce-Fedrow

, Maina

, Otiang

, Ade

, et al. Isolation of Candidatus Rickettsia asemboensis from Ctenocephalides Fleas. Vector Borne Zoonotic Dis, 2015; 15:268–277.

25.

Macaluso

, Sonenshine

, Ceraul

, Azad

. Rickettsial infection in Dermacentor variabilis (Acari: Ixodidae) inhibits transovarial transmission of a second Rickettsia. J Med Entomol, 2002; 39:809–813.

26.

Manangan

, Schweitzer

, Nibbelink

, Yabsley

, et al. Habitat factors influencing distributions of Anaplasma phagocytophilum and Ehrlichia chaffeensis in the Mississippi Alluvial Valley. Vector Borne Zoonotic Dis, 2007; 7:563–573.

27.

Masuoka

, Klein

, Kim

, Claborn

, et al. Modeling the distribution of Culex tritaeniorhynchus to predict Japanese encephalitis distribution in the Republic of Korea. Geospat Health, 2010; 5:45–57.

28.

McQuiston

, Zemtsova

, Perniciaro

, Hutson

, et al. Afebrile spotted fever group Rickettsia infection after a bite from a Dermacentor variabilis tick infected with Rickettsia montanensis . Vector Borne Zoonotic Dis, 2012; 12:1059–1061.

29.

McQuiston

, Wiedeman

, Singleton

, Carpenter

, et al. Inadequacy of IgM antibody tests for diagnosis of Rocky Mountain Spotted Fever. Am J Trop Med Hyg, 2014; 91:767–770.

30.

Moncayo

, Cohen

, Fritzen

, Huang

, et al. Absence of Rickettsia rickettsii and ocuurence of other spotted fever group Rickettsiae in ticks from Tennessee. Am J Trop Med Hyg, 2010; 83:653–657.

31.

Monello

, Gompper

. Biotic and abiotic predictors of tick (Dermacentor variabilis) abundance and engorgement on free-ranging raccoons (Procyon lotor). Parasitology, 2007; 134:2053–2062.

32.

Nadolny

, Wright

, Sonenshine

, Hynes

, et al. Ticks and spotted fever group rickettsiae of southeastern Virginia. Ticks Tick Borne Dis, 2014; 5:53–57.

33.

Paddock

, Childs

. Ehrlichia chaffeensis: a prototypical emerging pathogen. Clin Microbiol Rev, 2003; 16:37–64.

34.

Paddock

, Sumner

, Comer

, Zaki

, et al. Rickettsia parkeri: a newly recognized cause of spotted fever rickettsiosis in the United States. Clin Infect Dis, 2004; 38:805–811.

35.

Paddock

. Rickettsia parkeri as a paradigm for multiple causes of tick-borne spotted fever in the western hemisphere. Ann NY Acad Sci, 2005; 1063:315–326.

36.

Pagac

, Miller

, Mazzei

, Nielsen

, et al. Rickettsia parkeri and Rickettsia montanensis, Kentucky and Tennessee, USA. Emerg Infect Dis, 2014; 20:1750.

37.

Parola

, Paddock

, Socolovschi

, Labruna

, et al. Update on tick-borne Rickettsioses around the world: a geographic approach. Clin Microbiol Rev, 2013; 26:657–702.

38.

Peterson

. Uses and requirements of ecological niche models and related distributional models. Biodivers Inform, 2006; 3:59–72.

39.

Phillips

, Anderson

, Schapire

. Maximum entropy modeling of species geographic distributions. Ecol Model, 2006; 190:231–259.

40.

Phillips

, Dudík

. Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography, 2008; 31:161–175.

41.

Pretzman

, Daugherty

, Poetter

, Ralph

. The distribution and dynamics of Rickettsia in the tick population of Ohio. Ann NY Acad Sci, 1990; 590:227–236.

42.

Pryor

, Schoof

. Changes in the seasonality of precipitation over the contiguous USA. J Geophys Res, 2008; 113.

43.

Reif

, Macaluso

. Ecology of Rickettsia felis: a review. J Med Entomol, 2009; 46:723–736.

44.

Smith

, Ponnusamy

, Jiang

, Ayyash

, et al. Bacterial pathogens in ixodid ticks from a Piedmont County in North Carolina: prevalence of rickettsial organisms. Vector Borne Zoonotic Dis, 2010; 10:939–952.

45.

Socolovschi

, Mediannikov

, Raoult

, Parola

. The relationship between spotted fever group Rickettsiae and Ixodid ticks. Vet Res, 2009; 40:34.

46.

Sonenshine

, Atwood

, Lamb

Jr . The ecology of ticks transmitting rocky mountain spotted fever in a study area in Virginia. Ann Entomol Soc Am, 1966; 59:1234–1262.

47.

Stein

, Waterman

, Waldon

. The effects of vegetation density and habitat disturbance on the spatial distribution of ixodid ticks (Acari: Ixodidae). Geospat Health, 2008; 2:241–252.

48.

Stromdahl

, Evans

, O'Brien

, Gutierrez

. Prevalence of infection in ticks submitted to the human tick test kit program of the U.S. Army Center for Health Promotion and Preventive Medicine. J Med Entomol, 2001; 38:67–74.

49.

Stromdahl

, Jiang

, Vince

, Richards

. Infrequency of Rickettsia rickettsii in Dermacentor variabilis removed from humans, with comments on the role of other human-biting ticks associated with spotted fever group Rickettsiae in the United States. Vector Borne Zoonotic Dis, 2011; 11:969–977.

50.

Stromdahl

, Hickling

. Beyond lyme: aetiology of tick-borne human diseases with emphasis on the south-eastern United States. Zoonoses Public Health, 2012; 59 Suppl 2:48–64.

51.

Telford

III . Status of the “East Side Hypothesis” (Transovarial Interference) 25 years later. Ann NY Acad Sci, 2009; 1166:144–150.

52.

Whitman

, Richards

, Paddock

, Tamminga

, et al. Rickettsia parkeri infection after tick bite, Virginia. Emerg Infect Dis, 2007; 13:334–336.

53.

Wimberly

, Baer

, Yabsley

. Enhanced spatial models for predicting the geographic distributions of tick-borne pathogens. Int J Health Geogr, 2008a; 7:15.

54.

Wimberly

, Yabsley

, Baer

, Dugan

, et al. Spatial heterogeneity of climate and land-cover constraints on distributions of tick-borne pathogens. Global Ecol Biogeogr, 2008b; 17:189–202.

55.

Wood

, Artsob

. Spotted fever group rickettsiae: a brief review and a Canadian perspective. Zoonoses Public Health, 2012; 59 Suppl 2:65–79.