Study site
The study site is a 67-acre private property in Red Oak, Virginia, approximately 32 km from the North Carolina border (Supplementary Figure 1) [32]. The site is within the Piedmont region of Virginia, which is characterized by hot, humid summers and mild winters. Our sampling period coincided with the time of year that typically receives the most precipitation. The property is heterogeneous, with natural, heavily vegetated areas and routinely mowed areas. To evaluate landscape risk factors influencing A. americanum abundance, we stratified our sampling across four distinct areas of the property, herein referred to as habitat plots, each representing unique landscape characteristics: (1) a forested area dominated by white oak (Quercus alba), pines (Pinus spp.), hickories (Carya glabra), and a bushy understory of various shrub species; (2) an open field with long and medium grass (described below); (3) an open field with short and medium grass (described below); and (4) a domicile area consisting of a mix of unpaved gravel, forested underbrush, open-canopy grass patches, and ornamental vegetation. The domicile buildings were not used as primary residences and were only visited by the property owners once or twice a month. These structures included a storage shed, a bunkhouse, and an outhouse. The study site contained a diverse vertebrate community with evidence of black bears (Ursus americanus, Linnaeus), white-tailed deer, groundhogs (Marmota monax, Linnaeus), eastern cottontails (Sylvilagus floridanus, Bachman), and wild turkeys (Meleagris gallopavo, Linnaeus). In addition, several semi-ground-dwelling bird species were observed at the site, including northern cardinals (Cardinalis cardinalis, Linnaeus), turkey vultures (Cathartes aura, Linnaeus), and American crows (Corvus brachyrhynchos, Brehm), all of which may serve as hosts for ticks.
Field sampling design
We sampled the study site over two consecutive years. In year 1, all sampling was conducted over 10 days, from 19 June to 17 July 2023. In year 2, all sampling was conducted over 8 days from 27 June to 16 July 2024. All tick sampling was conducted between 10 a.m. and 2 p.m., the period when tick activity is expected to be highest [33]. We avoided days with excessive rain that would hinder cloth sampling methods for ticks (i.e., drags and flags). At each habitat plot, we established a 5-m2 sampling grid using wooden stakes (Supplementary Figure 2). This sampling scale was selected based on the limited host-seeking movement of adult ticks, as a previous study found that they move an average of 3.2 m/day [26]. We collected spatial coordinates (i.e., longitude, latitude) using a Garmin 66st (Garmin Ltd., Olathe, KS, USA), which has reported accuracy of up to 2.5 m. The number of squares varied by habitat plot, with the short grass field (SGF) having the largest sampling area at 224 squares, followed by the long grass field (LGF) with 46 squares, and the domicile with 38 squares. The number of squares sampled in the forested plot varied between years: 11 squares were sampled in year 1, while 59 squares were sampled in year 2. We expanded forest sampling in year 2 to increase sampling effort in this plot and to assess potential edge effects observed during the previous year. In year 1, each square was sampled four times, twice by dragging and twice by flagging. The only exception was in the forest, where only two dragging events were conducted because flagging was not feasible. In year 2, all squares, including those in the forest, were sampled once by dragging and once by flagging. Sampling for a given square was never conducted on the same day, with at least one full day, often more, between resampling events to avoid bias. For each sampling event, we recorded the number of times a square had been sampled and the sampling order to identify potential biases in tick collection.
Drag cloths were made of 1-m2 cotton flannel cloth attached to a wooden dowel with a diameter of 3 cm, with two 1-m sections of rope attached to a polyvinyl chloride (PVC) pipe used as a handle, following Centers for Disease Control and Prevention guidelines [34]. We conducted drag sampling by dragging the entire 5-m2 grid square. At the end of each square, we checked the bottom and top of the drag cloth for the presence or absence of ticks, collecting and recording any. Flagging involved sweeping the flag cloth across the square, checking for ticks on both sides of the flag cloth every 1 m2. Flag cloths were made of 1-m2 cotton flannel cloth attached to a wooden pole [34]. All collected ticks were preserved in 70% ethanol and stored at −20 °C. Adult ticks were stored individually in tubes, while nymphs and larvae were pooled by grid square, with up to 25 nymphs per pool; multiple tubes were used if necessary for a single square. All ticks were identified to species using published morphological keys, except for Ixodes spp. larvae, which we only identified to genus [35,36,37,38,39,40].
To evaluate landscape risk factors contributing to A. americanum abundance, we recorded the following characteristics: habitat plot type sampled (SGF, LGF, forest, domicile), vegetation length, presence of leaf litter, presence of rocks on unpaved roads, adjacency to human structures (e.g., buildings), planted ornamental vegetation, mowed vegetation, sparse vegetation (e.g., patches of vegetation in a heterogeneous square), shade, exposed soil, adjacency to forests (e.g., forest edge), adjacency (e.g., < 5 m) to human-made roads (e.g., road edge), adjacency to open fields (e.g., field edge), dried or desiccated vegetation, and proximity to a water source (Supplementary Table 1). We selected these characteristics because they are biologically relevant factors that may influence tick abundance and serve as potential targets for intervention in tick control. We categorized vegetation length into three groups: 0–4 cm (short vegetation), 5–30 cm (medium vegetation), and 31 cm or longer (long vegetation). Although vegetation was not identified to species level, the site supported a diverse and heterogeneous plant community. We recorded the number of times each square was sampled to assess potential oversampling effects. The following model predictors were treated as either categorical or continuous variables: habitat plot, humidity, temperature, sampling method, year, and sample order (Supplementary Table 1). All other predictors in the tested models were binary (i.e., presence/absence). We recorded temperature and humidity using iButtons (iButton Link Technology, Whitewater, WI, USA) every 15 min, which were attached inside 20 plastic cups, which were inverted and placed in squares within our habitat plots. Four cups, one in each habitat plot, were placed in a randomly selected square and remained in place for the entire temperature and humidity data collection period. The remaining 16 cups were rotated weekly to randomly selected squares, with the number of cups in each habitat plot proportional to the number of squares sampled in that plot, over 3 weeks. This resulted in a total of 18,144 readings for the SGF, 6048 readings for the LGF, 10,080 readings for the forested habitat, and 6048 readings for the domicile habitat. This approach yielded average temperature and humidity readings for each habitat plot over the sampling period. The corresponding habitat plot values were assigned to each grid square for each 15-min sampling interval.
Molecular diagnostics
We extracted DNA from all A. americanum adults and nymphs to screen for R. amblyommatis, Ehrlichia chaffeensis, and E. ewingii, as these microbial agents are known to be associated with A. americanum. We completed extractions using Qiagen DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer’s recommendations with slight modifications, described below. We bisected all adult ticks longitudinally, using one half for extraction, while the other half was stored in 70% ethanol in case of extraction failure. Unlike adults, we extracted DNA from pools of up to 25 nymphs. To increase DNA yield, we incubated samples overnight in proteinase K and macerated ticks with a pestle on the second day of extraction. We stored extracted DNA at −20 °C using 60 µL of elution buffer (Qiagen, Hilden, Germany) until real-time polymerase chain reaction (qPCR) screening. To prevent contamination, we (1) sanitized materials between tick samples, (2) used new pestles for each tick maceration, and (3) completed extractions in separate areas from qPCR assays. For both E. ewingii and E. chaffeensis qPCR assays, we used the same forward primer, ECH16S-1 [41]. For E. ewingii, we used the EEW16S-97 reverse primer and EEW16S-P probe [42]. For E. chaffeensis, we used the reverse ECH16S-97 primer and ECH16S-38 probe [41]. For R. amblyommatis detection, we used Ra477F and Ra618R primers, with the Ra532P probe [43]. All probes in our study utilized the FAM (6-carboxyfluorescein) fluorescent dye for DNA quantification. We completed all qPCR assays on a CFX Connect real-time PCR detection system (Bio-Rad, Hercules, CA, USA). We included positive and negative controls in all DNA extractions and qPCR assays to validate our results.
Autocorrelation analysis
All spatial maps and autocorrelation analyses were completed in ArcGIS Pro (Esri, Redlands, CA, USA; version 3.4). To test our hypothesis of spatial autocorrelation in A. americanum abundance, we used both the optimized hotspot analysis and outlier analyses. The optimized hotspot analysis uses the Getis-Ord Gi* statistic to identify statistically significant clusters of grid squares with high tick abundance [44]. The optimized outlier analysis is similar to hotspot analysis but uses the Anselin local Moran’s I statistic to identify hot and cold spots, as well as significant spatial outliers [45]. Analyses were conducted separately for adults and nymphs. In addition to analyzing tick abundance, we also conducted these analyses using tick pathogen data to test the hypothesis that these bacteria also exhibit spatial autocorrelation at fine scales. Due to smaller individual sample sizes, we grouped E. ewingii and E. chaffeensis-positive ticks for spatial analyses as Ehrlichia spp.
Statistical modeling and analysis
We used a generalized linear mixed model (GLMM) approach to assess the landscape characteristics influencing A. americanum abundance. Our data best fit a negative binomial distribution with a log link function due to the frequency of zero collection events among the sampled squares. We first assessed the grid characteristic variables for collinearity and removed any variables that showed significant collinearity, defined as a correlation coefficient > 0.60. Using the remaining variables, we performed a stepwise variable removal to identify the best-fit model. We evaluated both year and sampling day as random effects, but year fit better as a fixed effect and was retained as an explanatory variable. To compare and assess the impacts of our explanatory variables, we scaled and zero-centered our covariates (i.e., standardized). Like our spatial analyses, we created models for both adults and nymphs. Additionally, due to significant spatial autocorrelation observed in our results, we separately tested models that included an autocorrelation explanatory variable. This was based on Moran’s I statistic from our localized hotspot analysis. Due to the strong influence of outliers on early nymph models, we applied a threshold of 10 nymphs per collection event, capping any values above 10. This cutoff was chosen based on both biological and statistical rationale. Biologically, collection events exceeding 10 nymphs were rare and often suggested unusual local conditions, such as repeated vertebrate host activity at the same grid square. Statistically, 10 nymphs represented the 99th percentile of our data (n = 29 collections), indicating that values above this point were extreme and not representative of typical collection events. Given the significance of edge habitat variables in both models, we conducted a post hoc Chi-square test of independence between edge status and infection status for both Ehrlichia spp. and R. amblyommatis to assess whether edge habitat was associated with tick infections. We conducted all data analysis using R statistical computing software (R Foundation for Statistical Computing, Vienna, Austria; version 4.3.2) and utilized the glmmTMB package for our modeling [46]. We used an alpha value of P = 0.05 for judging statistical significance.