All animal sampling was conducted in strict accordance with the recommendations in the Guide for Care and Use of Laboratory Animals of the National Institutes of Health, USA. Collection and handling of birds were approved by the U.S. Geological Survey Bird Banding Laboratory (permit #24101), the Louisiana Department of Wildlife and Fisheries, the Alabama Department of Conservation and Natural Resources, and the Institutional Animal Care and Use Committee (IACUC) at the University of Southern Mississippi (protocol #17081101) and the U.S. Geological Survey (protocol #LFT 2019-05).

Songbirds were sampled at six stopover sites along the northern Gulf of Mexico during autumn (late August to early November) in 2018 and 2019 and seven stopover sites during spring (mid-March to mid-May) in 2018 and 2019. Birds were captured using 15–20 nylon mist nets (12 or 6 x 2.6 m; 30-mm mesh) per site from up to four sites in southwest Louisiana and three sites in southern Alabama ( Table 1 ). Netting was conducted daily with the exception of inclement weather. During spring, nets were operated from local sunrise until approximately five hours after sunrise; during three days of the week, nets were reopened from three hours before local sunset until ~30 minutes before sunset. During autumn, nets were operated from local sunrise until approximately six hours after sunrise. Each captured individual was ringed with a uniquely numbered U.S. Geological Survey metal leg band, had its condition assessed, morphometrics recorded, and, when time allowed, inspected for the presence of ectoparasites. Our classification of bird species included three categories: resident, short-distance migrants, and long-distance migrants, primarily determined by their predominant location during the non-breeding period. Resident species were present throughout the year at our study sites. Short-distance migratory species were observed spending the non-breeding season primarily north of the Tropic of Cancer, while long-distance migratory species predominantly inhabited regions south of the Tropic of Cancer during the stationary non-breeding period (DeGraaf and Rappole, 1995; Carlisle et al., 2004; Zenzal et al., 2018).

When a tick was discovered attached to a bird, we collected the tick in the field by carefully detaching it with fine-tipped forceps and preserved it in a vial of 70% ethanol. For each tick sample, we recorded the collection date, unique sample code, and the bird’s unique leg band number. Ticks collected from songbirds were mostly immature developmental stages. Taxonomic keys for immature stages of exotic or uncommon ticks are not available for all species. Furthermore, engorged immature specimens are often damaged to such an extent that diagnostic morphological characters are obliterated. Therefore, individual whole ticks were molecularly identified by using 12S DNA sequencing. A fragment of each tick’s mitochondrial 12S rRNA gene (~360 bp) sequence was amplified from 485 samples with T1B and T2A primers (Beati and Keirans, 2001) obtained from Invitrogen. Dream Taq Polymerase (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA) was used to amplify tick DNA. The amplicons were bidirectionally sequenced at Eurofins Genomics (Eurofins Genomics, Louisville, KY, USA). Complementary strands were visually examined and assembled by using Benchling. The sequences were compared to homologous nucleotide fragments in GenBank using BLAST (Mukherjee et al., 2014). Sequencing data was submitted to GenBank and the accession number for individual ticks was provided ( Supplementary Table S1 ).

To analyze tick infestation data, we used generalized linear mixed models (GLMMs) via the lme4 package in the R statistical language (Bates et al., 2015; R Core Team 2022). We used a two-step process to quantify species-level variation in tick outcomes and to assess additional temporal and migratory drivers of infestation. For binary tick positivity across all birds sampled during our study (n = 17,550), we first fit a GLMM with a fixed effect of bird species and a site-level random effect to account for regional artifacts. We then had another GLMM with an additional random effect for bird species that included fixed effects of year, migratory season, and their interaction. To analyze the intensity of tick parasitism on birds with identified ticks (n = 164), we likewise fit a GLMM with Poisson errors and an observation-level random effect nested within the site random effect to account for overdispersion (Harrison, 2014). We then fit another equivalent GLMM with an additional random effect for bird species that included fixed effects of year, migratory season, and migratory category. We excluded bird species with only one tick association (n = 7 species) from these intensity GLMMs. We assessed Poisson GLMMs for overdispersion using the performance package and adjusted any post-hoc comparisons from our models using the Benjamini–Hochberg correction. We derived model fit with marginal and condition R² via the MuMIn package (Nakagawa and Schielzeth, 2013). Statistical significance of fixed effects from GLMMs were tested using a Type II analysis-of-variance (i.e., Wald χ² statistics) with the car package.

Sixty-two of the 421 genomic DNA isolated from individual ticks did not pass quality control (QC) and the 359 that passed QC were used to generate libraries according to the library generation protocol by Illumina Indexing Methodology (Kumar et al., 2022). Briefly, a two-stage PCR amplification process was used to amplify the 16S rRNA V3-V4 region, followed by a dual indexing step that assigns unique index sequences to the V3-V4 amplicons. The concentration of the resulting PCR products were determined using qPCR and equal concentration of each sample was pooled and sequenced in a single run of an Illumina MiSeq sequencing instrument using reagent kit v2 (500 cycles) with 2 X 250 bp output at the University of Mississippi Medical Centre (UMMC) Genomics Core Facility. DNA extraction controls, PCR controls, and known mock bacterial communities (ZymoBIOMICS™ Microbial Community DNA Standard, Irvine, CA, USA) were simultaneously processed alongside the tick samples. All critical steps including determination of amplicon size, and amplicon purification following each PCR step were performed as described by Kumar et al. (2022).

Unless otherwise stated, all data preprocessing was done following the video tutorial of the Quantitative Insights into Microbial Ecology (QIIME2, v2022.2) pipeline (Bolyen et al., 2019). Briefly, demultiplexed fastq files were unzipped and the forward and reverse fastq files merged into a single fastq file using Casava (v1.8). The Atacama soil microbiome pipeline was incorporated to control demultiplexed paired-end reads using the DADA2 plugin as previously described (Callahan et al., 2016). Low-quality and Chimeric sequences were trimmed; subsequent merging of paired-end reads ensured 20 nucleotide overhangs between forward and reverse reads and the chimeric seqeuences removed from the sequence table. The non-chimeric representative sequences were aligned and construction of phylogenetic tree performed using MAFFT v. 7 and FasTree v. 2.1 plugins. Operational taxonomic assignment was performed using the qiime2 feature-classifier plugin v. 7.0, which was previously trained against the SILVA 138.1 database preclustered at 99% (Quast et al., 2013). Tables representing operational taxonomic units (OTUs) and representative taxonomy were exported and used for diversity metric analysis using the Microbiome Analyst web-based interface (Dhariwal et al., 2017; Chong et al., 2020). Raw sequences were submitted to the GenBank (BioProject # PRJNA1047655).

Microbiome Analyst, a web-based interface, was used for data visualization by employing taxonomy and metadata tables generated from data processing as input files (Dhariwal et al., 2017; Chong et al., 2020). Low count and prevalence data were filtered from the OTU table by setting values of 10 and 20, respectively. A filtered abundance table was exported and used in generating histograms of bacterial abundance in Microsoft Excel 2016 (Microsoft, 2018). Network correlation maps were constructed based on the sparse correlations for compositional data (SparCC) approach (Friedman and Alm, 2012). This approach uses the log-transformed values to carry out multiple iterations, which subsequently identifies taxa outliers to the correlation parameters (Chong et al., 2020). To compare the differences in the microbiome between tick groups based on measures of distance or dissimilarity, Bray-Curtis dissimilarity was used to calculate the differences in feature abundance and ordination of the plots was visualized using Principal Coordinates Analysis (PCoA).

The statistical analysis was carried out in Microbiome Analyst (Dhariwal et al., 2017; Chong et al., 2020). Kruskal–Wallis tests followed by Dunn’s multiple comparison tests were used to compare the differences in alpha diversity between all identified ticks at the genus and species level based on the observed OTU metric and Shannon’s diversity index. Permutational multivariate analysis of variance (PERMANOVA) was used to determine significant pairwise differences in the tick microbial communities by comparing the means across different tick genus and species. Statistically significant data were represented as P <0.05.

To identify spatial hotspots of potential tick and pathogen dispersal into North America, we mapped the distributions of bird species parasitized by ticks during spring migration, when detected ticks most likely originated from Central and South American non-breeding grounds. For those parasitized migratory birds, we aggregated species shapefiles from the International Union for Conservation of Nature using the rgdal and rgeos packages in the R statistical language (Baillie et al., 2004). Aggregated avian distributions were stratified by migration and breeding seasons to illustrate dispersal capacity to stopover sites and the breeding grounds respectively, for each identified tick species (Baillie et al., 2004). Lastly, we calculated distance between spring capture site and the centroid of the breeding range for each parasitized bird species using the geosphere package, for each unique combination of site, bird species, and tick species. We then used a GLM with a Gamma distribution to identify how mean dispersal capacity during spring migration varied across bird-infesting tick species, again using Type II analysis-of-variance to assess statistical significance.

A total of 164 individual birds of 28 species were found with attached ticks during our spring and autumn sampling periods ( Supplementary Table S1 ). When considering tick parasitism status across the 17,550 birds sampled, bird species did not significantly differ in their odds of infestation (χ² = 104.85, df = 103, P = 0.43, R² _m = R² _c = 0.22), likely owing to the overall low prevalence of parasitism (<1%). After accounting for site and species variation, tick positivity varied significantly by sampling year and migratory season (interaction: χ² = 5.96, P = 0.01, R² _m = 0.19, R² _c = 0.20). The odds of birds harboring ticks were greatest in spring 2019 compared to all other sampling periods (z > 5.2, P < 0.01).

Tick intensities among parasitized birds were highly variable across host species ( Figure 1 ; χ² = 33.93, df = 20, P = 0.03), which explained 17% of the variance (R² _m) compared to the individual- and site-level random effects (R² _c = 0.54). The Yellow-breasted Chat (Icteria virens) and Common Yellowthroat (Geothlypis trichas) had the greatest tick intensities (x̄ = 18 and 11), although only one individual per species was parasitized. Most ticks were contributed by the Hooded Warber (Setophaga citrina, 31%) and Swamp Sparrow (Melospiza georgiana, 16%), whereas the least ticks were contributed by the Red-eyed Vireo (Vireo olivaceus, <1%), Scarlet Tanager (Piranga olivacea, <1%), White-throated Sparrow (Zonotrichia albicollis, <1%), Winter Wren (Troglodytes hiemalis,<1%), American Redstart (Setophaga ruticilla, <1%), and Northern Waterthrush (Parkesia noveboracensis, <1%). Most ticks (54%) were collected at the stopover sites in Grand Chenier, LA.

Our secondary Poisson GLMM (R² _m = 0.08, R² _c = 0.52) revealed no significant variation in tick intensity by migratory season (χ² = 0.28, P = 0.60) or year (χ² = 2.38, P = 0.12), although birds tended to have more ticks in 2019 (x̄ = 3.09) and during autumn migration (x̄ = 3.00). Tick intensities did vary by migration category (χ² = 7.24, P = 0.03), which was driven by short-distance migrants having greater tick intensities on average than long-distance migrants (x̄ = 4.15 vs x̄ = 2.24; z = 2.36, P = 0.05), neither of which differed significantly from residents (x̄ = 3.44; short-distance migrants: z = 0.22, P = 0.82; long-distance migrants: z = 1.82, P = 0.10). Because short- and long-distance migratory birds were more commonly sampled, these species contributed 73% and 20% of all identified ticks ( Figure 1 ). Most ticks were collected in 2019 (71%) and during spring migration (70%). Neither of our tick intensity GLMMs showed significant overdispersion (χ² ≥ 18.74, P = 1).

Of the 421 tick samples collected from songbirds, 62 did not meet quality control standards for amplification and sequencing. The remaining 359 specimens were identified as Amblyomma americanum, A. maculatum, A. calcaratum, A. coelebs, A. geayi, A. longirostre, A. nodosum, A. ovale, A. parvum, A. sabanerae, A. triste, A. varium, Haemaphysalis leporispalustris, Ixodes brunneus, I. dentatus, and I. scapularis. Four tick species accounted for 81% of all ticks collected during the study period, with A. nodosum (29%; 120/421), A. longirostre (20%; 83/421), H. leporispalustris (16%; 66/421), and A. maculatum (16%; 69/421) representing the most tick species, respectively ( Supplementary Table S1 ).

Illumina 16S rRNA sequencing of all 359 tick DNA samples produced 11,064,738 raw forward and reverse reads. The distribution of the raw reads ranges from an average of 26,957 minimum reads per sample and a maximum of 89,223 reads per sample ( Supplementary Table S2 ). Rarefaction analysis of the individual samples from a sequencing depth of 500 to 5,000 raw reads confirmed there was adequate sequence coverage relative to the number of observed features and a plateau for individual tick samples ( Supplementary Figure S1 ).

Overall, 1,416 unique bacterial OTUs were identified across each tick genera ( Supplementary Table S3 ). Proteobacteria represented the most abundant phylum across the three genera of ticks sequenced ( Supplementary Table S4 ). The abundance of the phylum Proteobacteria ranges from 82% to 98% of the total bacterial OTU ( Supplementary Figure S2A ; Supplementary Table S4 ). Other identified OTUs were either classified in the phyla Actinobacteriota, Firmicutes, Cyanobacteria, or unassigned to any known bacteria phyla. At the genus level, regardless of the tick species, Francisella accounted for the highest abundance of bacteria, ranging from 60% to 75%. Notably, Francisella was the sole bacteria genus identified in A. varium ( Figure 2 , Supplementary Figure S2B ). Apart from Francisella, ten other bacteria genera had reads exceeding 1% in at least one tick ( Figure 2 , Supplementary Figure S2B ). Rickettsia had the second highest relative abundance (11-14%) in all ticks, while Spiroplasma and Cutibacterium were present in relatively high abundance in Haemaphysalis and Ixodes ticks respectively ( Figure 2 ; Supplementary Figure S2B ). Coxiella was present in A. coelebs (5%), A. calcaratum (13%), and A. geayi (60%), while Candidatus Midichloria was present in A. maculatum (5%) I. brunneus (37%).

Ticks from the genus Amblyomma exhibited the lowest bacterial alpha diversity, whereas Haemaphysalis exhibited the highest alpha diversity (Observed OTUs: p-value: 9.6369e-13; [Kruskal-Wallis] statistic: 58.995; Shannon: p-value: 1.0478e-13; [Kruskal-Wallis] statistic: 63.504; Figure 3A ). A. coelebs had the highest bacterial alpha diversity, while A. nodosum had the least diverse bacterial communities based on observed OTUs (p-value: 1.0466e-10; [Kruskal-Wallis] statistic: 65.716) and Shannon’s diversity index (p-value: 1.1177e-11; [Kruskal-Wallis] statistic: 70.685) when compared to all other Amblyomma species ( Figure 3B ). There was no significant difference in alpha diversity observed within the individual H. leporispalustris (p-value: 0.57981; [Kruskal-Wallis] statistic: 3.7915) and between Ixodes genera (Observed OTUs: p-value: 0.57037; [Kruskal-Wallis] statistic: 2.0098; Shannon’s Index: p-value: 0.28437; [Kruskal-Wallis] statistic: 3.7958; Simpson’s Index: p-value: 0.23367; [Kruskal-Wallis] statistic: 4.2708). ( Supplementary Figure S3 ).

Our analysis of microbial communities revealed significant differences in the bacterial communities between all tick genera (Amblyomma, Ixodes, and Haemaphysalis). The NMDS analysis of community structure using the Bray-Curtis distance matrix (PERMANOVA: F = 30.04, R² = 0.23, P < 0.001, [NMDS] Stress = 0.11) revealed overlapping clusters across all tick genera ( Figure 4A ; Supplementary Figure S4A ). Several Amblyomma and Haemaphysalis ticks were seen to cluster separately on the lower axis of the PCoA ( Figure 4A ).

Our observations indicated that tick species exhibiting close clustering on the PCoA plot demonstrated similar abundance levels for both Francisella and Rickettsia ( Supplementary Figures S4B, C ). While most of the Amblyomma species tend to share similar microbial communities based on the PERMANOVA (F = 10.33, R² = 0.24, P < 0.001), several Amblyomma ticks, in particular A. maculatum, A. calcaratum, A. geayi, and A. longirostre, were outliers, indicating significantly different community structure among Amblyomma species ( Figure 4B ). The observed clustering patterns in Amblyomma ticks, including the outlier clustering of certain Amblyomma species, can be attributed to variations in the abundance of Rickettsia and Francisella ( Supplementary Figures S5A, B ). Ixodes ticks differed significantly according to the PERMANOVA of the Bray-Curtis distances (F = 2.60, R²: 0.41, P < 0.004). Unlike the Amblyomma and Ixodes ticks, only one Haemaphysalis species was part of this study, and no unique clustering was observed in this species. Similarly, the PERMANOVA test of the Bray-Curtis distances was not statistically significant (F = 1.2002, R²: 0.09525, p-value: 0.228) ( Figure 4C ). Distinct clustering was observed with I. dentatus, I. scapularis and I. sp, while I. brunneus clustered separately as explained by 18.2% (Axis 1) and 44.8% (Axis 2) of the PCoA ( Figure 4D ; Supplementary Figures S5C, D ).

Our network analysis identified 16 significant partial correlations between 12 OTUs. 87% of interactions were identified as positive correlations ( Supplementary Table S4 ). We identified negative partial correlations between the prevalence of the Francisella genus and both Rickettsia and Cutibacterium. Notably, these were the exclusive bacteria that showed interactions with Francisella. Within the network interactions, four unique clusters were identified with the largest clusters consisting of bacteria commonly reported within the tick microbiome such as Coxiella, Francisella, Rickettsia, Staphylococcus and Cutibacterium ( Figure 5 . Log-transformed counts of reads assigned to Rickettsia (p-value: 0.001022, FDR: 0.012008), Francisella (p-value: 3.7672e-11, FDR: 5.3117e-9), and Coxiella (p-value: 0.005503, FDR: 0.045643) showed significantly higher abundance in Amblyomma ticks ( Supplementary Figure S6 ). At the genus level, 65 microbe-microbe interactions were observed in the Amblyomma network, 54 in the Ixodes network, and 48 in the Haemaphysalis network. Detailed microbe-microbe interactions observed amongst each tick genus can be found in the additional file section ( Supplementary Tables S5 - S9 ).

Sixty-two percent of the total microbial interactions from the Amblyomma dataset were identified as positively correlated ( Supplementary Table S5 ). Due to the differences in the numbers from each Amblyomma species, they were divided into samples from a large sample size (> 20) and a small sample size (< 20). In the network analysis performed on the datasets with large sample sizes (A. calcaratum, A. longirostre, A. maculatum, A. nodosum, A. parvum, and A. sabanerae), we observed 56 positive and 46 negative significant partial correlations ( Supplementary Figure S7 ; Supplementary Table S6 ). Several of the identified genera belonged to environmental groups and a few bacteria were identified in other arthropod microbiomes such as Wolbachia and Spiroplasma. Spiroplasma was positively correlated with both Wolbachia and Coxiella. However, we did observe negative partial correlations between Francisella and Rickettsia similar to the whole tick dataset. Candidatus Midichloria was positively correlated with Methylobacterium Methylorubrum and Escherichia Shigella, which were both from environmental sources ( Supplementary Figure S7 ). Of the datasets with small sample sizes (A. americanum, A. coelebs, A. geayi, A. ovale, A. triste, and A. varium), 25 microbial interactions were positive significant correlations, while 10 were negatively correlated ( Supplementary Figure S8 ; Supplementary Table S7 ). From the network analysis, 5 unique clusters were identified, with the largest involving interactions between 10 bacteria genera. In contrast to the Amblyomma datasets with large representation, Francisella was not detected within any of the network clusters. In addition, Coxiella was positively correlated with Sphingomonas, while Candidatus Midichloria was negatively correlated with Bacillus.

Of the 54 significant partial interactions detected in the dataset from Ixodes ticks, 42 were positive interactions, while 12 were negative interactions ( Supplementary Table S8 ). From the network analysis, three different clusters were identified. However, we did not detect Francisella and Rickettsia in the Ixodes network dataset, indicating that other microbes, which include several of the bacteria identified as environmental groups, are more important in driving the microbial interactions. We observed positive partial correlations between Anaerobacillus, Wolbachia and Bacillus and between Wolbachia and Microbacterium ( Supplementary Figure S9 ; Supplementary Table S8 ). Candidatus Midichloria only interacted with Lawsonella and both show partial positive correlations.

In the network analysis performed on Haemaphysalis ticks, we observed 48 significant partial positive correlations and 10 significant partial negative correlations, which were grouped into three unique network clusters ( Supplementary Figure S10 ; Supplementary Table S9 ). Two of the clusters involved commonly reported environmental group bacteria with previous associations to ticks and arthropods and were linked through direct or indirect interactions with Rickettsia ( Supplementary Figure S10 ). Rickettsia was positively correlated to Spiroplasma and Candidatus Midichloria but negatively correlated with Massilia. Candidatus Midichloria was positively correlated to Stenotrophomonas, while Stenotrophomonas was also positively correlated to Coxiella ( Supplementary Figure S10 , Supplementary Table S9 ).

A total of 17,550 birds (comprising 14,929 uniquely banded individuals) were captured and sampled for ticks, encompassing 101 bird species. 1,721 individual birds were recaptured at least once. Among our entire sample, 164 individual birds, representing 28 bird species, were found to have ticks. Of the 28 bird species, we focused our spatial analyses on the 21 species of migratory birds sampled in spring for which ticks could be molecularly identified to species. Geographies of plausible tick and pathogen dispersal varied widely among the 12 tick species ( Figure 6 ). This analysis suggested migratory birds using the Gulf of Mexico during spring migration could disperse A. calcaratum, A. coelebs, A. longirostre, A. maculatum, A. nodosum, A. parvum, A. sabanerae, and A. varium possibly as far as the Northwest Territories and Alaska, whereas A. geayi, A. ovale, A. triste, and H. leporispalustris had overall less dispersal capacity (e.g., north into the Great Lakes region) owing to the shorter-distance migrations of their avian hosts. Formal analysis of distances between spring migration capture sites and breeding range centroids for parasitized birds (n = 69 unique combinations of bird species, tick species, and sites) demonstrated dispersal capacity ranged from 421 to 5003 kilometers ( Figure 7 ). On average, A. geayi had the least dispersal capacity (x̄ = 925 km), whereas A. coelebs had the most extreme dispersal capacity (x̄ = 3,331 km). A GLM excluding the two tick species with only one parasitized bird species or capture site (A. triste and A. varium) suggested such dispersal capacity among tick species to be non-significant (χ² = 15.52, P = 0.08).