Skip to main content
Download PDF

Developmental transcriptomics throughout the embryonic developmental process of Rhipicephalus turanicus reveals stage-specific gene expression profiles

Parasites & Vectors volume 15, Article number: 89 (2022) Cite this article

Abstract

Background

Ticks are important vectors and transmit diverse pathogens, including protozoa, viruses, and bacteria. Tick-borne diseases can cause damage to both human health and the livestock industries. The control and prevention of ticks and tick-borne diseases has relied heavily on acaricides.

Methods

In the present study, using a high-throughput RNA sequencing (RNA-Seq) technique, we performed a comprehensive time-series transcriptomic analysis throughout the embryogenesis period of Rhipicephalus turanicus.

Results

Altogether, 127,157 unigenes were assembled and clustered. Gene expression differences among the embryonic stages demonstrated that the most differentially expressed genes (DEGs) were observed in the comparisons of early embryonic stages (RTE5 vs. RTE10, 9726 genes), and there were far fewer DEGs in later stages (RTE25 vs. RTE30, 2751 genes). Furthermore, 16 distinct gene modules were identified according to weighted gene co-expression network analysis (WGCNA), and genes in different modules displayed stage-specific characteristics. Gene Ontology (GO) annotations and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment suggested that some genes involved in organ and tissue formation were significantly upregulated in the early embryonic developmental stages, whereas metabolism-related pathways were more enriched in the later embryonic developmental stages.

Conclusions

These transcriptome studies revealed gene expression profiles at different stages of embryonic development, which would be useful for interrupting the embryonic development of ticks and disrupting the transmission of tick-borne diseases.

Graphical Abstract

Background

Ticks (Acari: Ixodidae) are obligate ectoparasitic arthropods widely distributed across the world. The life-cycle of hard ticks typically includes eggs, larvae, nymphs, and adults (male and female). Blood meals are needed for survival, growth to the next developmental stage, and for reproduction. Then, engorged larvae and nymphs molt into nymphs and adults, respectively. Blood-feeding is a prerequisite for oviposition of adult females, and thousands of eggs can be laid by a single female tick [1]. Bloodsucking of ticks not only harms hosts by causing tick paralysis, allergic reactions, damage to the skin, decreased productivity, and decreased immune function, but also transmits disease-causing pathogens, such as protozoa, viruses, and bacteria, from the infected hosts to other hosts in subsequent blood-feeding [2,3,4,5,6,7,8,9]. Ticks are currently considered to be second only to mosquitoes as vectors of infectious diseases to humans and animals [10].

Ticks have few natural enemies, and vaccines against tick-borne diseases remain unavailable. The control of ticks and tick-borne diseases has largely relied on chemical insecticides. However, long-term and intensive use of insecticides has resulted in the increased resistance of ticks, and chemical residues of insecticides also contaminate the environment and agricultural products [11,12,13,14,15]. Therefore, excessive reliance on pesticides for tick and tick-borne disease control has been proven unsustainable [3, 10, 16,17,18]. The increasing resistance of ticks to these agents and the environmental pollution problems caused by acaricides had led to urgent requirements for novel sustainable methods that will greatly improve both current and future ticks and tick-borne disease control strategies.

Interrupting the life-cycle of ticks might be a feasible strategy to control ticks and block the transmission of tick-borne diseases [19]. Among all of the developmental stages of hard ticks, the embryonic development (from the eggs being laid to the larvae hatching) phase would be an optimal choice. As this stage exists off-host and will last for approximately 30 days, blocking strategies targeting the eggs will not affect their hosts. Additionally, embryonic development is vital for the whole life-cycle of ticks and involves an exclusive sequence of cellular and molecular processes that are not present in other stages [19]. Therefore, the identification of crucial molecular targets for the interruption of embryonic development would contribute to reducing the size of tick populations and would be helpful for preventing the transmission of both transstadial (stage-to-stage, also called horizontal) and transovarial (female-to-egg, also called vertical) tick-borne pathogens.

Rhipicephalus turanicus is widely distributed in Africa, Asia, and the Mediterranean regions [20]. Hosts of this tick species include goats, dogs, cattle, sheep, lions, and humans [20, 21]. This tick species has been implicated as a vector of several human and veterinary pathogens, such as Rickettsia spp. [22, 23], Crimean–Congo hemorrhagic fever virus, West Nile virus, Babesia, Theileria, Anaplasma, and Hepatozoon [20, 24,25,26,27]. Due to the lack of genomic resources and systematic research, the molecular mechanisms of the embryonic development of R. turanicus remain largely unknown. In this study, based on a high-throughput RNA sequencing (RNA-Seq) technique, we performed comprehensive time-series transcriptomic analysis of eggs of R. turanicus to provide global insights into the dynamics of gene expression during the successive stages of embryonic development and to identify development-related genes across different stages.

Methods

Tick collection and rearing

Rhipicephalus turanicus specimens were collected from grasslands in Tai’an, Shandong Province, and have been maintained in the laboratory since 2018 (Additional file 1). Ticks were fed on the ears of domestic rabbits for blood meals and were reared in an incubator at 70–80% relative humidity (RH) under a 12/12 h light/dark (L/D) photoperiod at 25 ± 1 °C after they were detached from their hosts. The eggs used in this study were collected from engorged female ticks (descendant of the same female R. turanicus) that fed on the same rabbit. Engorged females were examined daily to ensure that eggs could be collected in time. Laid eggs were collected every day from the start of oviposition and were kept in different plastic tubes covered with cotton mesh (to reduce water loss from evaporation and maintain surface gas exchange). The newly collected eggs were marked as 1 day, eggs of 5 days (RTE5), 10 days (RTE10), 15 days (RTE15), 20 days (RTE20), 25 days (RTE25), and 30 days (RTE30) were chosen, and approximately 300 mg of eggs (in three biological replicates) for each day was used for transcriptome analyses.

RNA extraction, library construction, and sequencing

Total RNA was extracted using the RNAiso Plus reagent (Takara, Japan). RNA degradation and contamination were checked with denaturing agarose gel electrophoresis and RNA integrity was assessed using the RNA 6000 Nano assay kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Altogether 18 transcriptome sequencing libraries from six groups (with three biological replicates for each group) of eggs were constructed using a NEBNext® Ultra™ RNA Library Prep Kit for Illumina (NEB, USA) according to the manufacturer’s instructions. High-throughput transcriptome sequencing was performed on an Illumina HiSeq platform, and paired-end reads (raw reads) with a length of 150 base pairs (bp) were generated.

RNA-Seq data analysis

The quality and number of reads of each sample were assessed using FastQC version 0.11.4 [28]. Raw reads in FASTQ format were trimmed using Trimmomatic version 0.38 [29]. Clean data (clean reads) were obtained by removing adapter, poly-N and low-quality reads from the raw data. De novo transcriptome assembly was carried out using Trinity [30]. The Benchmarking Universal Single-Copy Orthologs (BUSCO) v3.0.2 program was utilized to assess the quality and annotation completeness of the genome assembly [31]. Then assembled unigenes were annotated by searching against the NCBI non-redundant protein sequences (Nr) and the Swiss-Prot database using BLAST alignment with an E-value cutoff of 1e−5. Gene functions were annotated by searching against the Pfam database [32] using the HMMER3 program [33]. Functional annotation by Gene Ontology (GO) terms was achieved by Blast2GO software (version 2.3.5, https://www.blast2go.com/) [34]. The unigenes were also aligned with the Cluster of Orthologous Groups (COG) database to classify and predict functions [35]. The best hits were used to determine the sequence orientations and coding sequences (CDSs) of the unigenes. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed using the KEGG Automatic Annotation Server (KAAS) [36].

Differential gene expression analyses

The number of fragments per kilobase of the transcript sequence per million mapped reads (FPKM) was used to estimate the relative expression levels of unigenes. Differential expression analysis between groups was performed using the DESeq2 R package (version 1.20.0) with default normalization. P-values of the results were adjusted using the Benjamini and Hochberg approach for controlling the false discovery rate. Besides the DESeq2 package, significant differentially expressed genes (DEGs) between two compared groups were also determined using the edgeR package (version 3.12.1) with a fold change threshold of 2 and a false discovery rate (FDR) threshold of 0.05. Genes with an adjusted P-value < 0.05 were assigned as DEGs. Results from both packages were used for DEGs analyses to ensure the reliability of DEG assignment. DEGs of all six groups were subjected to principal component analysis (PCA).

Weighted gene co-expression network analysis (WGCNA)

WGCNA (v1.47) was used to define co-expressed modules and hub genes. Based on pairwise correlations between genes, different modules can be divided into genes with similar expression profiles grouped into the same module [37, 38]. Each module was summarized by a single representative expression profile, which was referred to as the module eigengene. A hierarchical clustering tree was used for module classification according to the dynamic tree cut method. Then, modules with similar expression profiles were merged as merged dynamics. Gene connectivity was represented by edge weight and defined as the sum of weights across all edges of a node in the gene co-expression network analysis. Hub genes were defined based on the intra-modular connectivity which was identified by the CytoHubba plugin in Cytoscape v3.7.2.

Quantitative real-time PCR (qRT-PCR) validation

To validate the RNA-Seq data, the expression levels of eight DEGs involved in embryonic developmental processes were randomly chosen and confirmed by qRT-PCR. Specific primers for these selected genes were designed by Primer Premier 6.0 (Premier Biosoft International, CA, USA) and are listed in Additional file 10: Table S1. Total RNA was extracted from the same batch of samples as those used in RNA-Seq, and reverse transcription was performed using a PrimeScript RT Reagent Kit with gDNA Eraser (TaKaRa, Dalian, China). RT-PCR was performed on an ABI 7500 system (ABI, CA, USA). Individual reactions were prepared with 100 ng of complementary DNA (cDNA) and a ChamQ SYBR qPCR Master Mix Kit (Cat#Q311-02, Vazyme, Nanjing, China) in a final volume of 20 μl. qRT-PCR reactions were carried out for 15 min at 37 °C, followed by 5 s at 85 °C, then followed by 40 cycles of two-step PCR for 10 s at 95 °C and 30 s at 60 °C. All reactions were carried out in triplicate for each sample. β-actin was used as an endogenous control to normalize expression levels. Cycle threshold (Ct) values were normalized using the 2−ΔΔCt method [39], and one-way analysis of variance (ANOVA) was adopted to verify the significant difference (P ≤ 0.05).

Results

Sequencing, assembly, and annotation of the R. turanicus transcriptome

After quality trimming and filtering, 37,230,938 to 61,332,076 clean reads were obtained for 18 different libraries, with Q30 > 94% for all samples (Additional file 11: Table S2). Combined assembly based on total reads from all samples generated 127,157 unigenes and the N50 was 2356 bp. The average length of the unigenes was 1044 bp, and the shortest and longest sequences were 201 bp and 28,091 bp, respectively. For these unigenes, 33,151 were annotated, of which 33,107 (26%), 23,674 (18.6%), 13,354 (10.5%) and 13,793 (10.8%) unigenes were successfully annotated against known genes in the NR, KEGG, KOG, and SwissProt databases, respectively. A large number of mapped unigenes were matched with Ixodes scapularis (Additional file 2: Figure S1). BUSCO revealed 94.2% complete, 76% complete and single copy, 18.2% complete and duplicated, 3.9% fragmented, and 1.9% missing BUSCOs.

Unigene functional classification

Functional classification based on GO demonstrated that annotated unigenes were divided into 58,986 GO terms and classified into 58 groups in three major GO categories (Additional file 3: Figure S2). The cellular process (7443 unigenes) was predominant in the biological process ontology, followed by metabolic processes (6879 unigenes) and single-organism processes (6553 unigenes). In the cellular component ontology, cell (4611 unigenes) and cell part (4,611 unigenes) exhibited the same abundance. Catalytic activity (6536 unigenes) and binding (5152 unigenes) were dominant in the molecular function ontology. All unigenes obtained from COG annotation (13,354 unigenes) were classified into 25 families (Additional file 4: Figure S3). Signal transduction mechanisms (2246 unigenes) were the most enriched family, followed by posttranslational modification, protein turnover, chaperones (1266 unigenes), and transcription (954 unigenes). Then, all assembled unigenes were mapped to the KEGG database to identify pathways in which genes were involved. All 23,674 unigenes were aligned to 145 pathways, and the most representative pathways were metabolic pathways (159 unigenes), lysosome (322 unigenes), and RNA transport (229 unigenes) (Additional file 12: Table S3).

Comparison of gene expression profiles among different groups

Based on the PCA results, the samples were differentiated according to embryonic development status (Fig. 1). Groups 1, 2, and 3 clustered together, group 4 clustered separately, and groups 5 and Group 6 had more similar gene expression patterns than those of other groups.

Fig. 1
figure 1

Principal component analysis (PCA) plot displays distance for transcripts of different samples

Full size image

The FPKM method was used to assess the expression level and normalize unigene in each sample. Only genes that were differentially expressed in both edgeR and DESeq2 packages were selected for further analysis. DEGs between every two consecutive groups were identified. According to the comparative analysis, more genes were upregulated than downregulated in RTE10 vs. RTE15 (4,601/2,171 up-/downregulated), RTE15 vs. RTE20 (4,160/2,230 up-/downregulated), RTE20 vs. RTE25 (2,104/1,657 up-/downregulated), and RTE25 vs. RTE30 (1,523/1,228 up-/downregulated); while downregulated genes substantially exceeded upregulated genes in RTE5 vs. RTE10 (3,193/6,533 up-/downregulated) (Fig. 2). Comparisons between every two consecutive groups suggested that most of the unique DEGs were found in the comparison between RTE5 and RTE10 (5,949 genes), and only 696 unique DEGs were found between RTE25 and RTE30 (Fig. 3). Taken together, 108 DEGs were shared by all five comparisons of six time points.

Fig. 2
figure 2

Numbers of DEGs between different embryonic developmental stages of R. turanicus. Red bars represent upregulated and blue bars represent downregulated genes. The x-axis indicates comparison groups; the y-axis denotes the number of genes

Full size image
Fig. 3
figure 3

Venn diagram showing the number of DEGs expressed in different comparison groups

Full size image

GO enrichment analysis was carried out to investigate the functions of annotated DEGs. Among the annotated DEGs, the biological processes category contained the most GO annotations (7651 genes), followed by molecular function (7497 genes) and cellular components (4944 genes) (Additional file 13: Table S4). Enriched biological processes of RTE5/RTE10 were generally associated with DNA binding (GO: 0003677), nucleic acid binding transcription factor activity (GO: 0001071), and regulation of transcription from the RNA polymerase II promoter (GO: 0006357). Enriched biological processes of RTE10/RTE15 were DNA binding (GO: 0003677), gated channel activity (GO: 0022836), and ligand-gated ion channel activity (GO: 0015276). Receptor activity (GO: 0004872), transmembrane receptor activity (GO: 0099600), and molecular transducer activity (GO: 0060089) were the dominant biological processes in RTE15/RTE20. Fatty acid metabolic process (GO: 0006631) was the most significantly enriched GO term in both RTE20/RTE25 and RTE25/RTE30 (Additional file 13: Table S5).

The functional classifications and predictions were further analyzed by searching DEGs against the KEGG database to determine the involvement of biological pathways in the embryogenesis of R. turanicus. There were 140 KEGG pathways were predicted in RTE5/RTE10 and RTE10/RTE15, respectively (Additional file 14: Table S6). The most enriched pathways of RTE5/RTE10 included axon regeneration (ko04361), Wnt signaling pathway (ko04310), and MAPK signaling pathway-fly (ko04013) (Additional file 5: Figure S4; Additional file 14: Table S6). DNA replication (ko03030), neuroactive ligand–receptor interaction (ko04080), and pyrimidine metabolism (ko00240) were highly enriched in RTE10/RTE15 (140 pathways). The representative pathway in RTE15/RTE20 (136 pathways) was neuroactive ligand–receptor interaction (ko04080). Biosynthesis of unsaturated fatty acids (ko01040), fatty acid elongation (ko00062), and fatty acid metabolism (ko01212) were the top three highly enriched pathways in both RTE20/RTE25 (126 pathways) and RTE25/RTE30 (119 pathways). Two drug metabolism pathways were involved, including drug metabolism-other enzymes (ko00983) in RTE5/RTE10 and RTE15/RTE20 and drug metabolism-cytochrome P450 (ko00982) enrichment in RTE10/RTE15, RTE20/RTE25, and RTE25/RTE30 (Additional file 5: Figure S4; Additional file 14: Table S6).

Correlation of differentially expressed genes and WGCNA

WGCNA identified modules of co-expressed genes and candidate hub genes for each time point. After modules with similar expression profiles were merged, 16 distinct gene modules were identified (Additional file 6: Figure S5). The number of genes in each module were 4748 (blue), 3276 (brown), 3086 (green), 2255 (darkgrey), 1154 (magenta), 697 (darkgreen), 390 (darkred), 369 (tan), 296 (lightyellow), 244 (midnightblue), 231 (lightcyan), 185 (royalblue), 122 (orange), 114 (white), and 57 (steelblue). Unassigned genes (22 genes) were placed into the “grey” module (Additional file 7: Figure S6).

The genes in the module blue, brown, and darkgrey were upregulated in the early stages of embryogenesis (RTE5, RTE10), the lightcyan module consisted of genes upregulated in RTE10 and RTE20 of embryogenesis, while genes in the green and tan module were merely upregulated in RTE25 and RTE30 (Fig. 4). The upregulated expression patterns of the darkred, lightyellow, and royalblue modules were maintained in both RTE15 and RTE20. Additionally, the genes of the orange and steelblue modules were upregulated in RTE15 and RTE30, respectively. The genes of the darkgreen module were mainly upregulated in RTE20, while the genes of the white module were upregulated in RTE5 and RTE30, respectively. The magenta module genes were upregulated in RTE10 and RTE15, whereas the midnightblue genes were upregulated in RTE20 and RTE25 (Additional file 8: Figure S7).

Fig. 4
figure 4

Expression pattern (histogram) and corresponding KEGG enriched pathway (scatter plot) analysis of blue, brown, darkgrey, darkgreen, green, and tan modules. In the histogram, the x-axis represents the sample, the y-axis represents the expression profile of the eigengene; the x-axis and y-axis in the scatter plot are the rich factor and KEGG pathways. Rich factor refers to the ratio of the number of genes located in the KEGG pathway and the total number of genes in the KEGG pathway. The larger the rich factor, the greater enrichment.

Full size image

GO analysis and KEGG pathway enrichment were performed to explore potential biological processes associated with different embryogenic developmental stages. GO term enrichment analysis revealed that the genes of the blue, brown, darkgrey, darkgreen, tan, and green modules were related to three main categories (biological processes, molecular function, cellular components). The most enriched GO subcategories of these six modules were cellular process, metabolic process, single-organism process, cell, cell apart, catalytic activity, and binding, respectively (Additional file 9: Figure S8).

The most enriched KEGG pathways in the blue module included oxidative phosphorylation (ko00190), aminoacyl–transfer RNA (tRNA) biosynthesis (ko00970), and proteasome (ko03050). Spliceosome (ko03040), DNA replication (ko03030), and messenger RNA (mRNA) surveillance (ko03015) were the top three pathways in the brown module. The Wnt signaling pathway (ko04310), Notch signaling pathway (ko04330), and endocytosis (ko04144) were abundant in the darkgrey module. In the darkgreen module, lysosome (ko04142), glutathione metabolism (ko00480), biosynthesis of unsaturated fatty acids (ko01040), and peroxisome (ko04146) were significantly enriched. Biosynthesis of unsaturated fatty acids (ko01040), fatty acid elongation (ko00062), and fatty acid metabolism (ko01212) were predominant in the tan module. The green module was enriched with neuroactive ligand–receptor interaction (ko04080), calcium signaling pathway (ko04020), and nitrogen metabolism (ko00910) (Fig. 4).

Networks were constructed to explore relationships among genes, in which each node represents a gene and the edges between genes represent co-expression correlations. To capture a more meaningful correlation, genes with the highest connectivity were selected as the hub genes in each module (Fig. 5). The hub genes of the blue, brown, darkgrey, and tan modules were translation initiation factor 3 subunit (v1g168210), synaptic functional regulator FMR1 (Fxr1), rac GTPase-activating protein 1-like (RACGAP1), and tyrosine-protein phosphatase 10D (Ptp10D), respectively. However, the hub genes of the green and darkgreen modules are still unannotated (Additional file 15: Table S7).

Fig. 5
figure 5

Networks of genes in different modules. The top ten hub genes are shown in red in the blue, darkgrey, darkgreen, green, and tan modules; the top ten hub genes are shown in green in the brown module

Full size image

Quantitative real-time PCR analysis

The expression profiles of 8 randomly selected genes (Nop60B, unigen0086511, SRSF7, HSP90AA1, NOP58, rbm4.1, UBA1, CAM) were quantified using qRT-PCR to validate the accuracy and reproducibility of the expression profiles of RNA-Seq. The results from the qRT-PCR experiments were consistent with those observed in transcriptome analysis, and both of the strategies revealed similar trends in the upregulated or downregulated genes (Fig. 6).

Fig. 6
figure 6

Comparison of expression values of eight genes by RNA-Seq and real-time RT-PCR analyses. The x-axis represents the developmental time point; the y-axis represents the relative expression levels. Asterisks indicate statistical significance at the level of P < 0.01

Full size image

Discussion

Blood-sucking arthropods are responsible for major epidemics throughout the world, causing thousands of deaths and great harm to public health. Ticks affect both humans and animals, causing considerable losses to the cattle industry worldwide and incurring substantial damage to livestock [40]. At each blood meal, ticks can become integrated into the epidemiological chain of pathogen transmission by means of so-called transstadial or transovarial passages.

The cattle tick R. turanicus is an obligate hematophagous arthropod, and embryogenesis is one of the most important life stages, as it involves large changes to the RNA transcript profile [41]. Identifying transcriptional changes during embryogenesis is of crucial importance for unraveling the evolutionary, molecular, and cellular mechanisms that underpin patterning and morphogenesis [42]. Most of the knowledge about the developmental gene regulation of Arthropoda has been derived from the Drosophila paradigm [43]. However, the involvement of genes in specific developmental processes is usually influenced by many factors, including species and the expression patterns of genes. Developmental transcriptomics across all embryonic stages can provide gene expression divergence information spanning different developmental time points and offer novel insights into the genetics of embryogenesis.

Numerous efforts have been made in basic and applied tick and tick-borne disease control research [18, 44,45,46,47]. However, as one of the important stages of the tick life-cycle, research on embryos is limited, and available studies mainly focus on the morphological, histological, and ultrastructural characteristics of embryonic development [48, 49]. This paucity of early-stage information hinders the implementation of targeted approaches, such as RNA interference or characterization of vaccine candidates [50, 51]. Our study is among the first to provide a comprehensive and biologically relevant catalog of transcripts for future research aiming at controlling the population of ticks in the early stages of their development [52].

In the present study, the transcriptomes of six different developmental points that were determined according to the days after oviposition and covering the entire embryonic developmental process of R. turanicus were compared. PCA demonstrated that a few individuals from the same developmental stage were not closely clustered (Fig. 1), indicating that there were variations among individuals from the same embryonic stage [41]. Close clustering of Group 1, 2, and 3 indicated that the gene expression profiles were very similar during the early stages (RTE5 to RTE15) of embryogenesis, whereas the divergence of Group 4 from all other groups suggested that great changes in the gene expression patterns occurred at this time point (Fig. 1). According to the variations in gene expression, the changes between RTE5 and RTE10 were more intense than those between RTE25 and RTE30 (Fig. 2), confirming that more dramatic changes occurred during the earlier embryonic developmental periods.

Embryogenesis is a complex process involving an elaborate network of signaling pathways. Comparisons of gene expression among different stages revealed that a relatively high proportion of DEGs were present during early embryonic development stages, suggesting that more intense changes occur during the earlier developmental periods. The GO term and KEGG enrichment pathway analyses showed that the annotated DEGs participated in multiple physiological processes and that the expression levels of fundamental development-related genes were high during the early embryonic developmental stages (Additional file 3: Figure S2; Additional file 5: Figure S4; Additional File 15: Table S6). However, DESeq pairwise differential expression analysis might miss some expression profiles using only pairwise comparisons. WGCNA classifies the whole transcriptome into different expression patterns based on the comparison of all genes in all samples, which can help to obtain a whole picture of the expression profiles of embryogenesis. Genes involved in the same biological process were highly correlated among samples, and co-expressed genes throughout different embryonic developmental stages of R. turanicus clustered into 16 modules (Additional file 6: Figure S5), which provided more detailed gene expression dynamics and enrichment of pathway profiles according to embryonic developmental stages.

Limited overlap of these genes across embryonic developmental stages was observed, and this is consistent with previous phenotypic reports [48]. Stage-specific transcript characteristics were evidenced by the functional categories of the modules to which they belonged (Fig. 4). The main biological processes (e.g., RNA transport, DNA replication, and mRNA surveillance) in which the genes were upregulated in the blue and brown modules indicated that cell division is particularly intense during RTE5 and RTE10. These results are consistent with the results of morphogenetic processes, which suggested that intense cell proliferation occurred during the early stages of tick embryonic development [48]. Moreover, the significant enrichment of aminoacyl–tRNA biosynthesis, ribosome, and spliceosome gave a signature of intense protein production that could be related to secretory cells. Therefore, we can deduce that salivary glands might be formed in this stage.

The Wnt signaling pathway regulates a huge variety of embryonic developmental processes by inducing transcriptional or morphological changes in responding cells [53, 54]. Embryonic development analysis of Rhipicephalus microplus confirmed the expression of the Wnt signaling pathway during early embryonic development [48]. The enriched Wnt signaling pathway in the darkgrey module in this study indicated that the RTE10 and RTE15 stages may correspond to stages 4 and 5 of R. microplus, in which cumulus cells and axis establishment were detected. In addition, fundamental developmental genes related to the Notch [55], Hippo [56], Hedgehog [57, 58], transforming growth factor beta [59,60,61], and mitogen-activated protein kinase (MAPK) [62] signaling pathways showed higher expression during RTE10 and RTE15 stages than others (Fig. 4). Therefore, these stages are characterized by more active biological processes, including organ and tissue formation.

Lysosomes are membrane-delimited organelles that contain strongly hydrolytic enzymes and serve as the main digestive compartment of cells [63, 64]. Peroxisomes are essential organelles of eukaryotic cells that play a key role in the oxidation of fatty acids and the generation and removal of hydrogen peroxide [65]. Both lysosomes and peroxisomes are immune-related pathways, and the high enrichment of these pathways in the darkgreen module suggested that R. turanicus embryos have a relatively complete immune system at RTE20. Moreover, the high expression of immune-related pathways also indicated that immune defense is vital to this embryonic developmental stage. The embryogenesis of the tick is typically described as an energy-consuming process [40]. Energy metabolism is essential to supporting the molecular backbone needed for cell proliferation, differentiation, and embryonic growth. The catabolism of biomolecules, such as carbohydrates, is the main source of energy for the developing embryo [66]. These studies elucidated the highly enriched metabolism-related pathways during the later embryonic developmental stages (RTE25, RET30). Furthermore, the biosynthesis of unsaturated fatty acids and fatty acid metabolism, which are probably involved in exoskeleton structure formation (e.g., wax synthesis), were found in the later stages of embryonic development of ticks [48]. These results revealed that gene expression dynamics and patterns are highly varied among different embryonic developmental stages. However, highly enriched energy metabolism pathways, including oxidative phosphorylation in the blue module, sulfur metabolism in the tan module, and nitrogen metabolism in the green module, crossing multiple time points confirmed that embryogenesis of ticks is an energy-consuming process [40].

Additionally, drug metabolism-cytochrome P450 and drug metabolism-other enzymes maintained high expression during all stages (Fig. 4; Additional file 14: Table S6). Cytochrome P450s constitute a diverse and important gene superfamily in all organisms and are known to catalyze a diverse range of chemical reactions important for both developmental processes and the detoxification of exogenous compounds [67]. The overexpression of drug metabolism pathways suggests the high metabolic level of R. turanicus embryos to hazardous substances, while it is difficult to illustrate whether endogenous or exogenous substances induced the high expression levels of cytochrome P450 and related enzymes in the present study. Since all egg samples used in this study were kept in sterile plastic tubes and covered with cotton mesh, sterilized water was added onto the cotton mesh to keep the environment moist, as needed for the incubation of eggs. Additionally, ancestors of these eggs have been maintained in the laboratory since 2018, and they have not been in contact with any harmful substances during this time. Therefore, we speculate that endogenous products, which cannot be released in time due to the existence of wax in eggs, caused the high expression of drug metabolism pathways. Further research is needed to gain a better understanding of the role of drug metabolism pathways in the embryonic developmental stages of ticks and to identify possible new targets for the development of novel control measures against ticks.

Co-expression networks of each module were constructed and hub genes were identified accordingly (Fig. 5; Additional file 15: Table S7). However, as most of the candidate hub genes are functionally uncharacterized, it is currently impossible to correlate these hub genes with crucial biological functions in the embryonic developmental processes of R. turanicus. Directly linked genes usually have similar functions and may be involved in the same biological pathways related to embryonic development (Fig. 5; Additional file 15: Table S7). Therefore, the present research revealed new insights for further studies, as the hub genes may act as potential targets for further exploration of detailed biological functions to shed light on possibilities for the control of ticks and tick-borne diseases.

Conclusion

In this study, we recapitulated the time-series transcriptomic profile of embryonic development of R. turanicus. The results demonstrated that not all genes are expressed at once, but each has its own specific temporal expression pattern. WGCNA classified all DEGs into 16 modules, and co-expression networks provided clue candidate hub genes of each module. KEGG enrichment pathways suggested that the most significant expression changes occur during the early embryonic developmental stages, particularly in RTE10 and RTE15, in which a large set of genes involved in organ and tissue formation are upregulated. In the later embryonic developmental stages, more genes were involved in metabolism-related pathways. These different expression profiles should be taken into consideration in the further exploration of target genes to interrupt the embryonic development of ticks.

Availability of data and materials

Raw RNA-Seq sequences were deposited in NCBI Sequence Read Archive (SRA) database with accession number: PRJNA741544.

References

  1. Apanaskevich DA, Oliver JH. Life cycles and natural history of ticks. In: Sonenshine DE, Roe RM, editors. Biology of ticks, vol. 1. 2nd ed. New York: Oxford University Press; 2014. p. 59–73.

    Google Scholar 

  2. Dennis DT, Piesman JF. Overview of tick-bore infections of humans. In: Goodman JL, Dennis DT, Sonenshine DE, editors. Tick-borne diseases of humans. Washington D.C: ASM Press; 2005. p. 3–11.

    Google Scholar 

  3. Ghosh S, Azhahianambi P, Yadav MP. Upcoming and future strategies of tick control: a review. J Vector Borne Dis. 2007;44:79–89.

    CAS  PubMed  Google Scholar 

  4. Elango G, Rahuman AA. Evaluation of medicinal plant extracts against ticks and fluke. Parasitol Res. 2011;108:513–9.

    Article  PubMed  Google Scholar 

  5. Pfäffle M, Littwin N, Muders SV, Petney TN. The ecology of tick-borne diseases. Int J Parasitol. 2013;43:1059–77.

    Article  PubMed  Google Scholar 

  6. Sonenshine DE, Roe RM. Overview: ticks, people and animals. In: Sonenshine DE, Roe RM, editors. Biology of ticks, vol. 1. New York: Oxford University Press, New York; 2014. p. 3–16.

    Google Scholar 

  7. Cabezas-Cruz A, Mateos-Hernández L, Pérez-Cruz M, Valdés JJ, Mera IGF, Villar M, et al. Regulation of the immune response to α-gal and vector-borne diseases. Trends Parasitol. 2015;31:470–6.

    Article  CAS  PubMed  Google Scholar 

  8. Fang LQ, Liu K, Li XL, Liang S, Yang Y, Yao HW, et al. Emerging tick-borne infections in mainland China: an increasing public health threat. Lancet Infect Dis. 2015;15:1467–79.

    Article  PubMed  PubMed Central  Google Scholar 

  9. Medlock JM, Leach SA. Effect of climate change on vector-borne disease risk in the UK. Lancet Infect Dis. 2015;15:721–30.

    Article  PubMed  Google Scholar 

  10. de la Fuente J, Estrada-Pena A, Venzal JM, Kocan KM, Sonenshine DE. Overview: Ticks as vectors of pathogens that cause disease in humans and animals. Front Biosci. 2008;13:6938–46.

    Article  PubMed  Google Scholar 

  11. Graf JF, Gogolewski R, Leach-Bing N, Sabatini GA, Molento MB, Bordin EL, et al. Tick control: an industry point of view. Parasitol. 2004;129:S427–42.

    Article  Google Scholar 

  12. Rosario-Cruz R, Almazán C, Miller RJ, Dominquez-Garcia DI, Hernandez-Ortiz R, de la Fuente J. Genetic basis and impact of tick acaricide resistance. Front Biosci. 2009;14:2657–65.

    Article  CAS  Google Scholar 

  13. de la Fuente J, Kocan KM. Development of vaccines for control of tick infestations and interruption of pathogen transmission. In: Sonenshine DE, Roe RM, editors. Biology of ticks, vol. 2. 2nd ed. Oxford, UK: Oxford University Press; 2014. p. 333–52.

    Google Scholar 

  14. de León AAP, Teel PD, Auclair AN, Messenger MT, Guerrero FD, Schuster G, et al. Integrated strategy for sustainable cattle fever tick eradication in USA is required to mitigate the impact of global change. Front Physiol. 2012;3:195.

    Google Scholar 

  15. Abbas RZ, Zaman MA, Colwell DD, Gilleard J, Iqbal Z. Acaricide resistance in cattle ticks and approaches to its management: the state of play. Vet Parasitol. 2014;203:6–20.

    Article  CAS  PubMed  Google Scholar 

  16. Rodriguez-Vivas RI, Trees AJ, Rosado-Aguilar JA, Villegas-Perez SL, Hodgkinson JE. Evolution of acaricide resistance: phenotypic and genotypic changes in field populations of Rhipicephalus (Boophilus) microplus in response to pyrethroid selection pressure. Int J Parasitol. 2011;41:895–903.

    Article  CAS  PubMed  Google Scholar 

  17. Parizi LF, Githaka NW, Logullo C, Konnai S, Masuda A, Ohashi K, et al. The quest for a universal vaccine against ticks: cross-immunity insights. Vet J. 2012;194:158–65.

    Article  PubMed  Google Scholar 

  18. Merino O, Alberdi P, Pérez de la Lastra JM, De la Fuente J. Tick vaccines and the control of tick-borne pathogens. Front Cell Infect Microbiol. 2013;3:30.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Silveira AB, Castro-Santos J, Senna R, Logullo C, Fialho E, Silva-Neto MA. Tick vitellin is dephosphorylated by a protein tyrosine phosphatase during egg development: effect of dephosphorylation on VT proteolysis. Insect Biochem Mol Biol. 2006;36:200–9.

    Article  CAS  PubMed  Google Scholar 

  20. Walker JB, Keirans JE, Horak IG. The genus Rhipicephalus (Acari, Ixodidae): a guide to the brown ticks of the world. Cambridge: Cambridge Academic Press; 2000. p. 656.

    Book  Google Scholar 

  21. Horak IG, Heyne H, Williams R, Gallivan GJ, Spickett AM, Bezuidenhout JD, et al. The Ixodid ticks (Acari: Ixodidae) of Southern Africa. Pretoria: Springer; 2018.

    Book  Google Scholar 

  22. Matsumoto K, Ogawa M, Brouqui P, Raoult D, Parola P. Transmission of Rickettsia massiliae in the tick, Rhipicephalus turanicus. Med Vet Entomol. 2005;19:263–70.

    Article  CAS  PubMed  Google Scholar 

  23. Germanakis A, Chochlakis D, Angelakis E, Tselentis Y, Psaroulaki A. Rickettsia aeschlimannii infection in a man. Greece Emerg Infect Dis. 2013;19:1176–7.

    Article  PubMed  Google Scholar 

  24. Vitale G, Mansuelo S, Rolain J-M, Raoult D. Rickettsia massiliae human isolation. Emerg Infect Dis. 2006;12:174–5.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Mura A, Masala G, Tola S, Satta G, Fois F, Piras P, et al. First direct detection of rickettsial pathogens and a new rickettsia, 'Candidatus Rickettsia barbariae’, in ticks from Sardinia. Italy Clin Microbiol Infect. 2008;14:1028–33.

    Article  CAS  PubMed  Google Scholar 

  26. Giannelli A, Lia RP, Annoscia G, Buonavoglia C, Lorusso E, Dantas-Torres F, et al. Rhipicephalus turanicus, a new vector of Hepatozoon canis. Parasitology. 2017;144:730–7.

    Article  CAS  PubMed  Google Scholar 

  27. Belkahia H, Ben Said M, Ghribi R, Selmi R, Ben Asker A, Yahiaoui M, et al. Molecular detection, genotyping and phylogeny of Anaplasma spp. in Rhipicephalus ticks from Tunisia. Acta Trop. 2019;191:38–49.

    Article  PubMed  Google Scholar 

  28. Andrews S. FASTQC. A quality control tool for high throughput sequence data. Babraham Bioinformatics, Babraham Institute, Cambridge, United Kingdom. 2010.

  29. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29:644–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31:3210–2.

    Article  PubMed  Google Scholar 

  32. Finn RD, Coggill P, Eberhardt RY, Eddy SR, Mistry J, Mitchell AL, et al. The Pfam protein families database: towards a more sustainable future. Nucleic Acids Res. 2016;44:D279–85.

    Article  CAS  PubMed  Google Scholar 

  33. Eddy SR. A new generation of homology search tools based on probabilistic inference. Genome Inform. 2009;23:205–11.

    PubMed  Google Scholar 

  34. Aparicio G, Götz S, Conesa A, Segrelles D, Blanquer I, García JM, et al. Blast2GO goes grid: developing a grid-enabled prototype for functional genomics analysis. Stud Health Technol Inform. 2006;120:194–204.

    CAS  PubMed  Google Scholar 

  35. Tatusov RL, Fedorova ND, Jackson JD, Jacobs AR, Kiryutin B, Koonin EV, et al. The COG database: an updated version includes eukaryotes. BMC Bioinformatics. 2003;4:41.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 2007;35:W182–5.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Zhang B, Horvath S. A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol. 2005;4: Article17.

  38. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2 (-Delta Delta C (T)) Method. Methods. 2001;25:402–8.

    Article  CAS  PubMed  Google Scholar 

  40. Martins R, Ruiz N, Fonseca RND, Vaz Junior IDS, Logullo C. The dynamics of energy metabolism in the tick embryo. Rev Bras Parasitol Vet. 2018;27:259–66.

    CAS  PubMed  Google Scholar 

  41. Liu HQ, Li Y, Irwin DM, Zhang YP, Wu DD. Integrative analysis of young genes, positively selected genes and lncRNAs in the development of Drosophila melanogaster. BMC Evol Biol. 2014;14:241.

    Article  PubMed  PubMed Central  Google Scholar 

  42. Simon S, Sagasser S, Saccenti E, Brugler MR, Schranz ME, Hadrys H, et al. Comparative transcriptomics reveal developmental turning points during embryogenesis of a hemimetabolous insect, the damselfly Ischnura elegans. Sci Rep. 2017;7:13547.

    Article  PubMed  PubMed Central  Google Scholar 

  43. Peel AD, Chipman AD, Akam M. Arthropod segmentation: beyond the Drosophila paradigm. Nat Rev Genet. 2005;6:905–16.

    Article  CAS  PubMed  Google Scholar 

  44. Labuda M, Trimnell AR, Licková M, Kazimírová M, Davies GM, Lissina O, et al. An antivector vaccine protects against a lethal vector-borne pathogen. PLoS Pathog. 2006;2:e27.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Nuttall PA, Trimnell AR, Kazimirova M, Labuda M. Exposed and concealed antigens as vaccine targets for controlling ticks and tick-borne diseases. Parasite Immunol. 2006;28:155–63.

    Article  CAS  PubMed  Google Scholar 

  46. Sprong H, Trentelman J, Seemann I, Grubhoffer L, Rego ROM, Hajdušek O, et al. ANTIDotE: anti-tick vaccines to prevent tick-borne diseases in Europe. Parasit Vectors. 2014;7:77.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Rego ROM, Trentelman JJA, Anguita J, Nijhof AM, Sprong H, Klempa B, et al. Counterattacking the tick bite: towards a rational design of anti-tick vaccines targeting pathogen transmission. Parasit Vectors. 2019;12:229.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Santos VT, Ribeiro L, Fraga A, de Barros CM, Campos E, Moraes J, et al. The embryogenesis of the tick Rhipicephalus (Boophilus) microplus: the establishment of a new chelicerate model system. Genesis. 2013;51:803–18.

    Article  CAS  PubMed  Google Scholar 

  49. Friesen KJ, Dixon M, Lysyk TJ. Embryo development and morphology of the rocky mountain wood tick (Acari: Ixodidae). J Med Entomol. 2016;53:279–89.

    Article  CAS  PubMed  Google Scholar 

  50. Seixas A, Oliveira P, Termignoni C, Logullo C, Masuda A, da Silva VI. Rhipicephalus (Boophilus) microplus embryo proteins as target for tick vaccine. Vet Immunol Immunopathol. 2012;148:149–56.

    Article  CAS  PubMed  Google Scholar 

  51. Hussein HE, Johnson WC, Taus NS, Suarez CE, Scoles GA, Ueti MW. Silencing expression of the Rhipicephalus microplus vitellogenin receptor gene blocks Babesia bovis transmission and interferes with oocyte maturation. Parasit Vectors. 2019;12:7.

    Article  PubMed  PubMed Central  Google Scholar 

  52. Vechtova P, Fussy Z, Cegan R, Sterba J, Erhart J, Benes V, et al. Catalogue of stage-specific transcripts in Ixodes ricinus and their potential functions during the tick life-cycle. Parasit Vectors. 2020;13:311.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Huelsken J, Birchmeier W. New aspects of Wnt signaling pathways in higher vertebrates. Curr Opin Genet Dev. 2001;11:547–53.

    Article  CAS  PubMed  Google Scholar 

  54. Clevers H. Wnt/beta-catenin signaling in development and disease. Cell. 2006;127:469–80.

    Article  CAS  PubMed  Google Scholar 

  55. Artavanis-Tsakonas S, Rand MD, Lake RJ. Notch signaling: cell fate control and signal integration in development. Science. 1999;284:770–6.

    Article  CAS  PubMed  Google Scholar 

  56. Halder G, Johnson RL. Hippo signaling: growth control and beyond. Development. 2011;138:9–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Aikin RA, Ayers KL, Thérond PP. The role of kinases in the Hedgehog signalling pathway. EMBO Rep 2008; 9(4):330–6.

  58. Ingham PW, Nakano Y, Seger C. Mechanisms and functions of Hedgehog signalling across the metazoa. Nat Rev Genet. 2011;12:393–406.

    Article  CAS  PubMed  Google Scholar 

  59. Morita K, Chow KL, Ueno N. Regulation of body length and male tail ray pattern formation of Caenorhabditis elegans by a member of TGF-beta family. Development. 1999;126:1337–47.

    Article  CAS  PubMed  Google Scholar 

  60. Massagué J. How cells read TGF-beta signals. Nat Rev Mol Cell Biol. 2000;1:169–78.

    Article  PubMed  Google Scholar 

  61. Shi Y, Massagué J. Mechanisms of TGF-beta signaling from cell membrane to the nucleus. Cell. 2003;113:685–700.

    Article  CAS  PubMed  Google Scholar 

  62. Martín-Blanco E. p38 MAPK signalling cascades: ancient roles and new functions. BioEssays. 2000;22:637–45.

    Article  PubMed  Google Scholar 

  63. Flynn JL, Chan J. Immunology of tuberculosis. Annu Rev Immunol. 2001;19:93–129.

    Article  CAS  PubMed  Google Scholar 

  64. Eskelinen EL, Tanaka Y, Saftig P. At the acidic edge: emerging functions for lysosomal membrane proteins. Trends Cell Biol. 2003;13:137–45.

    Article  CAS  PubMed  Google Scholar 

  65. Rottensteiner H, Theodoulou FL. The ins and outs of peroxisomes: co-ordination of membrane transport and peroxisomal metabolism. Biochim Biophys Acta. 2006;1763:1527–40.

    Article  CAS  PubMed  Google Scholar 

  66. Moraes J, Galina A, Alvarenga PH, Rezende GL, Masuda A, da Silva Vaz Jr I, et al. Glucose metabolism during embryogenesis of the hard tick Boophilus microplus. Comp Biochem Physiol A Mol Integr Physiol. 2007;146:528–33.

    Article  PubMed  Google Scholar 

  67. Coelho A, Fraichard S, Le Goff G, Faure P, Artur Y, Ferveur JF, et al. Cytochrome P450-dependent metabolism of caffeine in Drosophila melanogaster. PLoS ONE. 2015;10:e0117328.

    Article  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We appreciate the contributions of Zhang Qian for her help in samples collection and data analysis.

Funding

This research was supported by the National Natural Sciences Foundation of China (No. 81871686).

Author information

Authors and Affiliations

  1. Collaborative Innovation Center for the Origin and Control of Emerging Infectious Diseases, Shandong First Medical University (Shandong Academy of Medical Sciences), Tai’an, China

    Zhang Ruiling, Liu Wenjuan, Zhang Kexin, Wang Xuejun & Zhang Zhong

  2. School of Basic Medical Sciences, Shandong First Medical University (Shandong Academy of Medical Sciences), Tai’an, China

    Zhang Ruiling, Liu Wenjuan, Zhang Kexin & Zhang Zhong

  3. Shandong Provincial Center for Disease Control and Prevention, Jinan, China

    Wang Xuejun

Authors
  1. Zhang Ruiling
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Liu Wenjuan
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Zhang Kexin
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Wang Xuejun
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Zhang Zhong
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

ZL, WJ, and ZZ conceived and designed the experiments, ZL, LJ, and ZX performed the experiments, analyzed the data ZL, and LJ wrote the manuscript. WJ and ZZ edited and reviewed the manuscript. All authors read and approved the manuscript.

Corresponding authors

Correspondence to Zhang Ruiling, Wang Xuejun or Zhang Zhong.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no conflict of interest.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1

. The partial mitochondrial genome (9585 bp) used to identify species and the BLAST results against NCBI.

Additional file 2: Figure S1

. Species distribution of Nr annotation results. The x-axis represents the name of species; the y-axis represents the number of unigenes annotated to different species.

Additional file 3: Figure S2

. GO functional classification of unigenes. The x-axis represents the GO term; the y-axis represents the number of genes classified into each GO term.

Additional file 4: Figure S3

. COG functional classification of unigenes. The x-axis represents the COG category; the y-axis represents the number of genes classified into each COG category.

Additional file 5: Figure S4

. KEGG enriched pathway analysis of genes between adjacent embryonic developmental groups. The x-axis represents the rich factor, the y-axis represents KEGG pathways. Rich factor refers to the ratio of the number of genes located in the KEGG pathway and the total number of genes in the KEGG pathway. The larger the rich factor, the greater the enrichment.

Additional file 6: Figure S5

. WGCNA of genes in different embryonic developmental stages. Hierarchical cluster tree showing co-expression modules identified by WGCNA. Each leaf on the tree represents one gene. All genes were grouped into 16 different color-coded modules according to their expression pattern.

Additional file 7: Figure S6

. Gene numbers of different modules.

Additional file 8: Figure S7

. Express pattern of lightcyan, darkred, magenta, lightyellow, orange, midnight, steelblue, royalblue, and white modules. The x-axis represents the sample; the y-axis represents the expression profile of the eigengene.

Additional file 9: Figure S8

. GO functional classification of the blue, darkgreen, brown, tan, darkgrey, and green modules. The x-axis represents GO terms; the y-axis represents the number of genes classified into each GO term.

Additional file 10: Table S1

. Primers used in qRT-PCR.

Additional file 11: Table S2

. Summary of output statistics from different samples.

Additional file 12: Table S3

. KEGG enriched pathways of all unigenes.

Additional file 13: Table S4

. Gene numbers classified into three GO catalogs. Table S5. Top 20 GO terms of different comparison groups.

Additional file 14: Table S6

. Enriched KEGG pathways of differentially expressed genes between adjacent embryonic developmental groups.

Additional file 15: Table S7

. The top ten hub genes of the blue, brown, darkgrey, darkgreen, green, and tan modules.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, 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 changes were made. 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/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Ruiling, Z., Wenjuan, L., Kexin, Z. et al. Developmental transcriptomics throughout the embryonic developmental process of Rhipicephalus turanicus reveals stage-specific gene expression profiles. Parasites Vectors 15, 89 (2022). https://doi.org/10.1186/s13071-022-05214-w

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13071-022-05214-w

Keywords

  • Tick
  • Embryonic development
  • Transcriptome
  • Stage-specific

Parasites & Vectors

ISSN: 1756-3305

Contact us