DNA Methyltransferases Contribute to Cold Tolerance in Ticks Dermacentor silvarum and Haemaphysalis longicornis (Acari: Ixodidae)

DNA methylation, mediated by DNA methyltransferases (Dnmts), is a typical epigenetic process that plays an important role in affecting organism acclimatization and adaptation to environmental changes. However, information about Dnmts and their associations with the cold tolerance of ticks remains meager. Hence, in the present study, the Dnmts in important vector ticks Dermacentor silvarum and Haemaphysalis longicornis were cloned and identified, and their functions in cold response were further explored. Results showed that the length of DsDnmt and DsDnmt1 in D. silvarum, and HlDnmt1 and HlDnmt in H. longicornis were 1,284, 549, 1,500, and 1,613 bp, respectively. Bioinformatics in protein analysis revealed that they were all unstable hydrophilic proteins and were mainly characterized with Dcm (DNA cytosine methyltransferase domain), Dnmt1-RFD (DNA methyltransferase replication foci domain), zf-CXXC (zinc finger-CXXC domain), and BAH (Bromo adjacent homology domain). The relative expression of these Dnmts was reduced after cold treatment for 3 days (P < 0.05), and increased with the extension of treatment. Western blot revealed that Dnmt1 decreased first and then increased significantly (P < 0.05) in both tick species, whereas other Dnmts fluctuated at varying degrees. RNA interference significantly silenced the genes Dnmts (P < 0.01), and mortality increased significantly (P < 0.05), when exposed to sub-lethal temperature, underscoring the important roles of Dnmts during the cold response of D. silvarum and H. longicornis. The above results lay the foundation for further understanding of the epigenetic regulation of DNA methylation in cold acclimatization and adaptation of ticks.


INTRODUCTION
Ticks are obligate blood-sucking ectoparasitic arthropods, which can spread a variety of pathogens and cause zoonotic diseases (1). They are considered as the most important vector of pathogens of domestic and wild animals, and only second to mosquitoes as vectors of human disease (2). With social and economic development and changes in the ecological environment, humans have more contact with ticks which has led to an increase in the incidence of tick-borne diseases (1). Ticks and tick-borne diseases have raised public health issues and posed a serious threat to animal husbandry development. Recently, ticks have been implicated as the vector of Songling virus (SGLV) (3) and Alongshan virus (ALSV), both associated with human febrile illness (4,5), and severe fever with thrombocytopenia syndrome virus (SFTSV), the causative agent of severe fever with thrombocytopenia syndrome (SFTS) (6).
The tick Dermacentor silvarum Olenev is widely distributed in Russia, Mongolia, and the north of China (7,8), and has been implicated as the vector and reservoir for the maintenance of tick-borne encephalitis virus through transstadial and transovarial transmission in China (9). Rickettsia raoultii, R. slovaca, and R. heilongjiangensis (spotted fever group rickettsiae), along with Ehrlichia chaffeensis (the causative agent of human monocyclic Ehrlichia), have been identified in D. silvarum (10,11). D. silvarum is also the primary vector of two babesiosiscausing parasites, Babesia equi and B. caballi; and the North Asian tick typhus pathogen known as R. sibirica (1). The tick Haemaphysalis longicornis Neuman, is one of the most predominant tick species in Asia countries and a serious pest of livestock in the Australasian and Western Pacific Regions (12,13). H. longicornis has expanded rapidly to many states in the USA since its first identification in New Jersey (14)(15)(16). It has been implicated as a vector of several pathogenic agents including the causative agent (Theileria sergenti and T. buffeli) of hemolytic anemia in wild and domestic cattle (17), causative agent (T. orientalis) of bovine theileriosis in Australia, and New Zealand (18), and causative agent (Babesia spp.) of livestock babesiosis (19). Moreover, it can transmit many spotted fever group rickettsiae (SFGR) such as R. japonica, R. heilongjiangensis, R. raoultii, and Candidatus Rickettsia tarasevichiae (12,(20)(21)(22). More recently, the severe fever with thrombocytopenia syndrome virus (SFTSV) was also identified in H. longicornis from South Korea, Japan, and China, which has caused great damages to livestock production and human health (6,23,24).
Many terrestrial arthropods are frequently faced with temperatures low and sustained enough to freeze their body fluids and therefore, have implemented some strategies to mitigate this risk (25,26). Some of these strategies include the modification of body fluid constituents to prevent ice crystallization (27), cryoprotective dehydration (28), freeze avoidance, behavioral avoidance, and freeze tolerance which involves the conversion of about 82% of body water into internal ice (26,29). Although freeze-tolerant arthropods utilize many potentially protective molecules, there is no obvious molecule that appears to be sufficient to support this coldtolerance strategy (26). Climate variables are important factors affecting tick distribution (30), which also influence tick survival, development, and activity during off-host periods (31,32). Hence, the ambient temperature, especially the low temperatures in winter, could play a significant role in determining the Abbreviations: Leu, Leucine; Ser, serine; Pro, proline; Arg, arginine; Asp, aspartic acid; Glu, glutamic acid; Lys, lysine; Gly, glycine; Ala, alanine; Trp, tryptophan; Cys, cysteine; Val, valine; Gln, glutamine; Phe, phenylalanine; His, histidine. host-seeking behavior of ticks. A complex physiological and biochemical process is involved in the response of ticks to low temperatures, including changes in the intracellular environment and substance synthesis, as well as the expression and regulation of some temperature-stress-related genes (33,34).
DNA methylation is an epigenetic modification that is heritable and reversible, which could play some roles in the mediation of rapid plastic responses of organisms to environmental stress (35,36). The methylation of DNA at the fifth position of cytosine (5 mC) and the sixth position of adenine (6 mA) influences gene expression by the regulation of the chromatin architecture and gene transcription (37,38). The capacity of organisms to adapt to environmental perturbations is increased by DNA methylation (39). This mechanism is mediated by a vital group of evolutionarily conserved enzymes known as DNA methyltransferases (Dnmts), which catalyze the methylation of DNA at cytosine residue to form 5-methylcytosine (40). Though mounting evidence has highlighted that variation in DNA methylation patterns could contribute to the acclimation and adaptation of organisms to different environmental stresses (41)(42)(43), the relationship between Dnmts and the cold tolerance of ticks remains poorly understood.
In the present study, the DNA methyltransferases in D. silvarum and H. longicornis were cloned and identified, and their functions in the cold tolerance of ticks were further explored via RNA interference. Western-blot methods were used to investigate the expression dynamics of DNA methyltransferases of two ticks under cold stress. The findings of this study provide insight into the regulatory role of DNA methylation in the cold response of these two tick species, and lay the foundation for further investigation on the epigenetic and molecular regulatory mechanisms underlying the cold tolerance of ticks.

Ticks
The unfed adult ticks of Dermacentor silvarum and Haemaphysalis longicornis were collected by flag-dragging from the vegetation of Xiaowutai National Nature Reserve, Zhangjiakou City, Hebei Province (114 • 47'−115 • 28'E, 39 • 50'-40 • 6'N). They were brought back to the laboratory and put on the ears of the New Zealand white rabbits for feeding, and cultured in an environmental chamber (26 ± 1 • C, relative humidity 75 ± 5%, 16 h of light: 8 h in darkness) during the non-parasitic period (44). The second-generation adult ticks that molted for two weeks were used in subsequent tests. All experiments were approved by the Animal Ethics Committee of Hebei Normal University (Protocol Number: IACUC-176810).

RNA Extraction and cDNA Synthesis
Total RNA was extracted from ticks (each group contains 10 D. silvarum or 20 H. longicornis) using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocol. RNA quantification was carried out using Nano Drop R ND-1000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), with a ratio of A260/A280 typically above 2.0, and the RNA quality was evaluated by 1% agarose gel electrophoresis. Then, 2.0 µg of each extracted RNA was used in a 20.0 µL reaction volume with an oligo (dT) 18 primer, and reverse transcription was performed to synthesize cDNA following the instructions provided with the TransScript R One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech Co., Ltd, Beijing, China). Primers (Table 1) were designed using DNAMAN (Lynnon Biosoft, San Ramon, CA, USA) and Primer Premier version 5.0 for Windows (Premier Biosoft International, Palo Alto, CA, USA), based on an alignment of several sequences obtained from the National Center of Bioinformatics Information (NCBI) website (http://www.ncbi. nlm.nih.Gov) and transcriptomic assembled sequences (Bio-Project PRJNA565692). The PCR conditions include an initial 2 min denaturation at 94 • C, 35 cycles of 30 s at 94 • C, 30 s at the melting temperature T m (58-60 • C) of the gene-specific primers, and 30 s at 72 • C, and then a final extension at 72 • C for 10 min carried out on an Applied Biosystems R Veriti R 96-Well Thermal Cycler (Life Technologies Ltd, Marsiling, Singapore). The amplified fragments were verified and separated on a 1% agarose gel. Bands with expected sizes were excised and purified by EG101-01 EasyPure Quick Gel Extraction Kit (TransGen, Beijing, China) according to the manufacturer's protocol. The purified products were sequenced and used for subsequent analysis.
To reconstruct the evolutionary relationships of the DNA methyltransferase genes, phylogenetic analysis was performed using Molecular Evolutionary Genetics Analysis (MEGA) software version 6.06 (47) based on the amino acid sequences retrieved from the NCBI protein database. Alignment of sequences was carried out with the ClustalW method, and the phylogenetic tree was obtained by the Neighbor-Joining method (48) with complete gap deletion and a bootstrap value of 100% in 1,000 replicates. Distances were computed using the Poisson correction method (49) and were in units of the number of amino acid substitutions per site.

Synthesis of dsRNA of Target Genes
The primers ( Table 1) for RNA interference of DNA methyltransferase genes were designed using DNAMAN (Lynnon Biosoft, San Ramon, CA, USA) and Primer Premier Version 5.0 for Windows (Premier Biosoft International, Palo Alto, CA, USA). Oligonucleotide primers comprising T7 promoter sequences at the 5 ′ end for in vitro transcription and synthesis of double-strand RNA (dsRNA) were used to amplify cDNA encoding the DNA methyltransferases. All the oligonucleotide primers used for this study were synthesized by Invitrogen TM (Beijing, China). In vitro dsRNA synthesis was carried out using the T7 RiboMAX TM Express RNAi System (Promega, Madison, USA). The transcription reaction contained 1 µL of cDNA as the transcription template, 7.0 µL of RNase-free water, and 2 µL of T7 enzyme mix in 10 µL Express 2x Buffer culminating to a final volume of 20.0 µL.
The reactions were carried out at 37 • C for 30 min, 70 • C for 10 min, and allowed to stand for 20 min at room temperature. After adding 2 µl of RQ1 RNase-Free DNase and 2 µL of diluted RNase A solution, the reaction was incubated at 37 • C for 30 min. Then, 4.4 µL Sodium Acetate and 110 µL 95% ethanol was added and incubated on ice for 5 min, followed by centrifugation at 16,000 g for 5 min at 4 • C. The precipitate was added 500 µL of 70% ethanol (pre-cooled at −20 • C), resuspended and centrifuged at 10,000 g for 5 min at 4 • C, and the precipitate was dried at room temperature for 15 min. Finally, the dsRNA pellet was resuspended in diethylpyrocarbonate (DEPC)-treated water.
The concentration and quality of dsRNA were determined by Nano Drop R ND-1000 Spectrophotometer and 1% agarose gel electrophoresis. The dsRNA was then stored at −80 • C until use.

RNA Interference
The procedure for the injection of dsRNA into 40 female ticks (20 D. silvarum and 20 H. longicornis) was performed according to (50). After surface-sterilization of the ticks, 2 µL of dsRNA solutions were injected into the hemocoel of unfed female ticks through the third and fourth coxa at a concentration of 1 µg/tick using a 10 µL Microliter TM Syringes (Hamilton, Nevada, USA). The control groups were injected with dsRNA-GFP. After injection, the ticks were placed in the environmental incubator (temperature 26 ± 1 • C, relative humidity 75 ± 5%, 16 h of light: 8 h in darkness) for 24 h. Each treatment was carried out in three replicates.
The total RNA was extracted from the ticks using TRIzol reagent, and their cDNA was synthesized as described previously. The mRNA expression of the target genes after knockdown with dsRNA was evaluated by quantitative real-time PCR (qRT-PCR), performed using an Mx3005P qPCR system (Agilent Technologies, Santa Clara, USA) using TransStart R Top Green qPCR SuperMix (TransGen Biotech) according to the manufacturer's instructions. The real-time PCR assays were prepared using a 96 well plate (Axygen, California, USA) with 10 µL of 2 × TransStart R Top Green qPCR SuperMix, 1.0  (40). The relative expressions of the DNA methyltransferase genes were calculated by a 2 − CT method (51).

Cold Tolerance of Ticks After RNAi
After confirmation of the knockdown of DNA methyltransferase genes, D. silvarum and H. longicornis ticks were exposed to a sub-lethal temperature of −22 and −14 • C, respectively, for 2 h in a low-temperature thermostatic bath (Meixiang Instrument Co. Ltd., Shanghai, China). The temperature of −22 • C was selected given that the discriminating temperatures (resulting in 20% survival) for the females and males were −21.7 and −22.6 • C, respectively (52). Then the ticks were removed from the thermostatic bath, and the mortality rate was calculated. Ticks were considered dead if they cannot coordinate their appendage after stimulated by CO 2 .

Relative Expression of DNA Methyltransferase Proteins Under Different Cold Treatments
The unfed adult D. silvarum or H. longicornis were divided into 6 groups with each group comprising 20 ticks. Three groups kept at the normal environment chamber (26 • C) served as control; whereas the other three groups (experimental group) were exposed to low temperature (4 • C) for 3, 6, and 9 days, respectively. After treatment, the ticks were grounded to powder in liquid nitrogen and transferred to a pre-cooled 1.5 mL Eppendorf (EP) tube containing 1 mL of sterile phosphatebuffered saline (PBS, pH 7.4) and 10 µL of phenylmethylsulfonyl fluoride (PMSF) protein inhibitor. Then, it was vortexed and centrifuged at 13,000, rpm for 10 min at 4 • C. The supernatant was transferred to a pre-cooled 1.5 mL EP tube. Protein concentration was measured using BCA (Bicinchoninic Acid) Protein Assay Kit (CWBIO, Haimen, Jiangsu, China), following the manufacturer's recommendations. The extracted protein samples were dissolved in 0.1 M Tris-HCl (pH 6.8) containing 2% SDS, 5% 2-mercaptoethanol, 10% glycerol, and 0.05 % bromophenol blue, followed by 2 min boiling. The SDS-PAGE was performed on 12% separating gels with 4% stacking gels containing 0.1% SDS using PAGE Gel Preparation Kit (Epizyme, USA). A 0.05 M Tris-glycine buffer (pH 8.3) containing 0.1% SDS was used to perform electrophoresis at a constant voltage of 120 V until the tracking dye reached the end of the gel. Subsequently, the electro-transfer onto polyvinylidene difluoride (PVDF) membrane was carried out at 22V for 25 min using a Trans-Blot SD apparatus (BIO-RAD, Hercules, CA, USA). The membrane was blocked for 2.5 h with 5% skim milk in 1xTBST (Tris-buffered saline with 0.1% Tween R 20 detergent) solution, washed three times with 1xTBST in 15 min, and then incubated overnight with a monoclonal anti-DNMT antibody (1: 2000) (GeneTex, TX, USA). After three washes with 1xTBST in 15 min, the membrane was incubated with horseradish peroxidase (HRP)-conjugated goat anti-mouse IgG for 2 h at room temperature, the signal was detected by ultra-sensitive enhanced chemiluminescent (ECL) substrate using SuperSignal TM West Femto Trial Kit (Thermo Scientific, Waltham, MA, USA), visualized and analyzed with an Image Lab Software-controlled Gel Doc TM XR+ System (BIO-RAD).

Statistical Analysis
The differences among different treatment groups were compared with one way ANOVA test, and Tukey's test was used for post hoc analysis. A significant difference is defined as P < 0.05. Data analysis was carried out using Graph Prism version 8.0.2 for Windows (GraphPad Software, Inc., San Diego, USA).

Cloning of DNA Methyltransferase Genes
The cDNA obtained from the total RNA of D. silvarum and H. longicornis were used as a template with specific primers for PCR amplification of the target DNA methyltransferase gene. The amplified products subjected to agarose gel electrophoresis showed single, clear, and bright bands (Supplementary Figure 1). The target fragments of D. silvarum DNA methyltransferase gene 1 were 493, 456, 552, and 539 bp, while the target fragment for DNA methyltransferase gene 2 was 549 bp. The target fragments of H. longicornis DNA methyltransferase gene 1 were 433, 851, and 526 bp, whereas the target fragments for gene 2 were 389, 517, 470, and 467 bp. After the blast, the two DNA methyltransferases of D. silvarum showed the highest sequence identity with the DNA cytosine methyltransferase (Dnmt) of Ixodes ricinus (65%) and DNA (cytosine-5)-methyltransferase (Dnmt1) of Ixodes scapularis (64%), whereas the two DNA methyltransferases of H. longicornis showed the highest sequence identity with the DNA methyltransferase 1 (Dnmt1) (81%) and DNA (cytosine-5)-methyltransferase PliMCI (Dnmt) (77%) of Ixodes scapularis. Thus, the DNA methyltransferase of D. silvarum and H. longicornis were named DsDnmt (1,284 bp), DsDnmt1 (549 bp), HlDnmt1 (1,500 bp), and HlDnmt (1,613 bp), respectively. The multiple sequence alignments indicate that the deduced amino acid sequences of DsDnmt, DsDnmt1, HlDnmt1, and HlDnmt showed high sequence identity to their corresponding sequences (Supplementary Figure 2).

Sequence Analysis
The DsDnmt had a relative molecular mass of 391957.00 Daltons with a base ratio of A+T at 5.56% and a base ratio of G+C at 54.44%. The predicted full length of its open reading frame (ORF) was 1,128 nucleotide (nt) ranging from 67 to 1,194 base pairs. DsDnmt1 had a relative molecular mass of 166573.00 Daltons with the base ratio of A+T <44.99% and G+C more than

Prediction of Amino Acid Hydrophobicity
The hydrophobic properties of the DNA methyltransferases were determined with the Hphob/Kyte&Doolittle model by a window size of 9 amino acids. The maximum hydrophobic value for DsDnmt was 2.078, located at the 50th and 121st amino acids which were Ala and Phe, respectively, whereas the minimum was −2.911 at the 111th amino acid (Gln). For DsDnmt1, the maximum hydrophobic value was 2.233, which was the 172nd amino acid (Ala), while the minimum was −2.489, the 92nd amino acid (His). The maximum hydrophobic value for HlDnmt1 was 1.722 at the 412th amino acid (Ala) and the minimum value was −3.844 at the 223rd amino acid (Arg). HlDnmt had 2.044 as its maximum hydrophobic value at the 445th amino acid (Leu) and −2.900 as its minimum hydrophobic value at the 338th amino acid (Arg). The total average hydrophobicity of DsDnmt, DsDnmt1, HlDnmt1, and HlDnmt proteins are −0.348, −0.239, −0.728, and −0.459, respectively, suggesting that the four DNA methyltransferases were hydrophilic proteins (Supplementary Figure 3).

Prediction of the Protein Physical and Chemical Properties
The physicochemical properties of DNA methyltransferase protein were predicted through the ProtParam Expasy website (http://web.expasy.org/protparam/) and DiANNA website (http://bioinformatics.bc.edu/clotelab/DiANNA/). All four putative proteins were predicted to have a nuclear localization ( Table 2). The instability index for all four studied proteins was higher than 40, indicating that they were unstable, whereas the predicted disulphide bonds of the proteins might help to maintain their molecular shape, stability, and substrate specificity ( Table 2).

Phylogenetic Relationship
The evolutionary relationship between the DNA methyltransferases of D. silvarum and other species as well as those of DNA methyltransferases of H. longicornis were constructed (Figure 1). Results revealed that the DNA methyltransferases of D. silvarum and H. longicornis clustered the arthropod species into their corresponding taxonomic classes. The putative DsDnmt and DsDnmt1 share the same node with the sequences of the closely related I. ricinus (QCB92136.1), I. scapularis (EEC18555.1), and I. ricinus (QCB92132.1), I. scapularis (EEC00472.1), respectively. Similarly, the putative HlDnmt1 and HlDnmt share the same node with the sequence of the closely related I. scapularis (XP029828587.1, XP029831274.1). Ticks belong to a relatively independent branch (arachnids) and are shown to be closely related to other arachnids than the insects. Vertebrates such as mammals, birds, and fishes were placed in an independent branch indicating that their DNA methyltransferases are less closely related to those of ticks than mollusks and other marine organisms.  (Figure 2). The open reading frames of DsDnmt1 and HlDnmt are probably incomplete and it is speculated that they lack other domains at the N-and C-terminal, respectively.

Protein Secondary Structure Prediction
The prediction of the secondary structure of the DNA methyltransferase proteins was carried out and the results showed that the four DNA methyltransferases have mixed secondary structures (Supplementary Table 1). α-helix and random coiled loop account for a larger proportion of the structure, whereas β-sheet strand accounts for a relatively small proportion, which provides an important structural basis for protein functional implementation and conformation. Furthermore, the DNA binding sites and protein binding sites of the DNA methyltransferase proteins were predicted (Supplementary Table 2).

Protein Tertiary Structure (3D Structure) Prediction
The tertiary structures of the four DNA methyltransferase proteins were predicted by Swiss-Model (http://swissmodel. expasy.org/). The reliability range of the global model quality estimation (GMQE) is from 0 to 1. The GMQE values of DsDnmt, DsDnmt1, HlDnmt1, and HlDnmt proteins were 0.65, 0.79, 0.10, and 0.74, respectively. The tertiary structure prediction models of DsDnmt, DsDnmt1, and HlDnmt have higher confidence, except for HlDnmt1, possibly due to the lower homology of their amino acid sequences in the software database. The predicted DsDnmt contains one Ca 2+ and one SAH (S-adenosyl-homocysteine) binding site, where SAH is the binding site of the methyl donor molecule SAM (S-adenosyl-L-methionine molecule); DsDnmt1 contains two Zn 2+ binding sites. HlDnmt1 contains one Ca 2+ binding site, whereas HlDnmt contains five Zn 2+ and one SAH binding site (Figure 3).

Relative Expression of DNA Methyltransferase Genes Under Cold Treatments
The relative expression levels of HlDnmt1, DsDnmt1, and DsDnmt decreased significantly after 3 days of cold treatment (P < 0.01), but increased with the extension of cold treatment time/days (Figures 4A-C). HlDnmt expression decreased significantly after 3 days of cold treatment (P < 0.05), but had a sharp significant increase (P < 0.01) after 9 days of cold treatment (Figure 4D).

The Efficiency of RNA Interference of DNA Methyltransferase
The relative expressions of the DNA methyltransferase gene were determined after RNAi. Results showed that the expression of the DNA methyltransferase genes was decreased significantly (P < 0.01) (Figure 5). Analysis of gene silencing efficiency showed that the average silencing efficiency of DsDnmt and DsDnmt1 genes were 74.79 and 73.78%, respectively, whereas that of HlDnmt1 and HlDnmt genes were 82.42 and 83.72%, respectively, indicating that the dsRNA injection of DNA methyltransferase gene had a significant silencing effect on the expression of the target gene.

Cold Tolerance of D. silvarum and H. longicornis After RNAi
After the interference of DNA methyltransferases, the changes in cold tolerance of D. silvarum and H. longicornis were analyzed, and the role of DNA methyltransferase in the cold response of ticks was determined. The mortality rates of both tick species increased significantly (P < 0.05) (Figure 6). The resultant mortality rate from the knockdown of DsDnmt, DsDnmt1, and HlDnmt was more than 60%, while the resultant mortality from the knockdown of HlDnmt1 was over 70%, indicating the important function of DNA methyltransferase in the cold tolerance of D. silvarum and H. longicornis.

Relative Expression of DNA Methyltransferase Proteins Under Different Cold Treatments
A total of 5 DNA methyltransferase-related proteins were selected (Dnmt1, Dnmt2, Dnmt3A, Dnmt3B, and Dnmt3L) as the target proteins which were determined depending on the commercially available specific antibody products. The results showed that the relative expression levels of various DNA methyltransferase proteins in the ticks increased or decreased to varying degrees after cold treatment. In D. silvarum, the protein expression of Dnmt1 and Dnmt2 changed significantly (P < 0.05). After cold treatment for 3 days, the protein expression was down-regulated, but with the extension of treatment, the expression was up-regulated. After cold treatment, the expression of Dnmt3A increased significantly (P < 0.05), and then gradually decreased. Dnmt3L showed a similar trend of change as Dnmt3A, but not to a significant level (Figure 7). In H. longicornis, the expression of Dnmt1 decreased after cold treatment for 3 days, and then gradually increased significantly (P < 0.05) after cold treatment for 9 days. Dnmt3B showed a similar trend as Dnmt1. The protein expression of Dnmt3L up-regulated initially after cold treatment before down-regulating with time, just as in D. silvarum ticks (Figure 8).

DISCUSSION
In the present study, four DNA methyltransferase genes were characterized from D. silvarum and H. longicornis, and phylogenetic analysis revealed that they were clustered with the arthropod species into their corresponding taxonomic classes. The putative DsDnmt and DsDnmt1 share the same node with the sequences of the closely related species (I. ricinus and I. scapularis). Similarly, the putative HlDnmt1 and HlDnmt share  DsDnmt1 (B), HlDnmt1 (C), HlDnmt (D) measured using qRT-PCR. The results are expressed as the means (n = 3) ± SEM. The asterisks above bars indicate significant difference based on one-way ANOVA followed by Turkey's test, *P < 0.05, ***P < 0.01. the same node with the sequences of closely related I. scapularis. Ticks belong to a relatively independent branch and are closely related to other arachnids than the insects. Vertebrates such as mammals, birds, and fishes were placed in an independent branch indicating that their DNA methyltransferases are less closely related to those of ticks than mollusks and other marine organisms. Overall, these results confirm that DNA methyltransferases of ticks are good candidates for phylogenetic analysis at higher levels. However, a previous study could not find a common ancestry between putative I. ricinus DNA methyltransferase (IrDAMT) and other arachnid species, but it was rather placed in a sister branch of the cluster of Maxillopoda, Branchiopoda, and insect (38). Additionally, HlDnmt results suggest that they could be less suitable for a detailed study of the   evolutionary relationship among arachnids/arthropods as none of its orthologs were found (apart from I. scapularis), presumably because of their non-existence in a publicly available database.
All the four putative proteins were predicted to have a nuclear localization, which is consistent with their hypothetical assignment to the methyltransferase group of enzymes. This was corroborated by the three putative I. ricinus methyltransferases-IrDNMT1, IrDNMT3, and IrDAMT (38). The instability index is a measure used to determine whether proteins will be stable in a test tube (53), and the index values above 40 are an indicator for unstable proteins, which is true for the four DNA methyltransferases in the present study. The conserved domain search revealed that putative DsDnmt and DsDnmt1 contain DNA cytosine methylase (Dcm) domain with multiple substrates-and DNA-binding sites. At the Dcm domain target sites (CCAGG and CCTGG), a methyl group is added to the 5th position of the cytosine base in the context of CpG dinucleotide by C-5 cytosine-specific DNA methylases (C5 Mtase); this plays a vital role in gene expression by either inhibiting the binding of transcription factors or the synthesis of methyl binding protein and related chromatin remodeling factors (54,55). HlDnmt1 contains the DMAP1 (DNA methyltransferase 1-associated protein) domain and SANT domain. DMAP1 possesses intrinsic transcriptional inhibitory activity and could play an important role in the development and maintenance of the integrity of methylation via interaction with Dnmt1 (56), whereas the SANT domain contains a putative DNAbinding site (57). Moreover, some interesting findings showed that the SANT domain was involved in the activities of ATPdependent chromatin remodeling enzymes, as well as in the histone acetyltransferase and deacetylase complexes (57)(58)(59). The putative HlDnmt has BAH domain, Dnmt1-RFD at the interval of 19-148 bp, and Zinc finger-CXXC (zf-CXXC) domain. BAH domain is responsible for linking DNA methylation, replication, and transcriptional regulation (37,60). The functions of the BAH domain in chromatin biology are versatile including FIGURE 8 | The relative protein expression of Dnmts in H. longicornis after cold treatment. The relative level of the expression of DNA methyltransferase proteins were calculated, and the expression in the control group served as the reference, which was set at 1. The results are expressed as the means (n = 3) ± SEM. The asterisks above bars indicate a significant difference based on one-way ANOVA followed by Turkey's test, *P < 0.05. taking part in the protein-protein interaction that mediates gene silencing, recognition of methylated histones, and direct interaction with nucleosome (61). Dnmt1-RFD is located at the longer N-terminal region of Dnmt1 which targets the replication foci (62). Specifically, it ensures that the Dnmt1 protein is targeted toward the replication foci in a noncatalytic fashion, which allows it to methylate the correct residues, and has been associated with the conversion of critical histone lysine moieties (63). Of note, Dnmt1-RFD binds to HDACs (histone deacetylases) and DMAP1 in a manner that can mediate transcriptional repression (63). The zf-CXXC domain contains eight conserved cysteine residues (present in the typical structure of Dnmts) which bind to Zn 2+ ions and non-methylated CpG dinucleotides (64). Although it was generally acknowledged that methylated CpGs (cytosine guanine dinucleotides) are specifically recognized by proteins that encode MBDs (methyl-CpG-binding domains) (65,66), it has also been proposed that non-methylated CpG dinucleotides could present a protein-binding site (64). This was made possible by the identification of non-methylated CGBP (CpG-binding protein) whose DNA-binding capacity is dependent on a cysteine-rich zf-CXXC domain (67).
Temperature is one of the most important factors that greatly influence the development, growth, geographical distribution, abundance, and physiology of organisms (68,69). Organisms respond to changes in environmental temperature by activating stress response mechanisms. Studies suggest that epigenetic enzymes play a vital role in arthropod cold hardiness (70). To the best of our knowledge, this is the first report of the response of DNA methyltransferases of D. silvarum and H. longicornis ticks to cold stress. Several studies with other species have demonstrated the relationship between DNA methylation, a vital environmentinduced epigenetic phenomenon, and variations in temperature (40,(71)(72)(73). For example, there was a significant change in the expression level of Dnmt3 in the Senegalese sole (fish) after shortterm exposure to embryonic temperatures (42). In the present study, the first 3 days of cold exposure were marked with a significant decrease in the expression levels of all the Dnmt genes but increased significantly with the prolongation of days of cold treatment. A similar trend was observed in the armyworm (Mythimna separata) where a decreased pattern of expression was observed after short exposures to various temperatures (73). A sharp significant increase in the expression levels of HlDnmt was observed after 9 days of cold exposure (4 • C). This was corroborated by the significant up-regulation of the expression levels of Bemisia tabaci DNA methyltransferase genes (BtDnmt3) after exposure to low temperatures ranging from 0 to 20 • C (74). A possible explanation could be that for the first few days of cold exposure, the ticks employ other physiological strategies such as cryoprotective dehydration (28), the modification of body fluid constituents (27), and the utilization of numerous potentially protective molecules in their body (26) before the full activation of the activities of the Dnmts as the cold stress prolongs. Additionally, the same pattern of expression was observed via western blot where the relative expression of various DNA methyltransferase proteins in the ticks decreased or increased to varying degrees after cold treatment. The protein expression level was down-regulated after 3 days of cold treatment but was up-regulated with the extension of treatment time for both the Dnmts of D. silvarum and H. longicornis. Similarly, a study observed that cold and subzero temperature exposure in Eurosta solidaginis (goldenrod gall fly) resulted in the upregulation of several DNA methyltransferases (70). The down-regulation and up-regulation due to temperature variations is a confirmation of the thermal plasticity of Dnmt expression which is reflected in the corresponding protein expression. Ticks are remarkably plastic and can inhabit diverse ecological niches, from the tropics to the polar regions of the planet (2,75). Like insects, ticks are ectotherms that are sensitive to external temperature variations (76). The cold hardiness of ticks that enables them to overwinter could be a function of epigenetic modifications such as DNA methylation that might occur rapidly in response to large-scale temperature changes (40).
Due to the vital role of DNA methyltransferases in different biological processes, the knockdown of DsDnmt, DsDnmt1, HlDnmt1, and HlDnmt was carried out to identify their functions by assessing the mortality rate after the exposure of D. silvarum and H. longicornis ticks to sublethal low temperatures of −22 and −14 • C, respectively, for 2 h. Interestingly, the mortality rates of the two tick species significantly increased which underscores the impact of DNA methyltransferases on the cold tolerance of ticks. This result was corroborated by the resultant significant mortality of armyworms (whose MsDnmt1 genes have been knocked down) after exposure to a low temperature (73). Previous studies have demonstrated the effect of extreme temperatures on arthropod development which include reduced fecundity, deformity, and mortality (77)(78)(79). It is noteworthy that ticks, which are widely distributed globally, can adapt to different climatic regions and they might possess unique adaptability mechanisms. One of these adaptability mechanisms could be mediated by DNA methyltransferases. DNA methylation, a vital epigenetic mechanism induced by the environment, can serve as a potential link between phenotypic variability and temperature variation, which is reprogrammed by Dnmts (40).

CONCLUSION
This study characterized DNA methyltransferases in D. silvarum and H. longicornis, and the results suggest that they contribute to the cold tolerance of ticks. The above findings highlighted the functional implications of DNA methyltransferases during cold stress of ticks and thus could advance our understanding of the epigenetic regulation of cold acclimation and adaptation of ticks.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.