Genome-wide identification of histone lysine methyltransferases and their implications in the epigenetic regulation of eggshell formation-related genes in a trematode parasite Clonorchis sinensis

In silico identification of SET-HKMT proteins

The genomic and proteomic databases of C. sinensis (Platyhelminthes: Trematoda: Opisthorchiida) in GenBank (https://www.ncbi.nlm.nih.gov/genbank/) were searched using the BLAST programs with the aa sequences of human SET-HKMT and DOT1 family members [7,22]. The sequences of C. sinensis proteins retrieved by the similarity searches were used as queries in a second series of BLAST tests against the C. sinensis databases. Homologous sequences were similarly isolated from the genomic databases of Schmidtea mediterranea (Platyhelminthes: Turbellaria: Tricladida) in the GenBank and Schistosoma mansoni (Platyhelminthes: Trematode: Strigeidida), Hymenolepis microstoma (Platyhelminthes: Cestoda: Hymenolepididae), Echinococcus granulosus (Platyhelminthes: Cestoda: Taeniidae), and Globodera pallida (Ecdysozoa: Nematoda: Rhabditida) in the GeneDB (https://www.genedb.org/) (E-value cutoff, 1×10⁻⁵). Functional domains including SET were annotated by analyzing the retrieved sequences using the Simple Modular Architecture Research Tool (SMART, http://smart.embl-heidelberg.de). The schematic diagrams of the preserved domains of each protein were constructed using the SMART program.

Phylogenetic analysis

The aa sequences of the SET domains from a total of 126 proteins were aligned and manually trimmed using the ClustalX (ver. 2.1) and GeneDoc programs, respectively. The alignment was applied in the construction of a maximum likelihood tree using the PhyML (ver. 3.1) [23] or MEGA X program [24]. The analytical parameters were as follows: Jones-Taylor-Thornton model for aa substitution, gamma distribution rates with invariant sites (4 gamma categories), and partial deletion of gaps/missing data (site coverage cutoff, 95.0%). The gamma rate and proportion of invariant sites were empirically estimated using the initial alignment. The statistical significance of internal branching nodes was evaluated using the non-parametric Shimodaira-Hasegawa-like approximate likelihood ratio test (SH-aLRT) or the standard bootstrapping test of 1,000 replicas. The resulting trees were displayed using the TreeView program.

Immunological detection of methylated histone derivatives

C. sinensis adults (42-day-old) were collected from experimentally infected rats, as described previously [25]. The worms (>30 worms) were minced with a razor blade and homogenized in Triton extraction buffer (TEB: PBS containing 0.5% (v/v) Triton X-100 and a protease inhibitor cocktail (complete; Roche Diagnostics GmbH, Mannheim, Germany)) using a Dounce tissue grinder with a loose pestle (Type A; Wheaton, Millville, NJ, USA). The homogenate was centrifuged at 10,000 g for 20 min at 4°C, and the supernatant was stored as the cytosolic protein extract. After washing with TEB, the pellet enriched with nuclei was resuspended in 0.2 N HCl and incubated overnight at 4°C. The solution was centrifuged at 20,000 g for 5 min at 4°C, resulting in the nuclear protein extract. The cytosolic and nuclear protein extracts were similarly prepared from HeLa cells and the liver tissue of a normal rat. Total protein in the protein extracts was quantified using the Pierce BCA assay kit (Thermo Fisher Scientific, Rockford, IL, USA). All the protein samples (10 μg) were electrophoresed on 15.0% SDS-PAGE gels under reducing conditions and transferred onto nitrocellulose membranes (Schleicher & Schuell Bioscience, Dassel, Germany). The membranes were reacted with antibodies specific to histone H3 or its methylated derivatives (Supplementary Table S1) and then, with horseradish peroxidase (HRP)-conjugated goat anti-rabbit or anti-mouse IgG antibodies (Bethyl Laboratories, Montgomery, TX, USA). Positive signals were visualized using an enhanced chemiluminescence detection system (GE Healthcare, Pittsburgh, PA, USA). Images were captured using a LAS-400 Lumino image analyzer (Fujifilm, Tokyo, Japan).

The whole bodies of C. sinensis adults were fixed in paraffin blocks and sliced onto electrostatically charged glass slides (Superfrost Plus, Manzel-Glaser, Germany). After deparaffinization and rehydration according to the manufacturer’s instructions (DAKO, Carpinteria, CA, USA), the slides were soaked in a 3.0% H₂O₂ solution for 5 min and then in TBS/T buffer with 3.0% BSA for 1 h. The sections were incubated overnight at 4°C in the BSA-TBS/T buffer with each of the antibodies listed in Supplementary Table S1 (1:200 dilution). The primary antibodies bound to the worm sections were detected using HRP-conjugated secondary antibodies (1:1,000 dilution; Cappel, West Chester, PA, USA) and the HIGHDEF blue chromogen substrate (Enzo Life Sciences, Farmingdale, NY, USA). A mixture of pre-immune rabbit serum (n=3, 1:200 dilution) was used as a negative control for the immunohistochemical staining.

Temporal expression profiles of C. sinensis genes

Total RNAs were extracted from C. sinensis worms at various stages (>30 worms/stage) using QIAzol solution and the RNeasy Mini kit (Qiagen, Hilden, Germany). After removing residual genomic DNAs using RNase-free DNase (New England Biolabs, Ipswich, MA), Complementary DNAs (cDNAs) were synthesized from the RNAs using the iScript cDNA synthesis kit (Bio-Rad, Munich, Germany). The numbers of epigenetic tool gene and reproduction-related gene transcripts were examined via quantitative real-time PCR (RT-qPCR) using the cDNA samples and gene-specific primers (Supplementary Table S2). qPCR was conducted using the SYBR Green Master Mix and the CFX96 detection system according to the manufacturer’s instructions (Bio-Rad Laboratories, Hercules, CA, USA). The presence of a single-target amplicon was confirmed by analyzing melting curves. Each reaction was conducted in triplicate, and the relative transcription levels against a reference β-actin gene (DF143505) were presented as the mean±SD (ΔCq calculation method) [26]. A Student’s t test was used for statistical analyses, and a probability level of P<0.05 was considered statistically significant.

Methyltransferase inhibitors and gene expressions

Live C. sinensis worms (42- and 140-day-old, 30 worms/group) were incubated in RPMI-1640 media (phenol red- and serum-free, pH 7.2) supplemented with 5 μm SGI-1027 (Selleck, Houston, TX, USA) or 0.5 μm 3-deazaneplanocin A (DZNep; Selleck) at 37°C in a 5.0%-CO₂ incubator, by refreshing the media daily. After a 5-day incubation, the total RNAs were extracted from the worms for gene expression analysis using RT-qPCR, as described above. The responsive gene expression levels in the experimental groups were calculated as fold changes against those in a control group, which was not treated with the chemicals (ΔΔCq calculation method) [26].

Identification of SET-HKMTs in C. sinensis and other helminths

Based on previous reports [7,22], several human SET-HKMTs were selected to represent each of the diverse HKMT families. The aa sequences of these human proteins were used as queries in BLAST searches against the C. sinensis proteomic database in GenBank (BioProject No. PRJNA563175). The GenBank database was further surveyed for isolating Clonorchis proteins named HKMT or HMT by the original authors. The retrieved aa sequences were compared to one another to remove redundant entries (cutoff, 95.0%), and a total of 32 proteins were finally identified as the non-allelic HKMT paralogs of C. sinensis (Table 1). Information on orthologous proteins was similarly obtained from other trematode (24 proteins in S. mansoni), cestode (28 proteins in E. granulosus and 41 proteins in H. microstoma), turbellarian (7 proteins in S. mediterranea), and nematode (20 proteins in G. pallida) species, for which genome sequences had been drafted, using the human and C. sinensis SET-HKMT sequences as queries (Supplementary Tables S3-S7).

Table 1

SET-HKMT family proteins idnetified in the proteomic database of Clonorchis sinensis

Structural characters of the helminth SET-HKMTs

The Clonorchis SET-HKMT proteins were assigned to multiple SET-HKMT families, including SUV39, SET1, SET2, and SMYD (Table 1). Screening of the C. sinensis transcriptomic database in GenBank with the SET-HKMT sequences demonstrated that the majority of C. sinensis proteins, except for GAA50403, GAA52264, and GAA50855, are expressed at least at the mRNA level. The SMART program recognized aa stretches preserved in various functional domains such as SET, plant homeodomain (PHD), and Set2-Rpb1 interaction (SRI) in their primary structures, although the main SET domain could not be defined in 11 sequences (Table 1). The highly polymorphic lengths (81–3,518 residues) and lack of the SET domain indicated that at least some of protein sequences had been incompletely predicted.

The Clonorchis SET domain sequences were aligned with one another, and their aa conservation patterns were examined by referring to previous reports [7,27]. The SET domain sequences were highly diverged, and the similarity values ranged from 2.0% (GAA49834 vs GAA52155) to 67.1% (GAA32467 vs GAA39455) (average±SD, 18.2%±0.10). Short conserved aa stretches, which are involved in the binding of SAM and/or histone Lys [7,27], were readily detected in the SET-N and SET-C domains of these C. sinensis proteins (Motifs I – IV in Fig. 1). However, critical amino acids involved in the formation of H-bonds/van der Waals interactions with the SAM molecule (Gly, Asn, and His indicated by open arrowheads) or the ε-amino group of histone Lys (Tyr indicated by solid arrowheads) were substituted by other amino acids in some of the C. sinensis proteins. Similar aa conservation patterns were observed in the SET domain sequences identified in other helminths (data not shown).

Fig. 1

Comparison of SET domain sequences conserved in the Clonorchis sinensis histone lysine methyltransferase (CsHKMT) proteins. The amino acid (aa) sequences of CsHKMT SET domains, which were defined in the full CsHKMT sequences using the SMART program, were aligned against one another for comparison of their primary structures. The degrees of the aa conservations are highlighted in different shades of grey, and gaps represented by dots were introduced in the alignment to increase the similarity values. The short aa signatures of the SET domain involved in the binding of S-adenosyl methionine (SAM) and histone Lys are presented above the alignment (Motifs I–IV). The solid triangle indicates the invariant Tyr forming H-bond/van der Waals interaction with the ε-amino group of histone Lys. Open arrowheads represent the positions of Gly, Asn, and His participating in the formation of the non-covalent bonds with SAM molecule.

Phylogeny of the helminth SET-HKMTs

The genomes of the helminth species encoded diverse SET-HKMTs, although some entries were not completely predicted (entries marked with asterisks in Table 1). The evolutionary relationships of these SET-HKMT orthologs were examined via phylogenetic analysis based on their SET domain sequences (21, 16, 18, 28, 4, 14, and 25 proteins from C. sinensis, S. mansoni, E. granulosus, H. microstoma, S. mediterranea, G. pallida, and human, respectively). As shown in Fig. 2, helminth proteins were clustered into tight clades together with human proteins, consistent with their respective names, whereas some clades, such as NSD2, SMYD4, SMYD5, and SETD4 comprised only helminth proteins. In addition, the genes encoding EHMT1 and SMYD4 had been specifically multiplied in G. pallida (solid arrowhead in Fig. 2) and H. microstoma (open arrowhead), respectively. Considering the fact that a considerable number of SET-HKMTs lacking the SET domain were excluded in the phylogenetic analysis, the major SET-HKMT clades were likely to have diverged early in bilaterian evolution, at least prior to the divergence of deuterostomes and protostomes.

Fig. 2

Phylogeny of Clonorchis sinensis histone lysine methyltransferases (HKMTs) and their platyhelminth and nematode orthologs. The amino acid sequences of HKMT SET domains identified in trematodes (C. sinensis [red] and Schistosoma mansoni [orange]), cestodes (Echinococcus granulosus [blue] and Hymenolepis microstoma [purple]), turbellarian (Schmidtea mediterranea [brown]), nematode (Globodera pallida [green]), and human (black) were aligned and used in the phylogenetic analysis. The unrooted maximum likelihood tree was constructed using PhyML. Branch support values obtained using the Shimodaira-Hasegawa-like approximate likelihood ratio test are indicated at the corresponding branching nodes. Proteins, coding genes of which appeared to be expanded specifically in G. pallida and H. microstoma, are indicate with solid and open arrowheads, respectively. The solid circles indicate C. sinensis proteins selected for the gene expression analyses.

Retrieving DOT1 protein sequences from the helminth proteomes

BLASTp searches using the aa sequence of human DOT1 (1,537 aa; GenBank accession number NP_115871.1) retrieved a single non-redundant protein from each of the trematode and cestode proteomic databases in the GenBank (E-value <2×10⁻¹¹). Similar to SET-HKMTs, the lengths of these helminth proteins also varied ranging from 218 (Schistocephalus solidus DOT1; VDL97768.1) to 1,756 aa (E. granulosus DOT1; XP_024351501.1), the majority of which were likely to be partial. In the C. sinensis database, 2 proteins with 1,243 (RJW63735.1) and 1,704 (GAA52952.1) overlapping aa sequences significantly matched to the human DOT1 protein (E-values, 1×10⁻¹⁰⁹ and 3×10⁻¹⁰⁵, respectively). No homologous protein was retrieved from the current S. mediterranea proteome. A similarity search of the entire genomic sequences of the turbellaria using tBLASTn algorithm isolated several short DNA segments, of which theoretically translated proteins were matched to the human protein (query coverage <14.0%). Given the low coverage and/or the locations of matched sequences in the 3′-terminus of respective contigs as well as the low SET-HKMT protein number, the retrieval failure occurred probably due to the incomplete assembly of the flat worm genome.

Comparative analysis of the helminth DOT1s

Although partial, the catalytic DOT1 domain (PFAM id. PF08123) was readily recognized in the primary structures of these helminth DOT1 homologs (Fig. 3A). Some of these DOT1 motif sequences were aligned to that of a previously characterized human DOT1 protein [28]. As shown in Fig. 3B, the full catalytic core sequences were relatively well-preserved among them and exhibited identity values ranging from 36.0% (EgDOT1 vs. BxDOT1) to 100.0% (CsDOT1 vs OvDOT1) (average±standard deviation, 62.9±0.14%). The aa conservation was much more prominent in the segmental motifs constituting the SAM-binding pocket (Motifs I, I′, II, D1, and D2). The aa residues implicated in either the direct H-bond interaction with the cofactor SAM (solid circles in Fig. 3B) or the constitution of an access channel for H3K79 (open circles) were tightly conserved between the helminth and human orthologs. Identical amino acids may be involved in the formation of a channel leading to the inside of the SAM-binding pocket (solid arrowheads) and the regulation of the Lys-binding channel gate (open arrowheads), as described in the human DOT1 protein [28]. In addition to the DOT1 domain, several SCOP domains, including d1quua1 (2 central spectrin-like repeats of the σ-actin family), d1cuna2 (repeats of the σ-spectrin family), d1fcha (a domain found in the peroxisomal targeting signal 1 receptor), and d1eq1a (a domain found in the apolipophorin-III family), were occasionally detected in the primary structures of some DOT1 homologs. The human DOT1 protein further preserves the SMART AT_hook domain (id. SM00384), which is a DNA-binding motif with a preference for A/T-rich regions (Fig. 3A).

Fig. 3

Structural characterization of Clonorchis sinensis DOT1 (CsDOT1) proteins. (A) The functional domains conserved in CsDOT1 were predicted using the SMART program and compared with those of its orthologs. The domain legends are shown at the bottom. (B) The SET domain sequences of DOT1 homologs shown in panel A were aligned and the degrees of amino acid (aa) conservation were differentially highlighted in the alignment. The secondary structure of the SET domain, which has been determined in human DOT1, is presented at the bottom of each alignment (rod, α-helix; arrow, β-sheet). The segmental motifs constituting the S-adenosyl methionine (SAM)-binding pocket (Motifs I, I′, II, D1, and D2) are shown on the top. Solid and open circles represent aa residues involved in the binding of SAM and Lys79 of histone H3, respectively. Amino acids involved in the formation of a channel leading to the inside of the SAM binding pocket and the regulation of the Lys-binding channel gate are denoted with solid and open arrowheads, respectively. The numerals in the alignment are the positions of corresponding residues in human DOT1 (NP_115871.1).

The aa sequences of platyhelminth DOT1 proteins were aligned with those of nematode and deuterostomian orthologs, as well as cnidarian proteins (Nematostella venctensis and Hydra vulgaris DOT1), which were selected to comprise an outgroup. Using the sequence alignment, a maximum likelihood tree was constructed to examine their evolutionary relationships. In the tree rooted with the cnidarian proteins, DOT1 homologs showed a clustering pattern corresponding to the phylogeny of their donor organisms (Fig. 4).

Fig. 4

A phylogenetic tree of platyhelminth and nematode DOT1. The amino acid sequences of the DOT1 domain were isolated from some representative platyhelminth and nematode species including C. sinensis. The sequences were used in the phylogenetic analysis, together with those of deuterostomian and cnidarian orthologs, using the PhyML program. The tree was rooted with the cnidarian DOT1s. Numerals at branching nodes indicate the bootstrap values (1,000 replicas).

Methylated histones in C. sinensis adults

The aa sequences of Clonorchis histones H3 and H4 were isolated from GenBank database by BLAST searches using the sequences of human histones H3.1 (NP_003522.1) and H4 (NP_003486.1). The aa sequences were highly conserved among these fluke and human homologs (identity values >98.0%), except for that of the C. sinensis histone H3.1 (KAG5455013.1, 69.2%) (Supplementary Fig. S1). Consistent with the tight aa conservation, pan-antibodies against histone H3 and its methylated derivatives were highly reactive with the nuclear extracts of C. sinensis adults (Supplementary Table S1; Supplementary Fig. S2). Using these antibodies, the histological distribution of the methylated histones was examined across the whole body of the fluke. The anti-histone H3 rabbit antibody stained the nuclei in the parenchymal region and sexual organs of C. sinensis adults, as well as those in intrauterine eggs. Antibodies against H3K4me3, H3K27me3, and H3K79me3 showed staining patterns similar to that of the anti-H3 antibody, whereas the relative intensity of the positive signal differed slightly among the antibodies and tissues/organs. However, the specific targets of anti-H3K36me3 and H4K20me3 antibodies were largely detected in the vitellaria and intrauterine eggs. Proteins in the nuclei of ovaries, testes, and seminal receptacle reacted minimally with these antibodies. Finally, a weak positive signal was observed in the intrauterine eggs but not in the other tissues/organs, in the immunohistochemical staining images obtained with anti-H3K9me3 antibody (Fig. 5).

Fig. 5

Immunological detection of methylated histones in Clonorchis sinensis adults. The whole-body sections of C. sinensis adults (42-day-old) were reacted with antibodies specific to histone H3 and its methylated derivatives, as indicated on the top, and the positive reactions were enzymatically stained with the HIGHDEF blue chromogen. Rabbit normal serum was used in the control staining reaction. Ov, ovary; Sr, seminal receptacle; Tg, tegument; Ts, testis; Vt, vitelline follicle.

Temporal profiling of SET-HKMT gene expressions

Six Clonorchis SET-HKMT genes (CsSET-HKMT: CsSETD4 for GAA51355, CsPRDM9 for GAA52155, CsEZH for GG55462, CsNSD1/2 for GAA51760, CsSUV420H for GAA53884, and CsSETD8 for GAA4871) were arbitrarily selected based on their phylogenetic positions (solid dots in Fig. 2), and their temporal expression profiles were examined in C. sinensis via the RT-qPCR analysis. The relative expression levels of these genes gradually increased during the development/maturation of the worm from metacercaria to 15-day-old worm, except for CsPRDM9 and CsSUV420H, expressions of which were reduced in 6-day-old juveniles compared to those in metaceracariae, and then significantly decreased in 25-day-old worms (P<0.001; Fig. 6). The relative expression of CsNSD1/2 and CsSETD8 was further reduced in 84-day-old worms, similar to other functional genes such as C. sinensis tyrosinases (CsTyrs) and myoglobins (CsMbs) [25,29]. However, the transcriptional activity of the other 3 genes was significantly restored in the later stages of C. sinensis adults (P values for CsSETD4, CsPRDM9, and CsSUV420H<0.01).

Fig. 6

Temporal expression profiles of Clonorchis sinensis histone lysine methyltransferase (CsHKMT) genes. Expression profiles of the C. sinensis genes were determined using quantitative reverse transcription PCR, in association with the development/maturation of the liver fluke from metacercaria to 84-day-old worm. The calculations are based on independent technical triplicates (n=3, mean±standard deviation). The significance of the expression changes at each stage was statistically estimated by comparing the gene expression levels with those of its direct previous stage using the Student’s t test. **P<0.01, ***P<0.001.

Effects of HKMT inhibitors on functional gene expression

In our previous studies, the temporal expressions of C. sinensis DNA methyltransferase (CsDNMT2) and methyl-CpG binding protein (CsMBD1 and 2) genes showed patterns similar to that of CsEZH (Supplementary Fig. S3). Taken together with the present results, these expression patterns suggest that epigenetic gene regulation via DNA and/or histone methylation increases in C. sinensis with age. To test this hypothesis, the changes in functional gene expressions (CsSET-HKMTs, CsMBDs, CsDNMT2, CsTyrs, and CsVits) were examined in 42- and 140-day-old C. sinensis worms treated with the methyl transferase inhibitors SGI-1027 and DZNep for 5 days. As shown in Fig. 7, the transcriptional activities of these C. sinensis genes were significantly affected in the 140-day-old worms (P values<0.05). Overall, SGI-1027 induced these functional genes, but DZNep down-regulated them in the older worms, whereas the degrees of the responsive changes differed substantially depending on the chemicals and genes (Fig. 7). In 42-day-old worms, both chemicals strongly activated the CsMBD2 and CsDNMT2 genes. However, CsTyr and CsVit expressions were not significantly affected by SGI-1027 or DZNep in the younger worms (Fig. 8). Interestingly, the individual genes encoding SET-HKMTs responded differently to these inhibitors in both younger and older worms (Figs. 7, 8).

Fig. 7

Effect of methyltransferase inhibitors on functional gene expressions in 140-day-old Clonorchis sinensis adults. Live C. sinensis worms were incubated in the RPMI-1640 medium supplemented with SGI-1027 (5 μm) and/or 3-deazaneplanocin A (DZNep, 0.5 μm) in a 5.0%-CO₂ incubator. After a 24 h incubation, total RNAs were extracted from the worms and used in quantitative reverse transcription PCR with primer sets specific to each of the C. sinensis genes. (A) HKMT genes; (B) methyl-CPG; (C) DNMT genes; (D) Tyrosinase genes; (E) eggshell precursor protin genes. Fold increases in the globin gene expression in each of the experimental groups were calculated against the media-only control group (n=3, mean±standard deviation). *P<0.05, **P<0.01, ***P<0.001.

Fig. 8

Effect of methyltransferase inhibitors on functional gene expressions in 42-day-old Clonorchis sinensis adults. Live C. sinensis worms were incubated in the RPMI-1640 medium supplemented with SGI-1027 (5 μm) and/or 3-deazaneplanocin A (DZNep, 0.5 μm) in a 5.0% CO₂ incubator. After a 24 h incubation, total RNAs were extracted from the worms and quantitative reverse transcription PCRs were performed using primer sets specific to each of the C. sinensis genes. (A) HKMT genes; (B) methyl-CPG; (C) DNMT genes; (D) Tyrosinase genes; (E) eggshell precursor protin genes. Fold increases in the globin gene expression in each of the experimental groups were calculated against the media-only control group (n=3, mean±standard deviation). *P<0.05, **P<0.01. ***P<0.001.