Abstract
The Haemaphysalis longicornis tick, an important arthropod vector native to East Asia, has expanded its range to Oceania and the American continents, where populations are exclusively parthenogenetic females. While the parthenogenetic population possesses a stable triploid genetic system, the absence of a triploid genome assembly has prohibited us from understanding their genetic evolution and biological traits. Here, we present a chromosome-level, haplotype-resolved triploid genome assembly of the parthenogenetic population by integrating ultra-deep long-read high-fidelity, short-read high-throughput chromosome conformation capture, and isoform sequencing technologies. The phylogenetic analysis suggested that triploid genome originated from an autopolyploid event rather than interspecific hybridization. We identified functional genes that enable parthenogenetic ticks to complete embryogenesis without fertilization. Parthenogenetic ticks showed higher tolerance to Rickettsia and virus infections, with more differential expressed genes compared to bisexual ticks. They exhibited lower levels of deleterious mutations across the genome and a broader global distribution predicted by ecological niche modeling.
Similar content being viewed by others
A chromosome-scale reference genome assembly for Triplophysa lixianensis
Unzipped chromosome-level genomes reveal allopolyploid nematode origin pattern as unreduced gamete hybridization
Chromosome-level genome assembly of the Stoliczka’s Asian trident bat (Aselliscus stoliczkanus)
Introduction
The native Asian tick1, Haemaphysalis longicornis, has expanded its range into novel territories, including Australia, New Zealand, and the USA2, where it primarily reproduces through parthenogenesis. The parthenogenetic mode of reproduction enables the establishment of large populations from a single individual, thereby enhancing their invasive potential. Since its initial detection in the United States in 2017, parthenogenetic populations has rapidly spread to 19 states within a few years3, showcasing robust ecological adaptation and population expansion potential4. Haemaphysalis longicornis is known to transmit various pathogens that affect both humans and animals, including severe fever with thrombocytopenia syndrome virus (SFTSV)5 as well as Rickettsia, Anaplasma, Ehrlichia, Borrelia and Babesia6. Given the wide range of animal hosts associated with H. longicornis4, there is an elevated risk of pathogen acquisition and transmission. Therefore, intensified research and surveillance efforts on parthenogenetic H. longicornis are imperative for comprehending their biological traits, population distribution, and for preventing and controlling their threats to human and animal health.
Haemaphysalis longicornis is the only reported tick species with widely distributed parthenogenetic populations. It is generally believed that parthenogenetic species lack genetic diversity, leading to the continuous accumulation of deleterious mutations, which ultimately results in a decline of parthenogenetic populations until extinction7. However, the geographic distribution of parthenogenetic H. longicornis population continues to expand3,8, contradicting the traditional theory of parthenogenetic extinction. The parthenogenetic H. longicornis has a stable triploid genetic system9, which may result from abnormal gamete fusion mechanisms10. Triploids are generally unstable and often sterile due to the challenges of pairing and equal segregation of three homologous chromosomes during meiosis and gametogenesis10. Understanding how the triploid tick successfully overcomes cellular, genetic, and embryonic developmental constraints could provide insights into the genetic basis of both sexual and asexual reproductive modes. Additionally, genome polyploidization can significantly affect gene expression regulation11, impacting not only reproductive modes but also environmental and immune adaptation capabilities12,13, which are closely related to the vector competence of H. longicornis. The genomic and transcriptomic data are needed to understand these evolutionary adaptions.
In this study, we performed genome sequencing of the parthenogenetic population reared in the laboratory using ultra-deep long-read high-fidelity (HiFi), short-read high-throughput chromosome conformation capture (Hi-C), and long-read isoform sequencing (Iso-Seq) technologies. This approach enabled the assembly of the three homologous chromosome sets (subgenomes), followed by phylogenetic analysis to determine their evolutionary origin. A comprehensive analysis of the genomic features revealed functional gene groups linked to parthenogenetic reproduction. Transcriptome analyses of bacterial and viral infections provided insights into pathogen-vector interactions and the enhanced reservoir capacity of parthenogenetic Haemaphysalis longicornis. Additionally, we evaluated population dispersion, reproductive fitness, and predicted the global suitable habitats for parthenogenetic H. longicornis (Fig. 1a).
a The workflow of this study: ① The laboratory-reared parthenogenetic and bisexual H. longicornis populations used in this study; ② Ploidy analysis of the parthenogenetic and bisexual populations using flow cytometry; ③ Genome sequencing and assembling of the parthenogenetic population using ultra-deep high-fidelity (HiFi) sequencing, chromosome conformation capture (Hi-C), isoform sequencing (Iso-Seq), and RNA sequencing (RNA-seq) technologies; ④ Functional gene analysis related to parthenogenesis; ⑤ Transcriptome analyses after bacterial and viral infections in both parthenogenetic and bisexual H. longicornis populations; ⑥ Population dispersion and predicted the global suitable habitats of both populations. b Polyploid analysis results of parthenogenetic and bisexual ticks using flow cytometry. The parthenogenetic population shows a peak at the 112 position, while the diploid bisexual population peaks at the 74 position, indicating the triploid karotype of parthenogenetic population. c K-mer spectra and fitted models generated by GenomeScope software. The triploid genome has three major coverage peaks. “abb” indicates two nucleotides are the same in a given position of the triploid genome, and the third is different, while “abc” indicates all three nucleotides are different. d Haplotype-resolved Hi-C interaction heatmap for the 33 chromosomes across the three subgenomes (A, B, and C). Using the numbered chromosomes as coordinates, each dot’s color indicates the logarithmic value of the interaction intensity between the corresponding genomic bin pairs. The interaction intensity ranges from white (low) to red (high). “Chr” denotes chromosomes. e Example of a gene (DDX5) with multiple alternative splicing events. The Iso-Seq long reads (light red) are mapped to the annotated gene model to identify distinct alternative splicing events, indicating the improvement in the annotation of alternative splicing events. A Sashimi plot is shown to summarize the splice junctions of the gene from aligned Iso-Seq data.
Results
High-quality triploid genome
We established a laboratory parthenogenetic H. longicornis population using field-collected ticks. Adults from the tenth generation of this in-bred population were used for genome assembly and annotation. Flow cytometry was employed to compare the ploidy levels between the parthenogenetic and bisexual female ticks. The results showed fluorescence intensity peaks at positions 112 and 74 for the parthenogenetic and bisexual populations, respectively, indicating that the parthenogenetic tick had a ploidy level 1.5 times higher than that of the bisexual individual, confirming the triploid karyotype of parthenogenetic H. longicornis (Fig. 1b). A genome survey of the parthenogenetic population, conducted using second-generation sequencing, revealed three distinct k-mer peaks at different depths (Fig. 1c), indicating its triploid genome. Subsequently, we tripled the sequencing depth by increasing the HiFi coverage from 30× to 95× and the Hi-C coverage from 100× to 300×, and finally yielded 210.42 Gb of HiFi data (Supplementary Table 1) and 660 Gb of Hi-C data (Supplementary Table 2).
Using a haplotype-resolved de novo assembler, we obtained contig sequences for the parthenogenetic H. longicornis genome, resulting in a total contig length of 7936.5 Mb, a contig N50 of 11.06 Mb. Genome completeness was evaluated using BUSCO (v5.4.2) with the arthropoda_odb10 dataset, yielding a completeness score of 98.6% (Table 1). Compared to the previously assembled bisexual H. longicornis genome, which measured 2554.5 Mb with a contig N50 of 0.74 Mb and a BUSCO completeness of 91.3%14, the genome of the parthenogenetic population obtained in this study showed substantial improvements in continuity and completeness. Leveraging Hi-C sequencing data from the same population, we assembled three chromosomal haplotypes totaling 7378.94 Mb, which were mapped to 33 chromosomes, covering 98.04% of the total assembled genome sequence (Fig. 1d). Based on this assembly, we used the genome of the bisexual population as a reference to phase the chromosomal haplotypes into three subgenomes according to evolutionary relationships (Supplementary Fig. 1). The sizes of these subgenomes were 2454.2 Mb, 2438.4 Mb, and 2486.5 Mb, respectively (GenBank accession number: PRJNA1133666).
To achieve comprehensive and precise gene annotation, we employed long-read Iso-Seq on the same parthenogenetic population, generating 147,239 consensus isoforms with an average length of 1834 bp (Supplementary Table 3). By integrating these long-read transcripts into the genome annotation, we identified 28,591, 28,584 and 28,735 genes across the three subgenomes, respectively. The average coding sequence (CDS) length of parthenogenetic genes reached 1415 bp, with a 59.0% increase over the bisexual annotation. Iso-Seq analysis of the parthenogenetic tick genome significantly improved the gene models, resulting in a 5.3-fold increase in 5’ UTR elements (n = 19,851), 7.1-fold increase in 3’ UTR elements (n = 22,015), and 5.5-fold increase in alternative splicing events (n = 6430) over the annotation without Iso-Seq data (Supplementary Fig. 2), as illustrated in the example of DDX5 gene (Fig. 1e). Major alternative splicing events included exon skipping (39.1%), intron retention (31.7%), alternative 3’ splice sites (14.9%), alternative 5’ splice sites (11.4%) and mutually exclusive exon (2.9%, Supplementary Fig. 3). Furthermore, 54.9% of the parthenogenetic genome sequence was annotated as repetitive elements, including Class I/LINE (7.8%), Class I/LTR (21.0%), and Class II/DNA transposons (26.2%), with similar proportions across the three subgenomes (Supplementary Fig. 4). Compared to existing genome assemblies of bisexual H. longicornis, our assembly provided three phased subgenomes with significantly enhanced quality, supporting the future exploration of the triploid genome.
Likely origin from an autopolyploid event
To investigate the origin of parthenogenetic triploids, we constructed a phylogenetic tree based on orthologous gene families using genome data from published tick species14,15. The three subgenomes of the parthenogenetic tick formed a monophyletic clade, showing the highest similarity among themselves, while the diploid bisexual H. longicornis appeared as the closest outgroup (Fig. 2a). The evolutionary placement of bisexual H. longicornis outside of the parthenogenetic triploids suggests that the formation of the parthenogenetic triploids resulted from an autopolyploid event or involvement of closely related species16. The estimated divergence time between parthenogenetic and bisexual H. longicornis ticks is approximately 9.8 million years ago (MYA). In a triploid genome, two forms of heterozygosity were observed across the three subgenomes: one where one allele differed from the other two (abb), and another where all three alleles were different (abc). K-mer analysis revealed that the “abb” form (4.01%) in parthenogenetic H. longicornis was 3.4 times more prevalent than the “abc” form (1.18%) (Fig. 2b), and comparable to the “ab” form (3.60–4.18%) observed in the diploid genome of bisexual H. longicornis (Supplementary Fig. 5). This finding further supports the hypothesis that the parthenogenetic triploid genome arose from the duplication of a single haploid, possibly due to the fusion of unreduced gametes generated under environmental stress with haploid sperm17.
a The phylogenetic tree based on protein sequences of single-copy gene families using maximum-likelihood method. The three subgenomes of the parthenogenetic tick form a monophyletic clade, exhibiting the greatest similarity to one another. The number to the right of each node indicates the estimated divergence time (in million years ago). PHl denotes parthenogenetic H. longicornis. b Structure of heterozygous sites in the triploid genome and a hypothesis model for the origin of parthenogenetic ticks. The solid lines in the homologous chromosomes correspond to mutations, and colors correspond to distinct origins. The percentages in parentheses represent the proportion of each heterozygous form (e.g., “ab” in diploid genomes, and “abb” or “abc” in triploid genomes), estimated using Genomescope 2.0. c Chromosome synteny between the three subgenomes of the parthenogenetic ticks and the bisexual ticks. The chromosome numbers are indicated on the top. Light coral and steel blue lines represent synteny blocks involved in chromosome fission and fusion events. d Gene ontology (GO) enrichment analysis. Gene families that are the significantly expanded at the ancestral node of the three parthenogenetic subgenomes are included. The size of each circle represents the number of the enriched GO term in the database. The color denotes the p-value of the enriched GO term using Fisher’s exact test.
To evaluate the selection pressures on the three subgenomes, we calculated the nonsynonymous substitution rate (Ka), synonymous substitution rate (Ks), and the Ka/Ks ratio for each homologous gene in the subgenome compared to its homolog in the bisexual tick genome. The Ka/Ks ratio for genes in subgenome A (median Ka/Ks = 0.254) was significantly higher than those in subgenome B (median Ka/Ks = 0.243, Wilcoxon signed rank test, p = 0.002) and subgenome C (median Ka/Ks = 0.239, Wilcoxon signed rank test, p = 1.9 × 10−5, Supplementary Fig. 6), indicating that genes in subgenome A were under more relaxed purifying selection compared to their homologs in subgenomes B and C. This pattern is similar to the divergent selection pressures observed among subgenomes in allopolyploid species like plants18. Chromosome chr01 exhibited the highest Ka/Ks ratio (median Ka/Ks = 0.34, Supplementary Fig. 7) compared to other chromosomes (median Ka/Ks = 0.23, Wilcoxon rank sum test, p < 2.2 × 10⁻¹⁶). Notably, chr01 was identified as a sex chromosome (Supplementary Fig. 8), suggesting a faster evolutionary rate for the sex chromosome in the parthenogenetic population. Similarly, in studies on parthenogenetic species with typical X0 systems, such as aphids19,20, the dN/dS ratio of the X chromosome is also significantly higher than that of autosomes, confirming the phenomenon of accelerated evolution of the X chromosome.
Synteny analysis between the parthenogenetic and bisexual genomes revealed that the parthenogenetic subgenomes shared 12,738–14,344 syntenic genes, whereas they shared 7750–9962 syntenic genes with the bisexual genome (Supplementary Fig. 9), indicating a high degree of homology among the three parthenogenetic subgenomes. Remarkably, we identified chromosome fission and fusion events between the genomes of parthenogenetic and bisexual ticks. Specifically, the parthenogenetic chromosome 11 originated entirely from the fission of bisexual chromosome 2, while the parthenogenetic chromosome 7 resulted from the complete fusion of bisexual chromosomes 7 and 11 (Fig. 2c). Gene function enrichment analysis showed that genes associated with chromosome fission and fusion were primarily involved in pathways such as apoptosis, fatty acid biosynthesis, and Hedgehog signaling (Supplementary Fig. 10). These findings suggest potential roles in the development, metabolism, and unique reproductive mode of the parthenogenetic population, warranting further investigation.
To investigate the evolution of gene families following whole-genome duplication, we conducted a comparative genomic analysis between the parthenogenetic triploid and the bisexual diploid. The birth-death model of gene families revealed that the three subgenomes of the parthenogenetic tick respectively harbored 131, 123, and 151 significantly expanded gene families, as well as 138, 98, and 98 significantly contracted gene families. In contrast, the bisexual genome exhibited only 71 significantly expanded and 76 significantly contracted gene families (Supplementary Fig. 11), suggesting that the rates of gene family expansion and contraction accelerated following polyploidization, potentially enhancing the ecological adaptability21. To pinpoint major alteration of gene families during the differentiation between parthenogenetic and bisexual ticks, we further analyzed the functions of 72 expanded gene families at the ancestral nodes of the three subgenomes in the parthenogenetic tick. Gene Ontology (GO) analysis revealed that these genes were significantly enriched in the xenobiotic catabolic process, cell cycle checkpoint signaling, and positive regulation of the reproductive process (Fig. 2d), indicating that triploidization led to functional enhancements in environmental adaptation, cell cycle regulation, and reproduction.
To assess gene loss and retention among the three subgenomes (A, B, and C), we analyzed their gene content using multiple approaches. Overall, the three subgenomes exhibit strong functional complementarity. For example, KEGG Orthology (KO) annotation identified 5117 KO families across the three subgenomes, of which 74.6% were shared by all three. When considering any two subgenomes together, the coverage increased to 98.3–98.8%, suggesting that genes absent in one subgenome are often retained in another. We further analyzed the functional categories of the missing KO families. Subgenome A exhibited greater gene loss in Genetic Information Processing (n = 69), particularly in pathways related to translation and DNA repair, compared to subgenomes B (n = 44) and C (n = 47). In contrast, subgenome B showed more extensive loss in Organismal Systems (n = 55), especially in categories related to the immune and nervous systems, compared to A (n = 28) and C (n = 45). These patterns likely reflect differential retention and pseudogenization following triploidization. Similar gene loss patterns were observed in BUSCO and OrthoFinder-based gene family analyses (Supplementary Fig. 12), consistently indicating that gene loss in individual subgenomes is largely compensated by retention in the others.
Genes linked to parthenogenesis
To explore the mechanisms related to triploid reproduction and embryonic development in parthenogenetic H. longicornis, we compared the ovary transcriptional profiles from our previously published study induced by blood-feeding between parthenogenetic and bisexual H. longicornis (GenBank accession numbers: SRX11972413-SRX11972430)22, taking advantage of the parthenogenetic triploid genome assembled in this study. We firstly identified 7154 1:1:1:1 homologous gene pair between the three subgenomes and the bisexual genome. Comparative analysis revealed that the expression levels of homologous genes within each parthenogenetic subgenome were significantly lower than those in bisexual ticks at all stages of blood-feeding—early, partially engorged, and fully engorged (Fig. 3a, Supplementary Fig. 13). Among the 7154 homologous genes, the expression patterns of the three parthenogenetic subgenomes were highly similar to those of the bisexual ticks (Supplementary Fig. 14; Spearman’s rank correlation test, all P < 0.001). The reduced but similar expression profiles in the triploid parthenogenetic ticks compared to the diploid bisexual ticks suggest the presence of a dosage compensation mechanism in parthenogenetic ticks, similar to observations in tetraploid fish11.
a Expression of homologous genes of three subgenomes and bisexual H. longicornis. Each cell represents the expression level of a homolog gene in terms of transcripts per million reads (TPM). Transcriptome data are obtained from ovaries dissected at early blood-feeding early stage (48 h), partially engorged stage (72 h), and engorged stage. b Histograms of the genome-wide expression divergence of homologous genes in ovarian tissues of parthenogenetic H. longicornis. The x axis indicates the ratio of the gene expression between homologous gene pairs. N values indicate the number of dominant genes (fold change >2) in each subgenome. The p-values indicate the significance of subgenome dominance by comparing the expression levels of homologous gene pairs from the two subgenomes based on Wilcoxon rank-sum test. c Ka/Ks ratios of genes in dominant and balanced groups. The gray boxes represent the genes in balanced group. The line in the middle of each box plot represents the median of the dataset; the upper and lower edges of box plot indicate the third quartile and first quartile, respectively; and the line extending from the edge is 1.5 times the interquartile range. The p-values are estimated using the Wilcoxon rank-sum test. d Genes related to parthenogenesis in the cell cycle and oocyte meiosis pathway. TPM values of genes that are differentially expressed in parthenogenetic tick ovary but not expressed in bisexual tick during blood feeding were plotted. e Functional genes that promote stable reproductive development in the parthenogenetic H. longicornis. Genes downregulated during the transition from the early blood-feeding stage to the partially engorged stage are highlighted in blue, while those with no significant expression changes are shown in gray. G1, S, G2, and M represent the four phases of the cell cycle.
We further analyzed the expression differences among the three subgenomes at various blood-feeding stages. In both the early blood-feeding and engorged stages of ovarian tissues, the expression levels of homologous genes in subgenome A were significantly higher than those in subgenomes B and C (Fig. 3b; Wilcoxon signed-rank test, all P < 0.001), while no significant difference in expression levels was observed between subgenomes B and C. During the partially engorged stage, no significant differences were detected in pairwise comparisons (Supplementary Fig. 15). Notably, the Ka/Ks ratios of differentially expressed genes between subgenomes were significantly higher than those of non-differentially expressed genes (Fig. 3c; Wilcoxon signed-rank test, all P < 0.001). These findings suggest that differential expression and accumulation of genetic mutations among these homologous genes may facilitate neo- or sub-functionalization, ultimately driving functional innovation in gene evolution23,24.
Next, we focused on analyzing gene expression related to embryonic development in parthenogenetic and bisexual ticks. Within two key reproduction-related pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG), we discovered 27 genes in the cell cycle pathway and 13 genes in the oocyte meiosis pathway that were differentially expressed across various blood-feeding stages in parthenogenetic H. longicornis, but were not expressed in the bisexual ticks (Fig. 3d). The ovarian tissue of H. longicornis undergoes rapid development during the blood-feeding stage, and these genes may play a critical role in regulating meiosis and promoting the transition of oocytes into embryos. Based on the functional research background of these 40 genes, we primarily categorized them into the cell cycle involved to cell cycle regulation, DNA replication, chromosomal segregation or dislocation, and the oocyte meiosis involved to meiotic exit or arrest (Supplementary Fig. 16). Notably, five genes (CCNA, CCNB3, CDC20, APC12, APC13) that promote cell cycle progression showed a specifically downregulated expression from the early blood-feeding to partially engorged stage (Fig. 3e). During this stage, the ovarian cells rapidly proliferate, requiring a substantial synthesis of these proteins25. This regulatory change in expression may represent an adaptive mechanism shaped by long-term parthenogenesis, contributing to the reproductive development and stable inheritance of triploid lineages26,27. While the evolutionary origin of triploidy may involve historical chromosomal events, the dynamic expression of these cell cycle and meiosis-related genes likely reflects ongoing adaptations to maintain genomic stability and embryonic viability in the exclusively parthenogenetic reproduction of present-day H. longicornis.
Higher vector capacity
To investigate the impact of polyploid evolution on the potential transmission capability for pathogens, we inoculated Rickettsia heilongjiangensis and severe fever with thrombocytopenia syndrome virus (SFTSV) into parthenogenetic and bisexual H. longicornis, respectively. R. heilongjiangensis exhibited significant higher proliferation levels in parthenogenetic than bisexual ticks throughout the observation period (Fig. 4a, repeated measures ANOVA, P < 0.001). SFTSV proliferation was significantly higher in parthenogenetic ticks during the first 20 days (repeated measures ANOVA, P < 0.001), although there was no significant difference in SFTSV loads between parthenogenetic and bisexual ticks 25 days post-infection. This indicates that Rickettsia and SFTSV have a stronger replication capacity in parthenogenetic H. longicornis.
a Replication curves of Rickettsia heilongjiangensis and severe fever with thrombocytopenia syndrome virus (SFTSV) in the parthenogenetic and bisexual ticks. The replications of both pathogens are measured at 5d, 10d, 15d, 20d and 25d post-injection, with three biological replicates at each time point. Shaded region indicates mean ± standard deviation. b Histograms of the genome-wide expression divergence of homologous genes in whole body of parthenogenetic H. longicornis 15 days after pathogen infections. The x axis indicates the ratio of the gene expression between homologous gene pairs. N values indicate the number of dominant genes (fold change > 2) in each subgenome. The p-values indicate the significance of subgenome dominance by comparing the expression levels of homologous gene pairs from the two subgenomes based on Wilcoxon rank-sum test. c The heatmap of differentially expressed genes after Rickettsia and SFTSV infection. Each cell in the heatmap represents the fold-change of the genes with same Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation, with upregulated genes in red on the left and downregulated genes in blue on the right. d Differentially expressed genes in the Toll and Imd pathways after R. heilongjiangensis infection. The upregulated genes are marked in rose pink, and the downregulated genes are marked in light blue. The rounded rectangles highlight genes with significant transcriptional changes, and dashed lines indicate genes that are not annotated by KEGG.
We analyzed the expression differences of homologous gene pairs across the three subgenomes after infection with two pathogens. In both the uninfected control group and the groups infected with Rickettsia or SFTSV, the expression levels of homologous genes in subgenomes B and C were significantly higher than in subgenome A (Fig. 4b; Wilcoxon signed-rank test, all P < 0.05). This suggests that subgenomes B and C may play a more prominent role in response to pathogen infections. Following infection with R. heilongjiangensis and SFTSV, parthenogenetic H. longicornis exhibited 1048 (459 upregulated and 589 downregulated) and 1120 (387 upregulated and 733 downregulated) differentially expressed genes (DEGs), respectively. In contrast, the bisexual ticks had fewer DEGs, with 797 (359 upregulated and 438 downregulated) and 746 (323 upregulated and 423 downregulated) for the respective pathogens (Fig. 4c). These findings suggest that parthenogenetic ticks may have a more complex gene regulatory response to pathogen infections.
Subsequently, to explore the immunological evolutionary characteristics of parthenogenetic H. longicornis, we compared the differentially expressed genes in the Toll and Imd pathways between parthenogenetic and bisexual populations. After Rickettsia infection, the expression levels of extracellular cytokines Spatzle (SPZ), transmembrane cytokine receptor (Toll), interleukin-1 receptor-associated kinase 1 (Pelle), and dual oxidase (DUOX) were upregulated in both parthenogenetic and bisexual ticks (Fig. 4d), emphasizing the important roles of these genes in immune signaling and antibacterial defense28. Notably, this observation aligns with prior studies documenting the profound reduction of the IMD pathway in ticks29, which may explain the limited number of differentially expressed genes (DETs) detected in this pathway. However, our findings suggest that remnant components of the IMD pathway, such as Tube and Bendless, retain pathogen-responsive functionality: parthenogenetic ticks specifically upregulated interleukin-1 receptor-associated kinase 4 (Tube) and E2 ubiquitin-conjugating enzyme (Bendless) post-infection, which are linked to pathogen defense and intracellular signaling regulation30,31,32. At the same time, it specifically downregulated c-Jun N-terminal kinase (JNKK), which may have contributed to the increased rickettsia infection33,34. Following infection with the SFTSV, the parthenogenetic H. longicornis specifically downregulated X-linked inhibitor of apoptosis protein (XIAP) and transforming growth factor-β activated kinase 1 (TAK1) (Supplementary Fig. 17), which may contribute to immune evasion35,36. As parthenogenetic H. longicornis exhibited a higher Rickettsia load and faster replication of SFTSV, it is speculated that its complex immune regulatory mechanisms enhance pathogen tolerance rather than eliminating pathogens37, thereby increasing its capacity to serve as a pathogen vector38.
To further investigate potential microbial factors contributing to the differences in vector competence, we compared the composition and activity of Coxiella symbionts in parthenogenetic and bisexual H. longicornis. Metagenomic analysis based on whole-genome sequencing of five individuals from each group revealed that Coxiella burnetii was the predominant symbiont in both populations, with strain-level similarity (average ANI: 99.64%). Notably, Candidatus Coxiella mudrowiae was exclusively detected in several bisexual individuals but absent in the parthenogenetic group. Furthermore, mapping ovary transcriptomes from three feeding stages (early, partially engorged, and fully engorged) to the Coxiella genome revealed significantly higher transcript abundance in parthenogenetic ovaries, especially during the fully engorged stage (Student’s t test, P < 0.05) (Supplementary Fig. 18). These results suggest that Coxiella symbionts may play a more active role in parthenogenetic ticks during blood feeding, potentially contributing to their enhanced tolerance to pathogens and higher vector capacity.
Greater dispersal potential
We further investigated the role of parthenogenetic population structure in the global dissemination of H. longicornis. Parthenogenetic H. longicornis is morphologically similar to the bisexual ticks, leading to a lack of distinction in most public records. To differentiate between the two, we conducted next-generation sequencing on five parthenogenetic and five bisexual laboratory ticks. Utilizing the GATK4 pipeline (Supplementary Fig. 19), we identified 192 single nucleotide polymorphisms (SNPs) in the mitochondrial genome, demonstrating that parthenogenetic and bisexual ticks can be distinguished based on these mitochondrial SNPs. Subsequently, we analyzed 465 mitochondrial genomes from H. longicornis samples collected throughout China, including our published39 and unpublished data, and 87 full-length mitochondrial sequences obtained from GenBank (Supplementary Table 4). A total of 127 sequences were identified as belonging to parthenogenetic ticks, including 119 from China, three from the USA, two from Australia, and one each from New Zealand, Japan, and South Korea (Fig. 5a). Phylogenetic analysis revealed that the parthenogenetic population was divided into three clades: clade 1, primarily from southwestern China (Chongqing, Sichuan, Guizhou provinces); clade 2, from central and eastern China; and clade 3 from eastern China and other countries (Supplementary Fig. 20). This suggests that clade 3 is crucial in the global geographic expansion of H. longicornis. Bisexual H. longicornis sequences from China and Japan formed five clades, with the closest branches to the parthenogenetic ticks found in Japan and various regions of China. Further geographical genetic distance analysis revealed that the geographic distribution range of parthenogenetic population was 4.3 times wider than that of bisexual population within the same genetic distance range (Fig. 5b). The mean dispersal index for the parthenogenetic population was 1.8 × 10⁷, significantly higher than 9.5 × 10⁵ for the bisexual population (Student’s t test, P < 0.001), indicating a greater dispersal capacity of parthenogenetic H. longicornis.
a Phylogenetic tree of parthenogenetic and bisexual H. longicornis across the world. Maximum-likelihood method are used to build the tree using mitochondrial genomes of ticks collected in China and published mitochondrial genomes from GenBank. Parthenogenetic H. longicornis is indicated in red, and bisexual H. longicornis in blue. b Phylogeographic analysis of parthenogenetic and bisexual H. longicornis. Each point represents a pairwise comparison between two tick samples in terms of geographic distance (x axis) and genetic distance based on mitochondrial sequence (y axis). c Differences in genome-wide mutations between parthenogenetic and bisexual H. longicornis. Each point represents the average number of mutations per subgenome of parthenogenetic H. longicornis (n = 12) and bisexual H. longicornis (n = 12) deep sequencing samples. Wilcoxon rank-sum test is used to test the significance of the difference between two populations for each type of mutations. The line in the middle of each box plot represents the median of the dataset; the upper and lower edges of box plot indicate the third quartile and first quartile, respectively; and the line extending from the edge is 1.5 times the interquartile range. d Morphological and reproductive characteristics of parthenogenetic and bisexual H. longicornis. Wilcoxon rank-sum test is used to test the significant differences between the two populations. The p-values < 0.01 are summarized with three asterisks, and p-values < 0.0001 are summarized with four asterisks. e Potential global distribution map of parthenogenetic H. longicornis and bisexual H. longicornis. The probabilities of habitat suitability are shown for parthenogenetic H. longicornis (in red) and bisexual H. longicornis (in blue) predicted by ecological niche modeling.
It is generally believed that parthenogenetic species face challenges in overcoming deleterious mutation accumulation due to reduced genetic diversity, leading to decreased environmental adaptability40. This perspective appears contradictory to the observed widespread expansion of parthenogenetic H. longicornis. To investigate this, we compared deleterious (loss-of-function) mutations between field-collected parthenogenetic and bisexual ticks. The number of deleterious mutations per subgenome in parthenogenetic ticks was 49 on average, significantly lower than the mean of 905 in bisexual ticks (Wilcoxon rank-sum test, P < 0.0001), challenging the conventional understanding that deleterious mutation accumulation reduces the adaptability of parthenogenetic species17. Moreover, non-synonymous mutations in parthenogenetic ticks constituted only 5.6%, and synonymous mutations were 7.1% of those found in bisexual populations, both with significant differences (Wilcoxon rank-sum test, both P < 0.001, Fig. 5c). These findings suggest a lower mutation rate in the parthenogenetic population, which might be attributed to the lack of genetic recombination in parthenogenesis17. However, without the capacity for genetic recombination, these ticks face significant constraints in incorporating new beneficial mutations. An alternative explanation may involve concerted evolution within these genomes, which could promote sequence homogenization across subgenomes and facilitate error correction of deleterious mutations41. Additionally, the number of positively selected genes in the three subgenomes of the parthenogenetic population was 63, 47, and 57, respectively, all lower than the 314 identified in the bisexual population (Supplementary Fig. 21), further supporting the lower mutation rate in the parthenogenetic population. Collectively, both the reduced mutation rate and lower accumulation of deleterious mutations underpin the current adaptive advantage of parthenogenetic H. longicornis.
To assess the phenotypic impact of the long-term accumulation of deleterious mutations, we compared the morphological and reproductive characteristics between parthenogenetic and bisexual populations (Fig. 5d, Supplementary Fig. 22). Compared to bisexual population, parthenogenetic ticks exhibited significantly larger dorsal shield (4.16 ± 0.16 vs 3.68 ± 0.25 mm²), genital pores (357.21 ± 25.49 vs 321.21 ± 27.52 μm), and eggs (76.60 ± 3.30 vs 72.35 ± 3.56 μg). Additionally, because mating is not required, the blood-feeding and egg incubation periods of parthenogenetic H. longicornis were significantly shorter than those of bisexual population. Polyploidization brings unique physiological changes to parthenogenetic populations and is crucial for cell growth and yolk biosynthesis. It may lead to larger cell volumes and faster growth rates42. The reduced hatching rate of parthenogenetic H. longicornis eggs (80.85 ± 7.72% vs 89.01 ± 7.72%) might reflect the instability of the triploidy on embryonic cell development. These phenotype results further indicate that parthenogenetic H. longicornis might have successfully overcome the reproductive and developmental issues imposed by the triploid genome.
To further support this view, we conducted ecological model prediction for the two populations (Fig. 5e). The model predicted that the suitable habitats of bisexual H. longicornis were confined within East Asia (Supplementary Figs. 23, 24). In contrast, ecological niche modeling predicted that parthenogenetic H. longicornis could inhabit a broad range of suitable habitats across multiple continents (Supplementary Figs. 25, 26), including East Asia (China, North Korea, South Korea, Japan), Oceania (Australia, New Zealand), and North America (e.g., the United States). Notably, South American countries such as Uruguay, Chile, and Brazil were also identified as ecologically suitable, despite limited or no current records. These findings highlight the potential for future spread into new regions, particularly in South America, if accidental introductions occur.
Discussion
We assembled and annotated a high-quality genome of parthenogenetic H. longicornis, achieving superior genome assembly continuity and completeness using ultra-depth HiFi long-read sequencing data. Subsequently, a chromosome-level assembly of the initial haplotypes was completed through Hi-C data, resulting in three high-quality subgenomes. Additionally, Iso-Seq was performed on parthenogenetic H. longicornis, enabling a comprehensive annotation of alternative splicing events. Most polyploid genome assemblies come from allopolyploid species, such as hexaploid gibel carp (AAABBB) and a tetraploid frog (LLSS)43,44, which typically include a single representative genome for multiple haplotypes of the same origin but rarely differentiate between them. In this study, we differentiated the three haplotypes of the same origin and obtained a fully-phased genome assembly, providing critical data investigate the complex genetic and evolutionary mechanisms underlying parthenogenesis in H. longicornis.
In the allopolyploid genomes, subgenomes can experience varying degree of selection pressures, biased gene retentions, and biased gene expression (i.e. subgenome dominance)11,18,45. Subgenome dominance has not been frequently reported in young autopolyploid species, likely due to lack of phased polyploid genome. By integrating the phased triploid genome and transcriptome data, it was observed that homologous genes in subgenome A of parthenogenetic H. longicornis exhibit higher expression levels during key stages of embryonic development, while subgenomes B and C show higher expression levels following pathogen infection, demonstrating functional divergence among the subgenomes in H. longicornis. Further analysis of gene expression data revealed the unique expression of multiple genes involved in the cell cycle and oocyte meiosis pathways within the ovary of blood-feeding parthenogenetic H. longicornis. These genes may play crucial roles in overcoming constraints on embryonic development and promoting genomic stability in parthenogenetic populations25. Additionally, parthenogenetic H. longicornis exhibited higher pathogen loads and more sophisticated immune manipulation, suggesting a heightened pathogen tolerance and an enhanced capacity to serve as a vector.
To assess the dispersal potential of parthenogenetic H. longicornis, mitochondrial genome data were used to distinguish parthenogenetic from bisexual H. longicornis populations around the world, which had not been systematically categorized in previous distribution studies5. It was found that bisexual populations are mainly confined to East Asia, while parthenogenetic populations have spread to the Americas and Australia, occupying a wider geographic range. A major challenge faced by parthenogenetic species is the accumulation of deleterious mutations in the absence of sexual recombination, which could lead to the extinction of parthenogenetic population40. However, we found that the genome of parthenogenetic H. longicornis exhibited lower levels of deleterious mutations compared to bisexual populations and did not show signs of mutational degradation, as seen in clonally reproducing fish46. Furthermore, the comparison of the morphological and reproductive characteristics of the two populations further demonstrated that parthenogenetic populations have successfully overcome the reproductive degeneration pressures caused by asexuality and triploid instability17.
In conclusion, our study suggests that the evolutionary origin of parthenogenetic H. longicornis likely stems from an autopolyploid event. The triploid genome and parthenogenetic mode of reproduction have provided several advantages to H. longicornis, including increased reproductive efficiency, heightened pathogen tolerance, enhanced vector capacity, and greater dispersal potential. Future studies should focus on the global surveilence and transmitted pathogens of parthenogenetic H. longicornis to better understand the public health implications of this unique tick population.
Methods
Ticks and pathogens
The parthenogenetic and bisexual Haemaphysalis longicornis were collected by flag-dragging and 50 ml ventilated centrifuge tube from Cangxi County (31◦ 44’ 35” N, 105◦ 49’ 04” E), Sichuan Province, and Yu County (40◦ 03’ 03” N, 115◦ 23’ 15” E), Hebei Province, China. The collection of wild ticks was carried out with the permission and assistance of the local forestry protection agency. The parthenogenetic population have undergone more than ten generations of in-bred laboratory culture, were maintained in an incubator at 26 ± 1°C, 85 ± 5% relative humidity, with a 12:12 h light-dark cycle during the non-blood-feeding phase. The strains of Rickettsia heilongjiangensis and severe fever with thrombocytopenia syndrome virus (SFTSV) used in this study were isolated from field-collected H. longicornis ticks and previously characterized in our laboratory. Ticks were placed on experimental rabbits’ ears during blood-feeding. The rabbits used for tick feeding were New Zealand White Closed Colony, female, 8 weeks old (about 2.2 kg) at the start of the experiment. They were individually housed in 60 × 50 × 45 cm stainless steel cages (wood shavings changed every 2 days), with the housing room maintained at 22–24 °C, 50–60% relative humidity, and a 12 h light/12 h dark cycle. All experimental procedures adhered to the protocols approved by the Institutional Animal Care and Use Committee (IACUC) of the Beijing Institute of Microbiology and Epidemiology. The animal ethics approval number is IACUC-DWZX-027-20. We have complied with all relevant ethical regulations for animal use.
Genome ploidy analysis
The differentiation between parthenogenetic and bisexual H. longicornis was achieved through karyotype analysis using Sysmex Partec CyStain DNA 1 Step Kit (Sysmex Partec, Germany)9. The basis capitulum of adult ticks from parthenogenetic and bisexual populations were carefully severed, and the visceral tissues were dissected using fine forceps. The tissues were placed in a plastic dish containing 500 μl of CyStain UV Ploidy reagent and finely minced using a sharp blade to fully release the nuclei, and soaked for 1 min. The resulting suspension was filtered through a 50 μm CellTrics filter into a sample tube. Subsequently, 2 ml of CyStain UV Ploidy reagent was added and incubated in the dark for 2 min. The nuclei suspension was then analyzed using the CyFlow Space Flow Cytometer (Sysmex Partec, Germany) and the corresponding FloMax software.
Genome size estimation
Prior to de novo assembling, the genome size of the parthenogenetic H. longicornis was estimated using an Illumina short-read DNA sequencing dataset. Jellyfish (v2.2.10)47 was employed to determine the k-mer frequency, and GenomeScope 2.048 was utilized to estimate the genome size and heterozygosity rate.
Genome sequencing
A total of 17.730 μg of genomic DNA was extracted from the offspring of adult ticks (n = 90) from a single mother for Pacific Biosciences CCS library construction. The DNA quality was verified using agarose gel electrophoresis, confirming the high integrity of the DNA molecules49. For Pacific Biosciences sequencing, the genomic DNA was sheared into ~15 Kb fragments to construct a long-read library according to the manufacturer’s instructions, and then the library was sequenced on the Pacific Biosciences Sequel II system utilizing an 8 M SMRT Cell with the Sequel II Sequencing Kit (Pacific Biosciences, USA).
Transcriptome sequencing
Transcriptome samples were categorized into three groups: parthenogenetic female (n = 5), bisexual female (n = 5), and bisexual male (n = 5). Each group had one Iso-Seq sample and three RNA-seq samples sequenced to optimize sequence accuracy and completeness. Total RNA was extracted and purified from 5 ticks whole body for each sample using miRNeasy Mini Kit (QIAGEN, Germany). Transcripts Per Million (TPM) was employed as the normalization unit for the quantitative analysis of transcriptomic data. TPM values for samples in each group have been organized into a table (Supplementary Data 1), encompassing gene expression data from all experimental groups and their biological replicate samples (n = 3). Ovarian transcriptome data were derived from our previously published study22. Ovaries were dissected from parthenogenetic and bisexually populations at three distinct blood-feeding stages: early (48 h post-feeding), partially (96 h), and engorged. For each experimental group, 30 ovarian tissues were collected, homogenized, and pooled into three RNA-seq samples, which were then subjected to paired-end RNA sequencing using the Illumina Hiseq platform (Illumina, USA).
Hi-C sequencing
Chromatin conformation capture data were obtained from adult ticks (n = 15) issued from the same mother as those used in genome sequencing. The ticks were dissected into several sections using scissors, and 2% formaldehyde was used to crosslink proteins with DNA and DNA with DNA fragments within the cells. The DNA was then digested with DpnII enzyme (New England Biolabs, USA), creating sticky ends on either side of the crosslinked fragments. A suitable labeling system was prepared to biotin-label the DNA ends. The blunt ends were ligated using T4 DNA ligase. The DNA was de-crosslinked and purified, followed by fragmentation of 1–4 µg samples using the Covaris S220 (Covaris, USA), with 300–700 bp fragments recovered. Streptavidin magnetic beads were used to capture the DNA fragments that retained interaction signals, followed by library construction. The concentration and insert size of the library were assessed using the Qubit 3.0 and GX platforms. Sequencing was performed with the Illumina platform using paired-end (PE) 150 base pair reads.
De novo assembly
The genome of H. longicornis was assembled using hifiasm (v0.19) with the obtained HiFi reads and parameters “l = 0, n = 4”. Genome assembly completeness was evaluated using BUSCO (v5.4.2)50 in protein mode (-m prot) with the arthropoda_odb10 dataset. The accuracy of the assembly was verified by mapping all Illumina paired-end reads to the assembled genome using BWA (v0.7.17)51 and HiFi reads using Minimap2 (v2.26)52. The mapping rate and genome coverage of sequencing reads were assessed using samtools (v1.15).
Genome annotation
Protein-coding genes were annotated using three approaches: ab initio prediction, homology prediction, and transcriptome prediction. Ab initio prediction was performed with Augustus (v3.1.0)53 and SNAP (2006-07-28)54. Homology-based prediction used GeMoMa (v1.7)55, referencing genes from H. longicornis (GCA_013339765), Ixodes scapularis (GCF_016920785), Rhipicephalus microplus (GCF_013339725), Dermacentor silvarum (GCF_013339745), and Drosophila melanogaster (GCF_000001215). For transcriptome prediction, three methods were used to assemble transcripts. The first method utilized Hisat (v2.1.0)56 and Stringtie (v2.1.4)57 for transcript assembly, followed by gene prediction with GeneMarkS-T (v5.1)57. The second method assembled transcripts using Trinity (v2.11)58 and then performed gene prediction with PASA (v2.4.1)59. The third method aligned Iso-Seq transcripts with gmap (2020-06-30) and predicted genes with PASA (v2.4.1)59. All the prediction results were integrated using EVM (v1.1.1)60 and refined with PASA (v2.4.1)59. Gene sequences were annotated against databases including NR, eggNOG, KEGG, SWISS-PROT, and Pfam61,62,63,64.
GO and KEGG enrichment analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the topGO and clusterProfiler65 packages in R. GO term enrichment was conducted using Fisher’s exact test in the topGO package, while KEGG pathway enrichment was performed with the enricher function in clusterProfiler. For both analyses, p values were adjusted for multiple testing using the Benjamini–Hochberg method implemented via the p.adjust function in R (method = “BH”). Terms with adjusted p values (FDR) < 0.05 were considered significantly enriched.
Chromosome assembly using Hi-C data
Following Hi-C library sequencing, paired-end reads were mapped to the draft assembly using BWA (v0.7.17). Valid interaction pairs were identified using HiC-Pro (v2.8.1)66, and the LACHESIS pipeline generated chromosome-level scaffolds. For the synteny analysis, synteny blocks were identified using MCScanX (version not specified)67 and JCVI (v0.9.13)68 with a cscore threshold of >0.5.
Repeat annotation
Transposable elements (TEs) were initially predicted using RepeatModeler2 (v2.0.1)69, which incorporated RECON (v1.0.8)70 and RepeatScout (v1.0.6)71. The predicted results were classified with RepeatClassifier using the known Dfam (v3.5) database. Subsequently, de novo prediction of LTRs was specifically performed with LTR_retriever (v2.9.0)72, relying on prediction results from LTRharvest (v1.5.10)73 and LTR_FINDER (v1.07)74. The de novo prediction results were then merged and deduplicated with known databases to create a species-specific repeat sequence database. Finally, RepeatMasker (v4.1.2)75 was employed to predict transposable elements in the genome based on the constructed repeat sequence database. Tandem repeat sequences were primarily predicted using the MIcroSAtellite identification tool (v2.1)76 and Tandem Repeat Finder (version 4.09, parameters: 2 7 7 80 10 50 500 -d -h)77.
Non-coding RNA annotation
Five types of non-coding RNA were annotated. tRNAs were identified using tRNAscan-SE (v1.3.1)78, rRNAs with barrnap (v0.9)79. miRNAs, snoRNAs, and snRNAs were predicted using the Rfam (v14.5)80 database via Infernal (v1.1)81.
Pathogen infection
Both R. heilongjiangensis and SFTSV were cultured in Vero81 cells using DMEM medium supplemented with 5% fetal bovine serum (FBS), maintained at 37 °C in a 5% CO₂ incubator. Pathogen injection was performed on both parthenogenetic and bisexual populations using the Eppendorf FemtoJet 4× microinjector (Eppendorf, Germany). Each tick was injected with 0.5 μL of purified pathogen culture via the lateral margin ventral of spiracular plate, with three adult ticks per group. After injection, ticks were not offered blood meals and were maintained in a thermostatic incubator (26 ± 1°C, relative humidity 75%, L: D = 12 h: 12 h). Pathogen proliferation in ticks was monitored at 5d, 10d, 15d, 20d and 25d post-infection using the TB Green® Premix Ex Taq™ kit (Takara, Japan) to extract RNA from whole ticks. The loads of R. heilongjiangensis and SFTSV were detected using quantitative real-time PCR with specific primers (R. heilongjiangensis, targeting Sca1 gene: Forward 5’-GTTTGTGGATGCGTGGTA-3’, Reverse 5’-AACCCGATAGTAGCAC-3’; SFTSV, targeting L segment: Forward 5’-ACATTCCATGCCATCTCAGG-3’, Reverse 5’-GTCCTCTTCCATCTCCACRATC-3’). Upon establishing the pathogen proliferation curves within H. longicornis, RNA-seq sequencing was performed on both infected parthenogenetic and bisexual populations at 20 d post-infection. Uninfected adult ticks served as the control group.
Differential expression analysis
We collected samples of parthenogenetic and bisexual H. longicornis both before and after infection for RNA extraction and transcriptome sequencing. The RNA-Seq reads were aligned to the reference genome using HISAT2 (v2.1.0)56. We used featureCounts82 to count the number of mapped reads on each gene and transcript, and estimated gene expression levels using transcripts per kilobase million mapped reads (TPM). Statistically differentially expressed genes were identified using edgeR (v3.28.1)83 with a false discovery rate (FDR) ≤ 0.05 (Benjamini–Hochberg method) and |log2(fold change)| ≥ 1. Enriched KEGG pathways and GO terms were identified using the clusterProfiler (v3.6.1)84 package (FDR < 0.05).
Iso-Seq analysis pipeline
The adult female samples of parthenogenetic was used for RNA extraction and library construction following standard protocol of Iso-Seq from Pacific Bioscience. Each library was sequenced on the PacBio Sequel-II Systems (Pacific Bioscience, CA). The isoseq3 pipeline (v3.8.2) (https://github.com/PacificBiosciences/IsoSeq) characterized the full-length transcripts through three major steps: generating CCS reads, classifying full-length (FL) reads, and clustering full-length non-chimeric (FLnc) reads. First, CCS reads were generated using CCS (v6.4.0) from the subreads BAM files with a minimum quality of 0.9 (–min-rq 0.9) and minimum number of FL subreads (n = 3). Second, Lima (v2.7.1) was used to remove the cDNA primers. The polyA tail and artificial concatemers were removed to produce FLnc reads. Third, FLnc reads were clustered by Isoseq3 to generate high quality FL consensus sequences. The high-quality FL reads were then aligned to the reference genome using pbmm2, and redundant transcripts were collapsed based on exonic structures. Finally, pigeon (v1.0.0) was used to classify transcripts against annotations and filter out potential artifacts.
Gene family and phylogenetic analysis
To elucidate the evolutionary history of H. longicornis, a maximum-likelihood phylogenetic tree was constructed using protein sequences from the three subgenomes of parthenogenetic H. longicornis, bisexual H. longicornis (China National Center for Bioinformation: PRJCA040748) and seven other tick species14,15,85 and two mites (GCF_002443255.1 and GCF_000255335.2). OrthoFinder (v2.5.4)86 clustered protein sequences into orthogroups, alignments were built using MUSCLE (v3.6)87, and a maximum-likelihood tree was constructed with IQ-TREE (v2.2.0-beta)88. Divergence time of each node was estimated using MCMCTREE (v4.9) in PAML (v4.10.6)89 using an independent-rates model (clock = 2). The calibration points obtained from the TimeTree database (http://timetree.org/)90 included: (i) 336.0–453.4 million years ago for Parasitiformes; (i) the divergence of Ixodidae (62.1–316.5 Mya); (ii) and 55.9–134.0 million years ago for Metastriata; (ii) the split between Rhipicephalus sanguineus and Rhipicephalus microplus (37.0–52.0 Mya), and (iii) the split between R. sanguineus and Dermacentor silvarum (46.6–118.3 Mya). Additionally, we applied a soft maximum bound of 33.0 Mya for the divergence between BHL and PHL, based on the estimated divergence time between Haemaphysalis longicornis and its closest relative, Haemaphysalis flava. As there is no direct divergence estimate between Ixodes scapularis and Ixodes persulcatus in the TimeTree database—and given that I. scapularis is more distantly related to Ixodes persulcatus than I. ricinus—we set a minimum bound of 47 Mya for the divergence between I. persulcatus and I. scapularis, based on the estimated divergence time between I. persulcatus and I. ricinus.
The CAFÉ (v5.0)91 was used to analyze gene family expansion and contraction by comparing family sizes between parent and child nodes in the phylogenetic tree, using a p-value threshold of 0.05 to identify significant changes. For gene presence/absence analyses based on KO annotation, particularly for functionally important gene families, we adopted a conservative strategy. In addition to the HLB genome from Jia et al.14, we incorporated a second independently annotated bisexual genome65, which has a BUSCO completeness of 87.1%. A gene family was considered truly absent only if it was missing from both assemblies.
Positive selection analysis
For the genes in each orthogroup identified in the “Gene Family and Phylogenetic Analysis” section, protein sequences were aligned using MAFFT (v7.505)92 and gene trees were constructed with FastTree (v2.1.11)93. Paralogous genes were pruned from these alignments using PhyloPYPruner (v1.2.4)94, retaining only one-to-one orthologs present in at least seven species or subgenomes. The corresponding CDS sequences were then aligned using Prank (v.170427)95 based on the protein sequence alignments. CODEML from the PAML package90 was used to detect signals of positive selection under branch-site models along the predefined branch. A likelihood ratio test was performed to assess the significance of the branch-site model against the null model. The p-values were calculated using Chi-square statistics, and an FDR < 0.05 (BH method) was applied to eliminate potential false positives. Pairwise Ka/Ks ratio between the orthologous sequence pairs from each parthenogenetic subgenome and bisexual H. longicornis genome were estimated using “wgd ksd –one_v_one” command (v1.1)96.
Mitochondrial phylogenetic analysis
Based on the second-generation sequencing data from laboratory-reared parthenogenetic and bisexual populations, the full-length mitogenomes were de novo assembled using the “all” module from MitoZ (v2.3)97. We then identified the variants in the full-length mitochondrial genomes of both populations using Geneious Prime (version 2023.2)98, confirming that the mitochondrial genome can differentiate between the two populations. Next, we analyzed 465 mitochondrial genomes from our sequencing database of H. longicornis collected in the field for SNP calling. Variants in the mitochondrial genome were called using the GATK4 (v4.1.2.0) joint genotype pipeline99. Variants were filtered to include those with a read depth of ≥3 and a quality score of ≥30. Additionally, individuals with a genotype rate of >0.9 were included to generate the mitochondrial consensus sequence for each sample using bcftools (v1.9)100. We then aligned the sequences using MAFFT (v7.505)92 and constructed the phylogenetic tree with IQ-TREE (v2.2.0-beta)88, utilizing the best-fit substitution model and ultrafast bootstrap with 1000 replicates.
Geographical genetic distance analysis
In this study, nucleotide sequences were processed to remove redundancy using the CD-HIT-EST program with a threshold of 0.99100. The genetic distances for each sequence were calculated using MEGA11101. Geographic distances were determined based on latitude and longitude using the distm function from the geosphere R package. We utilized the dispersal index (I) to indicate dispersal velocity and performed an independent sample t-test to compare the dispersal rates of bisexual and parthenogenetic ticks. The dispersal index and t-test calculations were performed using custom scripts written in R version 4.4.
Variant calling and mutation impact prediction
The whole genome Illumina reads of the parthenogenetic H. longicornis and bisexual H. longicornis individuals were mapped to their respective genomes using BWA (Version 0.7.17-r1188)51. Duplicate reads were marked using Picard tools (v2.18.23, http://broadinstitute.github.io/picard/). SNPs were called using FreeBayes (v1.3.6)102 and filtered based on the following thresholds: 1) minimum read coverage ≥7 for parthenogenetic H. longicornis and ≥10 for bisexual H. longicornis; 2) for heterozygous variant calls, the allele frequency of alternative alleles ≥2/7; 3) for homozygous variant calls, the allele frequency of alternative alleles ≥6/7; 4) variant quality score ≥20, genotype quality ≥20, mapping quality of reads supporting the variant ≥20, and fraction of properly paired reads ≥0.9. The potential impact of called variants, including loss of function, synonymous, and non-synonymous mutations, were evaluated using SnpEff (v5.2a)103 with default parameters.
Phenotype analysis
Physiological and reproductive characteristics were recorded for both parthenogenetic and bisexual females104. Measurements included dorsal area (n = 30), genital pore length (n = 30), unfed weight (n = 20), engorged quantity (n = 20), egg weight (n = 20), egg mass (n = 20), blood-feeding weight (n = 20), egg perimeter (n = 60), blood-feeding stage (n = 20), preoviposition stage (n = 20), oviposition stage (n = 20), incubation stage (n = 20), hatchability (n = 20), blood-sucking index (n = 20) and reproductive index (n = 20). Statistical comparisons between the two populations were performed using Wilcoxon rank-sum test.
Modeling the distribution of parthenogenetic H. longicornis
To model the distribution of parthenogenetic and bisexual H. longicornis, the exact locations (longitude and latitude coordinates) of the collection points for these two populations were extracted from previous reports and our field survey data. In our field survey data, the identification of parthenogenetic H. longicornis or bisexual H. longicornis was based on the assembled mitochondrial sequences. All reported occurrences of H. longicornis in Australia, New Zealand, and the USA were treated as parthenogenetic H. longicornis since parthenogenetic H. longicornis are the dominant population in these areas105,106. The reported occurrences of male H. longicornis were treated as bisexual H. longicornis.
The environmental and meteorological factors were collected to predict the potential distributions of H. longicornis. The temperature and precipitation data from the WorldClim database (https://worldclim.org). The elevation, vegetation (Percent Tree Cover) and the global land cover (GLCNMO) data were downloaded from the Global Map data archives website (https://globalmaps.github.io/). Slope aspect and slope degree data was extracted based on elevation data using ArcMap v10.7. Highly correlated factors and duplicated distribution points were removed using ENMTools (v1.4.4)107. Then an ecological niche model was applied to predict potential distribution of each population, and Maxent (v3.4.4)108 was applied to train the models. Parameter calibration, model evaluation, and selection of candidate models were accomplished using kuenm package in R (version 3.6.3)109.
Statistics and reproducibility
All statistical analyses were performed using R software (v4.4.0). The statistical methods and corresponding p-values for each test are provided in the main text or figure legends. Statistical significance is indicated as follows: p < 0.05 (*), p < 0.01 (**), p < 0.001 (***), and p < 0.0001 (****). Data are presented as the mean ± standard deviation (SD). The number of biological replicates (n) for each experiment is specified in the figure legends.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Data availability
The genome assemblies and annotations generated in this study are available at GenBank with project ID: PRJNA1133666. We have deposited the raw sequencing data in the China National Center for Bioinformation (CNCB) under accession number CRA031217. This repository includes all raw HiFi, Iso-seq, and Hi-C data as required. Additionally, the numerical source data for all figures have been deposited in Figshare (https://doi.org/10.6084/m9.figshare.30373816). The RNA-seq data of uninfected ticks and ticks infected with two pathogens at GenBank with accession numbers SRR30735266-SRR30735286. The metagenomic sequencing data of tick whole body are available at GenBank with accession numbers SRR30788678-SRR30788687.
References
Hoogstraal, H. Review of Haemaphysalis (Kaiseriana) longicornis Neumann of Australia, New Zealand, New Caledonia, Fiji, Japan, Korea, and northeastern China and USSR, and its parthenogenetic and bisexual populations (Ixodoidea, Ixodidae). J. Parasitol. 54, 1197–1213 (1968).
Google Scholar
Rainey, T., Occi, J. L., Robbins, R. G. & Egizi, A. Discovery of Haemaphysalis longicornis (Ixodida: Ixodidae) parasitizing a sheep in New Jersey, United States. J. Med. Entomol. 55, 757–759 (2018).
Google Scholar
United States Department of Agriculture. National Haemaphysalis longicornis (Asian longhorned tick) situation report. USDA Animal Plant Health Inspection Service (United States Department of Agriculture, 2023).
Zhao, L. et al. Distribution of Haemaphysalis longicornis and associated pathogens: analysis of pooled data from a China field survey and global published data. Lancet Planet. Health 4, e320–e329 (2020).
Google Scholar
Zhuang, L. et al. Transmission of severe fever with thrombocytopenia syndrome virus by Haemaphysalis longicornis Ticks, China. Emerg. Infect. Dis. 24, 868–871 (2018).
Google Scholar
Fang, L. Q. et al. Emerging tick-borne infections in mainland China: an increasing public health threat. Lancet Infect. Dis. 15, 1467–1479 (2015).
Google Scholar
Butlin, R. Evolution of sex: the costs and benefits of sex: New insights from old asexual lineages. Nat. Rev. Genet. 3, 311–317 (2002).
Google Scholar
Zhang, X. et al. Rapid spread of severe fever with thrombocytopenia syndrome virus by parthenogenetic Asian longhorned ticks. Emerg. Infect. Dis. 28, 363–372 (2022).
Google Scholar
Chen, X. et al. Multiple lines of evidence on the genetic relatedness of the parthenogenetic and bisexual Haemaphysalis longicornis (Acari: Ixodidae). Infect. Genet. Evol. 21, 308–314 (2014).
Google Scholar
Comai, L. The advantages and disadvantages of being polyploid. Nat. Rev. Genet. 6, 836–846 (2005).
Google Scholar
Li, J. T. et al. Parallel subgenome structure and divergent expression evolution of allo-tetraploid common carp and goldfish. Nat. Genet. 53, 1493–1503 (2021).
Google Scholar
Anatskaya, O. V. & Vinogradov, A. E. Polyploidy and Myc proto-oncogenes promote stress adaptation via epigenetic plasticity and gene regulatory network rewiring. Int. J. Mol. Sci. 23, 9691 (2022).
Google Scholar
David, K. T. Global gradients in the distribution of animal polyploids. Proc. Natl. Acad. Sci. USA 119, e2214070119 (2022).
Google Scholar
Jia, N. et al. Large-scale comparative analyses of tick genomes elucidate their genetic diversity and vector capacities. Cell 182, 1328–1340.e13 (2020).
Google Scholar
De, S. et al. A high-quality Ixodes scapularis genome advances tick science. Nat. Genet. 55, 301–311 (2023).
Google Scholar
Lv, Z., Nyarko, C. A., Ramtekey, V., Behn, H. & Mason, A. S. Defining autopolyploidy: cytology, genetics, and taxonomy. Am. J. Bot. 111, e16292 (2024).
Google Scholar
Jaron, K. S. et al. Genomic features of parthenogenetic animals. J. Hered. 112, 19–33 (2021).
Google Scholar
Cheng, F. et al. Gene retention, fractionation and subgenome differences in polyploid plants. Nat. Plants 4, 258–268 (2018).
Google Scholar
Jaquiéry, J. et al. Disentangling the causes for Faster-X evolution in aphids. Genome Biol. Evol. 10, 507–520 (2018).
Google Scholar
Li, Z. et al. Phylloxera and aphids show distinct features of genome evolution despite similar reproductive modes. Mol. Biol. Evol. 40, msad271 (2023).
Google Scholar
Van de Peer, Y., Mizrachi, E. & Marchal, K. The evolutionary significance of polyploidy. Nat. Rev. Genet. 18, 411–424 (2017).
Google Scholar
Wang, T. H. et al. The ovarian development genes of bisexual and parthenogenetic Haemaphysalis longicornis evaluated by transcriptomics and proteomics. Front. Vet. Sci. 8, 783404 (2021).
Google Scholar
Schnable, J. C., Springer, N. M. & Freeling, M. Differentiation of the maize subgenomes by genome dominance and both ancient and ongoing gene loss. Proc. Natl. Acad. Sci. USA 108, 4069–4074 (2011).
Google Scholar
Ramírez-González, R. H. et al. The transcriptional landscape of polyploid wheat. Science 361, eaar6089 (2018).
Google Scholar
Vardy, L., Pesin, J. A. & Orr-Weaver, T. L. Regulation of Cyclin A protein in meiosis and early embryogenesis. Proc. Natl. Acad. Sci. USA 106, 1838–1843 (2009).
Google Scholar
Chen, J. et al. Comparative proteomic analysis provides new insights into the molecular basis of thermal-induced parthenogenesis in silkworm (Bombyx mori). Insects 14, 134 (2023).
Google Scholar
Snyman, M. & Xu, S. Transcriptomics and the origin of obligate parthenogenesis. Heredity 131, 119–129 (2023).
Google Scholar
Shaw, D. K. et al. Vector immunity and evolutionary ecology: the harmonious dissonance. Trends Immunol. 39, 862–873 (2018).
Google Scholar
Jalovecka, M. et al. Activation of the tick Toll pathway to control infection of Ixodes ricinus by the apicomplexan parasite Babesia microti. PLoS Pathog. 20, e1012743 (2024).
Google Scholar
Shaw, D. K. et al. Infection-derived lipids elicit an immune deficiency circuit in arthropods. Nat. Commun. 8, 14401 (2017).
Google Scholar
Rosa, R. D. et al. Exploring the immune signalling pathway-related genes of the cattle tick Rhipicephalus microplus: from molecular characterization to transcriptional profile upon microbial challenge. Dev. Comp. Immunol. 59, 1–14 (2016).
Google Scholar
Fogaça, A. C. et al. Tick immune system: what is known, the interconnections, the gaps, and the challenges. Front. Immunol. 12, 628054 (2021).
Google Scholar
Rana, V. S., Kitsou, C., Dumler, J. S. & Pal, U. Immune evasion strategies of major tick-transmitted bacterial pathogens. Trends Microbiol 31, 62–75 (2023).
Google Scholar
Garver, L. S., de Almeida Oliveira, G. & Barillas-Mury, C. The JNK pathway is a key mediator of Anopheles gambiae antiplasmodial immunity. PLoS Pathog. 9, e1003622 (2013).
Google Scholar
Severo, M. S. et al. The E3 ubiquitin ligase XIAP restricts Anaplasma phagocytophilum colonization of Ixodes scapularis ticks. J. Infect. Dis. 208, 1830–1840 (2013).
Google Scholar
Paquette, N. et al. Serine/threonine acetylation of TGFβ-activated kinase (TAK1) by Yersinia pestis YopJ inhibits innate immune signaling. Proc. Natl. Acad. Sci. USA 109, 12710–12715 (2012).
Google Scholar
Schneider, D. S. & Ayres, J. S. Two ways to survive infection: what resistance and tolerance can teach us about treating infectious diseases. Nat. Rev. Immunol. 8, 889–895 (2008).
Google Scholar
Dharmarajan, G., Walker, K. D. & Lehmann, T. Variation in tolerance to parasites affects vectorial capacity of natural Asian tiger mosquito populations. Curr. Biol. 29, 3946–3952.e5 (2019).
Google Scholar
Ye, R. Z. et al. Virome diversity shaped by genetic evolution and ecological landscape of Haemaphysalis longicornis. Microbiome 12, 35 (2024).
Google Scholar
Kearney, M. R. et al. Parthenogenesis without costs in a grasshopper with hybrid origins. Science 376, 1110–1114 (2022).
Google Scholar
Jiang, X., Song, Q., Ye, W. & Chen, Z. J. Concerted genomic and epigenomic changes accompany stabilization of Arabidopsis allopolyploids. Nat. Ecol. Evol. 5, 1382–1393 (2021).
Google Scholar
Chen, Q. et al. Every gain comes with loss: ecological and physiological shifts associated with polyploidization in a pygmy frog. Mol. Biol. Evol. 42, msaf037 (2025).
Google Scholar
Wang, Y. et al. Comparative genome anatomy reveals evolutionary insights into a unique amphitriploid fish. Nat. Ecol. Evol. 6, 1354–1366 (2022).
Google Scholar
Session, A. M. et al. Genome evolution in the allotetraploid frog Xenopus laevis. Nature 538, 336–343 (2016).
Google Scholar
Renny-Byfield, S. et al. Ancient gene duplicates in Gossypium (cotton) exhibit near-complete expression divergence. Genome Biol. Evol. 6, 559–571 (2014).
Google Scholar
Kočí, J. et al. No evidence for accumulation of deleterious mutations and fitness degradation in clonal fish hybrids: Abandoning sex without regrets. Mol. Ecol. 29, 3038–3055 (2020).
Google Scholar
Marçais, G. & Kingsford, C. A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics 27, 764–770 (2011).
Google Scholar
Ranallo-Benavidez, T. R., Jaron, K. S. & Schatz, M. C. GenomeScope 2.0 and Smudgeplot for reference-free profiling of polyploid genomes. Nat. Commun. 11, 1432 (2020).
Google Scholar
Abu Almakarem, A. S. et al. Extraction of DNA from plant and fungus tissues in situ. BMC Res. Notes 5, 59 (2012).
Google Scholar
Burton, J. N. et al. Chromosome-scale scaffolding of de novo genome assemblies based on chromatin interactions. Nat. Biotechnol. 31, 1119–1125 (2013).
Google Scholar
Li, H. et al. The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079 (2009).
Google Scholar
Li, H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34, 3094–3100 (2018).
Google Scholar
Stanke, M., Diekhans, M., Baertsch, R. & Haussler, D. Using native and syntenically mapped cDNA alignments to improve de novo gene finding. Bioinformatics 24, 637–644 (2008).
Google Scholar
Korf, I. Gene finding in novel genomes. BMC Bioinforma. 5, 59 (2004).
Google Scholar
Keilwagen, J., Wenk, M., Erickson, J. L., Schattat, M. H. & Grosse, I. Using intron position conservation for homology-based gene prediction. Nucleic Acids Res. 44, e89 (2016).
Google Scholar
Kim, D., Langmead, B. & Salzberg, S. L. HISAT: a fast spliced aligner with low memory requirements. Nat. Methods 12, 357–360 (2015).
Google Scholar
Pertea, M. et al. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295 (2015).
Google Scholar
Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z. & Amit, I. Trinity: reconstructing a full-length transcriptome without a genome from RNA-seq data. Nat. Biotechnol. 29, 644 (2013).
Google Scholar
Haas, B. J. et al. Improving the Arabidopsis genome annotation using maximal transcript alignment assemblies. Nucleic Acids Res. 31, 5654–5666 (2003).
Google Scholar
Haas, B. J. et al. Automated eukaryotic gene structure annotation using EVidenceModeler and the Program to Assemble Spliced Alignments. Genome Biol. 9, R7 (2008).
Google Scholar
Huerta-Cepas, J. et al. eggNOG 5.0: a hierarchical, functionally and phylogenetically annotated orthology resource based on 5090 organisms and 2502 viruses. Nucleic Acids Res. 46, D308–D314 (2018).
Kanehisa, M., Sato, Y., Kawashima, M., Furumichi, M. & Tanabe, M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 44, D457–D462 (2016).
Google Scholar
Boeckmann, B. et al. The SWISS-PROT protein knowledgebase and its supplement TrEMBL in 2003. Nucleic Acids Res. 31, 365–370 (2003).
Google Scholar
Finn, R. D. et al. Pfam: Clans, web tools and services. Nucleic Acids Res. 34, D247–D251 (2006).
Google Scholar
Yu, Z. et al. The new Haemaphysalis longicornis genome provides insights into its requisite biological traits. Genomics 114, 110317 (2022).
Google Scholar
Simao, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V. & Zdobnov, E. M. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31, 3210–3212 (2015).
Google Scholar
Wang, Y. P. et al. MCScanX: a toolkit for detection and evolutionary analysis of gene synteny and collinearity. Nucleic Acids Res. 40, e49 (2012).
Google Scholar
Tang, H. B. et al. jcvi: JCVI Utility Libraries (Zenodo, 2015).
Flynn, J. M. et al. RepeatModeler2 for automated genomic discovery of transposable element families. Proc. Natl. Acad. Sci. USA 117, 9451–9457 (2020).
Google Scholar
Bao, Z. & Eddy, S. R. Automated de novo identification of repeat sequence families in sequenced genomes. Genome Res. 12, 1269–1276 (2002).
Google Scholar
Price, A. L., Jones, N. C. & Pevzner, P. A. De novo identification of repeat families in large genomes. Bioinformatics 21, i351–i358 (2005).
Google Scholar
Ou, S. J. & Jiang, N. LTR_retriever: a highly accurate and sensitive program for identification of long terminal-repeat retrotransposons. Plant Physiol. 176, 1411–1422 (2017).
Ellinghaus, D., Kurtz, S. & Willhoeft, U. LTRharvest, an efficient and flexible software for de novo detection of LTR retrotransposons. BMC Bioinforma. 9, 18 (2008).
Google Scholar
Xu, Z. & Wang, H. LTR_FINDER: an efficient tool for the prediction of full-length LTR retrotransposons. Nucleic Acids Res. 35, W265–W268 (2007).
Google Scholar
Tarailo-Graovac, M. & Chen, N. Using RepeatMasker to identify repetitive elements in genomic sequences. Curr. Protoc. Bioinformatics 25, 10.11–4.10.14 (2009).
Beier, S., Thiel, T., Münch, T., Scholz, U. & Mascher, M. MISA-web: a web server for microsatellite prediction. Bioinformatics 33, 2583–2585 (2017).
Google Scholar
Benson, G. Tandem repeats finder: a program to analyze DNA sequences. Nucleic Acids Res. 27, 573–580 (1999).
Google Scholar
Lowe, T. M. & Eddy, S. R. tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res. 25, 0955–0964 (1997).
Google Scholar
Loman, T. A novel method for predicting ribosomal RNA genes in prokaryotic genomes. BINP30 20161, 8914064 (2017).
Griffiths-Jones, S. et al. Rfam: annotating non-coding RNAs in complete genomes. Nucleic Acids Res. 33, D121–D124 (2005).
Google Scholar
Nawrocki, E. P. & Eddy, S. R. Infernal 1.1: 100-fold faster RNA homology searches. Bioinformatics 29, 2933–2935 (2013).
Google Scholar
Liao, Y., Smyth, G. K. & Shi, W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923–930 (2014).
Google Scholar
Robinson, M. D., McCarthy, D. J. & Smyth, G. K. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140 (2010).
Google Scholar
Yu, G. C., Wang, L. G., Han, Y. Y. & He, Q. Y. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics 16, 284–287 (2012).
Google Scholar
Huo, S. M. et al. Comparative genome and transcriptome analyses reveal innate differences in response to host plants by two color forms of the two-spotted spider mite Tetranychus urticae. BMC Genomics 22, 569 (2021).
Google Scholar
Emms, D. M. & Kelly, S. OrthoFinder: phylogenetic orthology inference for comparative genomics. Genome Biol. 20, 238 (2019).
Google Scholar
Edgar, R. C. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32, 1792–1797 (2004).
Google Scholar
Nguyen, L. T., Schmidt, H. A., von Haeseler, A. & Minh, B. Q. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol. Biol. Evol. 32, 268–274 (2015).
Google Scholar
Yang, Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 24, 1586–1591 (2007).
Google Scholar
Kumar, S., Stecher, G., Suleski, M. & Hedges, S. B. TimeTree: a resource for timelines, timetrees, and divergence times. Mol. Biol. Evol. 34, 1812–1819 (2017).
Google Scholar
Mendes, F. K., Vanderpool, D., Fulton, B. & Hahn, M. W. CAFE 5 models variation in evolutionary rates among gene families. Bioinformatics 36, 5516–5518 (2021).
Google Scholar
Katoh, K. & Standley, D. M. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol. Biol. Evol. 30, 772–780 (2013).
Google Scholar
Price, M. N., Dehal, P. S. & Arkin, A. P. FastTree 2-approximately maximum-likelihood trees for large alignments. PLoS One 5, e9490 (2010).
Google Scholar
Kocot, K. M. et al. PhyloTreePruner: a phylogenetic tree-based approach for selection of orthologous sequences for phylogenomics. Evol. Bioinform. 9, 411–420 (2013).
Google Scholar
Löytynoja, A. Phylogeny-aware alignment with PRANK. Methods Mol. Biol. 1079, 155–170 (2014).
Google Scholar
Zwaenepoel, A. & Van de Peer, Y. wgd-simple command line tools for the analysis of ancient whole-genome duplications. Bioinformatics 35, 2153–2155 (2019).
Google Scholar
Meng, G., Li, Y., Yang, C. & Liu, S. MitoZ: a toolkit for animal mitochondrial genome assembly, annotation and visualization. Nucleic Acids Res. 47, e63 (2019).
Google Scholar
Kearse, M. et al. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28, 1647–1649 (2012).
Google Scholar
Van der Auwera, G. A. et al. From FastQ data to high confidence variant calls: The Genome Analysis Toolkit best practices pipeline. Curr. Protoc. Bioinform. 43, 11.10.1–11.10.33 (2013).
Li, W. & Godzik, A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 22, 1658–1659 (2006).
Google Scholar
Tamura, K., Stecher, G. & Kumar, S. MEGA11: molecular evolutionary genetics analysis version 11. Mol. Biol. Evol. 38, 3022–3027 (2021).
Google Scholar
Garrison, E. & Marth, G. Haplotype-based variant detection from short-read sequencing. arXiv https://arxiv.org/abs/1207.3907 (2012).
Cingolani, P. et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly 6, 80–92 (2012).
Google Scholar
Chen, Z., Yang, X. J., Bu, F. J., Yang, X. L. & Liu, J. Z. Morphological, biological and molecular characteristics of bisexual and parthenogenetic Haemaphysalis longicornis. Vet. Parasitol. 189, 344–352 (2012).
Google Scholar
Greay, T. L. et al. A survey of ticks (Acari: Ixodidae) of companion animals in Australia. Parasit. Vectors 9, 207 (2016).
Google Scholar
Egizi, A. et al. First glimpse into the origin and spread of the Asian longhorned tick, Haemaphysalis longicornis, in the United States. Zoonoses Public Health 67, 637–650 (2020).
Google Scholar
Warren, D. L., Glor, R. E. & Turelli, M. ENMTools: a toolbox for comparative studies of environmental niche models. Ecography 33, 607–611 (2010).
Google Scholar
Phillips, S., Anderson, R., Dudík, M., Schapire, R., & Blair, M. E. Opening the black box: an open-source release of Maxent. Ecography 40, 887–893 (2017).
Cobos, M. E., Peterson, A. T., Barve, N. & Osorio-Olvera, L. kuenm: an R package for detailed development of ecological niche models using Maxent. PeerJ 7, e6281 (2019).
Google Scholar
Acknowledgements
This work was supported by grants from the Natural Science Foundation of China (U24A20363, 32201274), the National Key Research and Development Program of China (2023YFC2305901), the Hebei Province Higher Education Science and Technology Research Project (BJK2023017).
Author information
Authors and Affiliations
Contributions
W.C.C., J.L. and T.W. conceptualized the project and designed the study. W.S., Y.F.W. X.M.C. and L.Z. collected the data and conducted the data analysis. L.F.L., S.J.S., W.J.Z. and W.Y.G. collected samples and carried out experiments. Y.F.W., D.Y.Z. and X.Y.S. prepared materials for sequencing. W.S., J.J.Z., H.W., Z.Y. and performed genome analysis and interpretation. T.W. and W.S. wrote the original draft of the manuscript. W.C.C., N.J. and J.F.J. revised the manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Peer review
Peer review information
Communications Biology thanks Jan Perner and the other, anonymous, reviewers for their contribution to the peer review of this work. Primary Handling Editors: John Mulley and Johannes Stortz.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Supplementary Information
Description of Additional Supplementary Files
Supplementary Data 1
Reporting-summary
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
Reprints and permissions
About this article
Cite this article
Wang, T., Shi, W., Wang, YF. et al. Triploidy drives vector capacity and expansion of the parthenogenetic tick Haemaphysalis longicornis.
Commun Biol 8, 1775 (2025). https://doi.org/10.1038/s42003-025-09153-x
Received:
Accepted:
Published:
Version of record:
DOI: https://doi.org/10.1038/s42003-025-09153-x
Source: Ecology - nature.com
