Ticks and tick-borne pathogens represent the greatest vector-borne disease threat in the US. Blacklegged ticks are responsible for most human cases, yet the disease burden is unevenly distributed across the northern and southern US. Understanding the genetic characteristics influencing phenotypic differences in tick vectors is critical to elucidating disparities in tick-borne pathogen transmission dynamics. Applying evolutionary analyses to molecular variation in natural tick populations across ecological gradients will help identify signatures of local adaptation, which will improve control and mitigation strategies. In this study, we performed whole genome nanopore sequencing of individual (n=1) blacklegged ticks across their geographical range (Minnesota, Pennsylvania, and Texas) to evaluate genetic divergence among populations. Our integrated analyses identified genetic variants associated with numerous biological processes and molecular functions that segregated across populations. Notably, northern populations displayed genetic variants in genes linked to xenobiotic detoxification, transmembrane transport, and sulfation that may underpin key phenotypes influencing tick dispersal, host associations, and vectorial capacity. Nanopore sequencing further allowed the recovery of complete mitochondrial and commensal endosymbiont genomes. Our study provides further evidence of genetic divergence in epidemiologically relevant gene families among blacklegged tick clades. This report emphasizes the need to elucidate the genetic basis driving divergence among conspecific blacklegged tick clades in the US. 

Sample Collection

Questing ticks were collected by dragging a 1-m² cloth along natural vegetation, stopping every ~10 m to collect attached ticks, from four locations throughout the United States: Minnesota (Anoka County, MN, USA), Pennsylvania (Allegheny County, PA, USA), Texas (Polk County, TX, USA), and California (Riverside County, CA, USA). Collected ticks were preserved in 70% EtOH or RNAlater and sent to the University of Minnesota for processing. Ticks were morphologically identified using keys from Keirans and Clifford (1978), and Cooley and Kohls (1945). Ticks were rinsed with molecular grade water and transferred to Zymo DNA/RNA shield until DNA extraction.

DNA Extraction and Sequencing

DNA was extracted using a Qiagen MagAttract kit (Cat. 67563; Aarhus, Denmark) according to the manufacturer's instructions, with a final elution step extended to 1 h at 37°C. DNA was sheared by 20 passes through a standard 30-gauge insulin syringe. Library prep was performed with an SQK-LSK-114 native ligation sequencing kit (ONT, Oxford, UK). Sequencing was performed on a P2 Solo instrument and basecalled using Nvidia 4090 GPUs. Sequencing was performed over 3 days with nuclease flush and library reload every 24 h. Individual adult females were used for sequencing of I. scapularis and I. pacificus (n = 1) from each location.

Raw data from all runs were basecalled together using Dorado v0.8.1 with the “super accuracy” model dna_r10.4.1_e8.2_400bps_supv4.2.0. Read quality was assessed using Nanoq v0.10.0 (Steinig and Coin 2022).

Variant Analysis

Mitogenome Assembly, Annotation, and Phylogenetics

The mitogenomes of the three I. scapularis (MN, PA, TX) and I. pacificus (CA) individuals were extracted from their respective assemblies and characterized with MitoHiFi v3.2.2 (Uliano-Silva et al. 2023). MitoHiFi identifies mitochondrial contigs from whole genome sequencing datasets, assembles the genome, and annotates it for protein and RNA genes through external tools (e.g., minimap2, Hifiasm, MitoFinder, etc.) (Uliano-Silva et al. 2023). MitoHiFi was explicitly designed to handle long-read sequencing data. To ensure proper assembly and annotation, the consensus FASTA generated from MitoHiFi was polished using Medaka v2.0.1 (https://github.com/nanoporetech/medaka) and annotated using Mitos2 v2.1.9 (Bernt et al. 2013) with specifications for metazoan RefSeq and invertebrate genetic code. Manual annotation correction was performed to correct frameshift mutations resulting in premature stop codons within highly conserved protein-coding genes to ensure consistency with conserved mitochondrial gene structures. All frameshift mutations were found in homopolymer regions (i.e., stretches of consecutive adenine and thymine), a known limitation of nanopore sequencing. Manual correction was also performed to confirm tRNA annotations or added if any annotations were absent. Any tRNA annotation additions were performed using Geneious Prime software v2024.0.7 (https://www.geneious.com) by aligning the gene in question from reference mitogenomes with our polished mitogenome assembly using minimap2. Polished and assembled mitogenomes were visualized using OGDRAW v1.3.1 (Greiner et al. 2019).

After the mitogenomes were annotated, each mitogenome, cytochrome c oxidase subunit I (COI), and 16S gene sequences were used in phylogenetic analyses to determine their relatedness to congeneric Ixodes spp. To build the phylogenies, COI, 16S, and mitogenome sequences were downloaded from GenBank for each available North American Ixodes spp. The genus Ixodes (Acari: Ixodidae) is monophyletic and shares a most recent common ancestor with Metastriata, which resides within the same family, Ixodidae. As such, the metastriate species Amblyomma americanum was chosen as the outgroup and included in each analysis to root the phylogeny. Sequences were aligned using CLUSTALO v.1.2.3 (Thompson et al. 1994), and alignments were used in maximum likelihood analyses in RAxML v8.0.0 (Stamatakis 2014) with the GTRGAMMAI model of evolution and 1000 bootstrap replicates. The final tree was visualized in FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).