Chromosome-scale genome assembly

We present a chromosome-scale 2.47 Gbp assembly of the I. scapularis genome with differentiated X and Y pseudochromosomes. A total of 5.38 million PacBio HiFi reads were generated from high molecular weight (HMW) DNA extracted from adult male and female ticks. These long-read sequencing data were used to create an initial contig-level assembly with a sequencing coverage of 26.3X (Table S1). The resulting contig-level assembly was split into 69,985 contigs with N50 of 51.06 kb (Table S1). This assembly spanned 2.95 Gbp of sequence, 476 Mbp more than the final assembly size (Tables S2 and S3). The longer length of the initial contig-level assembly reflects the substantial allelic variation and repeat content of the I. scapularis genome. To overcome these challenges, we used the Khaper algorithm, which effectively solved a highly heterozygous diploid tea plant genome (Zhang et al, 2021), to differentiate primary contigs from allelic contig pairs generated from loci with unique haplotypes. The contigs belonging to the X and Y pseudochromosomes were determined using a read-depth strategy that compared the alignment rate of sequencing reads generated from male and female ticks. The Y pseudochromosome was assembled with contigs that were exclusively mapped by reads generated from the male tick-derived DNA. The Illumina paired-end sequencing reads were then used to identify and correct sequencing errors.

We used the HiRise pipeline (Putnam et al, 2016) to generate a chromosome-level assembly from the chromatin confirmation data to improve the contig-level assembly. The CHiCAGO (in vitro chromatin assay) and Hi-C (in situ chromatin assay) sequencing reads were generated from egg batches of three individual ticks. Each library produced 153 million reads of 2 × 151 bp length sequences. Together, these CHiCAGO library reads provided 25.72 x physical coverage of the genome. The Hi-C libraries generated an average of 142.7 million reads with a 686.8X sequencing depth. The CHiCAGO-based HiRise assembly resulted in an N50 of 419 kb. A total of 50,261 contigs, or 83.76% of the contig sequences, were successfully anchored using CHiCAGO and Hi-C analysis. The final assembly produced an impressive scaffold N50 value of 207.9 Mbp (Table 1), with pseudochromosome 1 representing the largest scaffold at 299.2 Mbp (Tables 1 and S2). A total of 15 pseudochromosomes represent the 13 autosomes and the X and Y sex pseudochromosomes (Fig 1A and B and Tables S2 and S3). Genome-wide analysis of chromatin interactions shows well-organized sequences, supporting a high-quality genome assembly (Fig 1B). The final assembly shows a significant reduction in the number of scaffolds compared with existing I. scapularis assemblies (Table S1 and Fig 2A–D). Analysis of the pseudo-sex chromosomes revealed that the 116.4 Mbp predicted X pseudochromosomes is 1.7-fold larger than the 69.3 Mbp predicted Y pseudochromosome, which is the shortest of the I. scapularis chromosomes (Table S2). These data agree with previous cytogenetics work in I. scapularis, which identified the Y chromosome as the shortest element in the karyotype and suggested that the X chromosomes are much larger than the Y (Chen et al, 1994).

The assembly quality was assessed by the proportion of identified BUSCOs (Benchmarking Universal Single-Copy Orthologs) (Simão et al, 2015). A total of 97.5% of the 2,934 Arachnida lineage BUSCOs were identified, indicating a high-quality genome assembly and annotation (Fig 2E and Table S4). Whereas the number of complete BUSCOs is slightly lower than the annotation of De et al (2023) (ASM1692078v2, 98.9%), both annotations are highly complete and represent a marked improvement from the previous I. scapularis assemblies (IscaW1 and ISE6), which contain a much larger proportion of fragmented and missing BUSCOs (Fig 2E). The cumulative Basic Local Alignment Search Tool (BLAST) hit coverage of predicted proteins against the mosquito, Aedes aegypti, and vinegar fly, Drosophila melanogaster, supports the completeness of tick genome assemblies (Fig 2F and G). Each genome assembly was assessed for the performance of RNA-Seq read mapping (Fig 2H) and assignment (Fig 2I). RNA-Seq reads from multiple I. scapularis life stages (eggs, larvae, nymphs, and adult males and females) and tissues (midgut and epidermis) were aligned to each reference genome with STAR (Dobin et al, 2013) and GeneWise (Birney et al, 2004). Read counts were computed with FeatureCounts (Liao et al, 2014) to compare the performance of each genome. In both read alignment and read assignment, IscGN outperformed all other genome assemblies (Fig 2H and I), including the De et al (2023) annotation (Fig 2I). On average, 12% and 4% more RNA-seq reads map to the IscGN gene set annotation than IscaW1 and De et al (2023), respectively (Table S5). Therefore, the IscGN gene set is considerably more complete and correct than previous versions available on VectorBase.

Transposable elements

Highly contiguous genome assemblies more accurately reveal the transposable element (TE) content of genomes. Our analysis revealed 1.41 Gbp of repetitive sequence, comprising 57.27% of the assembly (Table S6). The TE content of the individual pseudochromosomes ranged from 54.0–61.42%. The shortest chromosome, pseudochromosome Y, has the highest TE content (Fig 1 and Table S6). On the whole-genome scale, the most prevalent elements were the Gypsy LTR retrotransposons, comprising 15.16% of the genome (Fig 3 and Table S7). The Copia LTR elements are much scarcer and only constitute 0.02% of the total haploid genome, displaying differential evolutionary pressures on the two LTR superfamilies. A total of 3.21% of the identified LTR lacked a classification, supporting the hypothesis that the arthropod mobilome is more expansive than that of better characterized groups, such as vertebrates and plants, and requires further classification (Petersen et al, 2019). The I. scapularis genome contains 1.8-fold more DNA transposon sequences than retrotransposon sequences, with hAT elements prevailing at 11.6% of the genome (Table S7).

Predicted protein-coding genes

Complete and correct gene models are essential for studying tick biology. We used the MAKER genome annotation pipeline (Cantarel et al, 2008) to produce an annotation for the IscGN assembly, followed by manual curation of core gene families. A total of 33,236 predicted gene models and 35,041 transcripts were identified in the IscGN genome, comprising 675.7 kb of the total genome (Fig 4 and Table S8). The IscGN assembly formed the basis for a comprehensive quantification of transcript abundance in eight developmental stage-specific, nine midgut timecourse, and three epidermis time-course RNA-seq libraries (Bioproject numbers PRJNA856331, PRJNA1001997; Table S8).

Curation of multi-gene families

Large, multi-gene families are challenging to assemble and correctly annotate because recently duplicated genes typically share high sequence similarities or can be misclassified as alleles of a single gene. We curated genes in large multi-gene families that encode Hox cluster genes, chemosensory genes, peptidases and peptidase inhibitors, and cuticular proteins (Figs 5 and 6).

Sample collection

Fully engorged I. scapularis females were obtained from the Oklahoma State University Tick Rearing Facility, and a colony was maintained at 20°C and 98% relative humidity (RH) in the Gulia-Nuss laboratory. Three batches of eggs of ∼3,000 eggs each were collected from three individual females and immediately flash-frozen in liquid nitrogen. The samples were then shipped to Dovetail Genomics for CHiCAGO and Hi-C library preparation. For the PacBio HiFi sequencing, we used the I. scapularis Oklahoma strain reared in our laboratory at UNR. High molecular weight (HMW) DNA was extracted from either a single female or a pool of five males using an HMW extraction protocol adapted from Miller et al (1988).

CHiCAGO library preparation and sequencing

Three CHiCAGO libraries were prepared as described previously (Putnam et al, 2016). Briefly, ∼500 ng of DNA with a mean fragment length of 100 bp was reconstituted into chromatin in vitro and fixed with formaldehyde. The fixed chromatin was digested with DpnII, the 5′ overhangs filled in with biotinylated nucleotides, and the free blunt ends were ligated. After ligation, crosslinks were reversed, and the DNA was purified from the protein. Purified DNA was treated to remove biotin that was not internal to ligated fragments. The DNA was then sheared to a mean fragment size of ∼350 bp and sequencing libraries were generated using NEB Next Ultra enzymes and Illumina-compatible adapters. Biotin-containing fragments were isolated using streptavidin beads before PCR enrichment of each library. The libraries were sequenced on the Illumina HiSeq X platform using the rapid run mode.

Hi-C library preparation and sequencing

Three Hi-C libraries were prepared as described previously (Lieberman-Aiden et al, 2009). Briefly, egg tissue was treated with formaldehyde to fix the chromatin in place and the crosslinked chromatin was extracted. The chromatin samples were further prepared and sequenced as for CHiCAGO sequencing, above.

Contig-level genome assembly

The initial assembly at the contig level was generated using the hifiasm algorithm (Cheng et al, 2021). By implementing this algorithm, we successfully generated a comprehensive set of contig sequences, resulting in a total of 5.3 Gbp of sequences. Upon analyzing the contig sequences, it became apparent that they exceeded the estimated genome size by a significant margin. This discrepancy indicated the presence of a substantial number of redundant sequences within the initial assembly. To address this issue and mitigate the redundancy within the assembly, we used the Khaper algorithm (Zhang et al, 2021). This algorithm is specifically designed to identify and subsequently eliminate redundant sequences using a non-repeat k-mer-based approach. The algorithm selects primary contigs from allelic contig pairs based on their contig length, following the method described by Zhang et al (2021).

Identification of X and Y-linked contigs

Contigs associated with the Y-chromosome were identified through a read-depth-based strategy, with the objective of identifying contigs present in male-derived reads but absent in female-derived reads. This process involved mapping PacBio HiFi reads from female and male samples against the assembled contigs, facilitating the computation of genome-wide read coverage. In addition, the copy number for each contig was determined using popCNV (https://github.com/sc-zhang/popCNV), a software that uses a concept similar to a previously established method (Bickhart et al, 2012).

Briefly, the software mosdepth (Pedersen & Quinlan, 2018) was used to compute read coverages, which were subsequently adjusted using a locally estimated scatterplot smoothing strategy. This adjustment served to minimize the potentially adverse effects of regions with high or low GC content within each sliding window. The copy numbers were determined based on the ratio of the corrected read depth to the overall sequencing depth of the input fastq file. Contigs present as a single copy in male samples but not detectable in female samples were classified as Y-linked contigs.

To identify X-linked contigs, a similar read-depth-based strategy as that used for Y-linked contigs was used, but with the additional integration of comparative analysis between male and female samples. The aim was to identify contigs with higher copy numbers in female-derived reads as compared with male-derived reads. For the X-linked contigs, the copy numbers were determined based on the ratio of the corrected read depth in female samples to that in male samples. The contigs with corrected read depth ratios indicating an increased copy number in female samples but not deviating significantly from the expected single or fewer copies in male samples were categorized as X-linked contigs.

Scaffolding the genome assembly with CHiCAGO and Hi-C sequencing

The CHiCAGO and Hi-C read pairs were aligned to the de novo IscGN assembly using a modified SNAP read mapper to identify and correct errors in the scaffolds of the initial (contig-level) I. scapularis genome assembly (Zaharia et al, 2011 Preprint) (https://www.microsoft.com/en-us/research/project/snap/). The CHiCAGO read pairs that were individually mapped within the draft scaffolds were analyzed by HiRise to produce a likelihood model describing the genomic distance between the read pairs (Putnam et al, 2016). The resulting model was used to identify and break putative misjoins, score prospective joins, and make joins above a specified threshold. After scaffolding by CHiCAGO data, the Hi-C library read pairs were used for the second round of scaffolding. The scaffolds were ordered and oriented using the ALLHiC pipeline, which uses the contact frequency of loci to determine their proximity in the genome (Zhang et al, 2019). Scaffolds less than 100 kb were excluded from the analysis. After scaffolding, shotgun sequences were used to close gaps between contigs. Finally, the contigs and scaffolds that correspond to overlapping regions of the genome were merged by identifying pairs of scaffolds that exhibit both strong sequence homology and strong similarity in long-range contact patterns.

Polishing draft genome sequences

The finalized assembly of the IscGN genome consisted of 13 pseudoautosomes, in addition to the sex-determining X and Y pseudochromosome scaffolds. An iterative refinement of this chromosomal-scale genome assembly was undertaken with Pilon software (Walker et al, 2014). Before this, the Illumina sequencing reads were subjected to quality control and trimming via Trimmomatic (Bolger et al, 2014). These trimmed reads were then mapped to the preliminary draft of the genome using the Burrows-Wheeler Aligner (Li & Durbin, 2009). The subsequent alignment file was processed through SAMtools (Danecek et al, 2021) and subsequently served as the input for further refinement through Pilon software (Walker et al, 2014).

Genome annotation

Illumina paired-end RNA-Seq reads were trimmed for low-quality and adaptor sequences using Trimmomatic (Bolger et al, 2014). Subsequently, we aligned the trimmed reads to the genome assembly using HISAT2 (Kim et al, 2019). The mapped reads were provided as input to the BRAKER2 pipeline (Brůna et al, 2021). Within this pipeline, we trained the ab initio gene predictors AUGUSTUS (Stanke et al, 2006) and GeneMark-ET (Besemer & Borodovsky, 2005) to predict complete gene models. In addition, we performed a de novo transcriptome assembly using Trinity (Haas et al, 2013) and RNAseq data.

The MAKER pipeline was run in three iterations to generate a comprehensive genome annotation. First, MAKER was run with ab initio models from BRAKER2, ab initio models generated using FGENESH with a pre-trained I. scapularis model (Solovyev et al, 2006), and arthropod protein sequences from SwissProt (Boutet et al, 2016) and TrEMBL (Apweiler et al, 2004) (Taxon ID 6656). Translated proteins derived from the de novo I. scapularis transcriptome assembly were also included to provide empirical protein evidence. The resulting gene models were filtered to include those with an Annotation Edit Distance greater than 0.7, which contained a significant (E-value < 1 × 10⁻¹⁰) PFAM domain identified by HMMER (Eddy, 2011). The filtered gene models were then used to train and predict gene models using SNAP (Korf, 2004) and AUGUSTUS (Stanke et al, 2006). The MAKER pipeline was run a second time using AUGUSTUS, SNAP, FGENESH, and the SwissProt and TrEMBL arthropod protein sequences. Optimal full-length transcripts were subsequently predicted from the combined results of MAKER, AUGUSTUS, and FGENESH using the MIKADO pipeline. This set of full-length transcripts was further filtered to select complete ORFs using Transdecoder (https://github.com/TransDecoder/TransDecoder). A third round of the MAKER pipeline was then run using the MIKADO/Transdecoder protein sequences alongside the arthropod protein sequences from SwissProt and TrEMBL. The resulting gene models were filtered to retain those that were supported by at least 30% AUGUSTUS annotation, expressed at least three FPKM, included a significant (E-value < 1 × 10⁻¹⁰) PFAM domain, or had 50% BLASTP alignment coverage in either SwissProt or TrEMBL (Taxon ID 6656).

For the newly annotated genes, we adhered to the gene model nomenclature standards for ticks as outlined by VectorBase. These standards involve using three abbreviations derived from the species name, followed by an assigned number consisting of five digits. To differentiate the genes from the previous version, we made a specific modification by removing a single capital letter “W” that was previously present in the gene number. This letter “W” denoted the Wikel strain of I. scapularis. Therefore, in the current version, the gene names will be represented as ISCGN instead of ISCW, allowing for a clear distinction.

BUSCO analysis was performed to assess the completeness of gene annotations (Simão et al, 2015). BUSCO v 5.2.2 was run using the protein mode and arachnida_odb10 and arthropoda_odb10 lineage datasets. To provide gene descriptions, Automated Assignment of Human Readable Descriptions (Hallab et al, 2020) was run using independent BLASTp results from SwissProt, TrEMBL, and Anopheles gambiae. Known protein domains were identified with InterProScan (Jones et al, 2014), and Gene Ontology (GO) terms were derived from the Assignment of Human Readable Descriptions results.

Tick samples for RNAseq

Pathogen-free, unfed I. scapularis ticks were kept in an incubator at 95% RH and 20°C in our laboratory. Larval and nymphal ticks were blood-fed on mice, and adult ticks were fed on New Zealand white rabbits. All procedures were approved by the Institutional Animal Care and Use Committee at the University of Nevada, Reno (# 21-01-1118; Institutional Animal Care and Use Committee).

RNA extraction

For midgut transcriptomes, adult females were collected at different time points: unfed (UF), removed from the host 5 d after attachment (partially engorged, PE), and at 1, 7, and 14 d (D) after voluntary host drop-off (post blood meal, PBM). Ticks were surface cleaned with 70% ethanol. Whole midguts were dissected in cold PBS from three ticks and pooled (six for unfed samples) at each time point. Intact midguts were washed in cold PBS to remove blood. Once cleaned of blood, midguts were immediately transferred to a cold 1.7 ml tube containing 200 μl of Trizol and stored at −80°C until processed.

For epidermis/cuticle samples, adult females were collected at UF and PE (3 and 6 d post-attachment time points). Epidermis/cuticles were collected by removing all other tissues and scraping off the epidermal layer from the hard outer cuticle. Three females were pooled per replicate (six for unfed samples). All samples were collected in triplicates.

For developmental stages, 7 d old eggs (∼500 per pool), larvae (UF and 2 D PBM, ∼100 per sample), nymphs (UF, 2D PBM, 20 per sample), adult females (UF, 2 per sample), 1 D PBM (1 per sample), and males (UF, 5 per sample) were collected, surfaced cleaned, and processed in triplicates.

Total RNA was extracted using a Trizol reagent (Invitrogen) and a Zymo Direct-zol kit (Zymo Research).

Illumina sequencing

Total RNA samples were submitted to either Novogene or Genewiz Inc. for Illumina RNA library construction and sequencing. The mRNA-enriched and amplified library fragments were purified and checked for quality and final concentration using an Agilent 2100 Bioanalyzer with a High Sensitivity DNA Chip (Agilent Technologies). The final quantified libraries were pooled in equimolar amounts for sequencing on an Illumina HiSeq 2500 DNA sequencer, using a 150-bp paired-end sequencing flow cell with a HiSeq Reagent Kit (Illumina).

Transcriptome analysis

Illumina paired-end RNA-Seq data was acquired and retrieved using the fastq-dump tool (Leinonen et al, 2011). Pre-processing steps were used to enhance data quality, including the trimming of low-quality sequences and adapter sequences using Trimmomatic (Bolger et al, 2014). Subsequently, the trimmed reads were aligned to the reference genome of each respective species using Hisat2 (Kim et al, 2019). The alignment results were then subjected to sorting procedures using SAMtools (Danecek et al, 2021), ensuring appropriate arrangements for downstream analysis. To quantify the gene expression levels, featureCounts (Liao et al, 2014) was used to convert the aligned reads into raw read counts, representing the number of reads mapped to each annotated genomic feature. To enable comparison across samples, these raw read counts were normalized to transcripts per million values, which account for transcript length and library size variations.

Chemosensory gene sequences

Phylogenetic trees were constructed using predicted chemosensory genes. To analyze the evolutionary relationships among these genes, the corresponding peptide sequences within the orthogroups were aligned using MUSCLE v3.8.31 (Edgar, 2004). Subsequently, the CDSs were aligned onto the amino acid alignments using PAL2NAL v13 (Suyama et al, 2006), allowing for a comparison between nucleotide and protein sequences. Alignments were filtered using TRIMAL (Capella-Gutiérrez et al, 2009), applying two criteria. First, columns in the alignments were removed if gaps were present in more than 90% of the sequences (rows in the alignment matrix). Second, transcript translations were excluded if their coverage accounted for less than 30% of the total alignment length for the gene family. To estimate the gene trees, maximum likelihood (ML) analysis was performed using RAxML v7.3.0 (Stamatakis, 2006) with the general time reversible + Γ substitution model. Bootstrap replicates (n = 1,000) were used to assess the reliability of the inferred phylogenetic relationships.

Analysis of cuticular proteins

Cuticle proteins (CPs) were identified using eggNOG-mapper (Cantalapiedra et al, 2021) in conjunction with the DIAMOND aligner (Buchfink et al, 2015) to match protein sequences to their functional groups and families. The eggNOG 5.0 database (Huerta-Cepas et al, 2019) provided the Arthropoda orthology and phylogeny group. The additional database included HMMER V3.3.2 (Finn et al, 2011; Potter et al, 2018) to tag for hidden Markov models for identification of chitin-binding (PF00379: CPR) (Mistry et al, 2021) and ExPasy database for chitinase protein identification (EC3.2.1.14) (Gasteiger et al, 2003). FASTA amino acid sequences of the proteins identified as CPs from the IscGN assembly and IscPal assembly (VEuPathDB) (Amos et al, 2022; De et al, 2023) were extracted using the Biopython package (https://github.com/LittleRibosome/GN-Genome-I.-scap-.git). These sequences were analyzed for CPAPs, RR3, RR2, and RR1 extended motifs using the arthropod cuticle database (CutProtFam-Pred [uoa.gr]) (Karouzou et al, 2007; Ioannidou et al, 2014). The cutoff score was 22 or an E-value of 5 × 10⁻⁷ for the respective protein annotation. Cuticle protein annotation from the eggNOG 5.0 (Huerta-Cepas et al, 2019), HMMER, and Expasy databases were added to their respective IscGN or IscPal gene ID. Gene annotations of IscaW and IscGN assemblies were based on NCBI-BLAST+ run using Linux (https://github.com/LittleRibosome/GN-Genome-I.-scap-.git). CPs from the Apl laboratory assembly were similarly matched to the corresponding IscGN using NCBI-BLAST+. Using a cutoff identity score of 75% or bit score <300. VEuPathDB’s annotation, curation, and identifiers tool was used to add an annotation for all IscaW proteins.

Proteases and protease inhibitors

The amino acid sequence file was analyzed using NCBI-BLAST+ (Camacho et al, 2009) and the MEROPS database (Rawlings et al, 2018) “pepunit.lib” to identify peptidase and inhibitors. MEROPS peptidase and peptidase inhibitor accession numbers identified were matched to their respective peptidase/inhibitor family and clan. All matches were also checked in UniProt (UniProt Consortium, 2021) and HMMER (http://hmmer.org/) for identification of protease and protease inhibitors.

Genomic data

The following genome sequence and annotation data were downloaded from public repositories (Table S1): I. scapularis (release ASM1692078v2) from NCBI, I. scapularis (release IscaI1.0 and release IscaW1.7) from VectorBase, H. asiaticum (release ASM1333968v1), H. longicornis (release ASM1333976v1), R. microplus (release ASM1333972v1), Ixodes persulcatus (release BMI_IPER_1.0), R. sanguineous (release ASM1333969v1), Drosophila silvarum (release ASM1333974v1) from NCBI, A. aegypti (release AaegL5.3) from VectorBase, and D. melanogaster (release dmel_r6.44) from FlyBase.