Metacercariae of M. orientalis were isolated from infected Pseudorasbora parva from the Wuyuer river basin (47.53°N, 124.44°E), Heilongjiang Province, China. Sheldrakes were orally infection with 100 metacercariae isolated from P. parva and then euthanized 6 weeks later. Adult M. orientalis were obtained from the liver and gallbladder of the ducks. These organs were thoroughly washed in sterile saline solution and frozen in liquid nitrogen until use. The metacercariae and adults were identified to the species level according to existing keys and descriptions (Sohn, 2009), immediately frozen in liquid nitrogen, and stored at -80°C until use.

Total RNA was extracted using TRIzol reagent (Invitrogen Life Technologies, Carlsbad, CA, USA) according to the manufacturer’s protocol in three biological replicates of each M. orientalis development stage (pool of adults comprising of n = 50 and pool of metacercariae comprising n = 2000). Total RNA of M. orientalis adults or metacercariae was stored at -80°C until use. Library construction was performed according to the Illumina sample preparation for RNA-sequencing (RNA-seq) protocol. Brifely, The Oligo (dT) was used to isolate poly (A) mRNA from total RNA. Then the mRNA is fragmented into short fragments by mixing with fragmentation buffer. The cDNA was synthesized using the mRNA fragments as templates. Short fragments were purified and dissolved with EB buffer for end reparation and then connected with adapters. The suitable fragments are selected for the PCR amplification as templates. During the quality control steps, Agilent 2100 Bioanaylzer and Qubit^® RNA Assay Kit in Qubit^® 2.0 Flurometer (Life Technologies, CA, USA) were used in quantification and qualification. Sequencing of the library preparations was performed by an Illumina Hiseq X Ten platform to obtain paired-end reads.

Raw reads were subjected to quality control to obtain clean reads by removing reads with adaptors, reads containing > 10% ‘N’ residues, and low-quality reads containing > 50% bases. Clean reads were assembled into unigenes base on the default settings of the Trinity program (Grabherr et al., 2011). Unigene sequences were aligned with the NCBI non-redundant nucleotide (NT) database (Liu et al., 2016) by BLASTn (Zhang et al., 2017), and aligned with the NCBI non-redundant protein (NR) database (Liu et al., 2016), Swiss protein (Swiss-Prot) database (Liu et al., 2016), Cluster of Orthologous Groups of proteins (COG) database (Liu et al., 2016) and Kyoto Encyclopedia of genes and genomes (KEGG) database (Liu et al., 2016) by BLASTx (Zhang et al., 2017) to assign the predicted function. Hmmscan version 3.3.2 (Finn et al., 2011) was employed to match the established HMM model of protein structure domain among the Pfam database (Zhang et al., 2017). ESTScan version 3.0.3 (Iseli et al., 1999) was employed to predict protein coding sequences (CDS) with default setting. Blast2GO (Götz et al., 2008) was employed to classify unigenes to Interpro and Gene Ontology (GO) terms including molecular function, biological processes, and cellular components (Conesa et al., 2005) and analyzed the distribution of M. orientalis gene functions at the macro-level. The clean reads were deposited in the Sequence Read Archive database of NCBI (accession no. PRJNA474572), with sra run accessions numbers SRR7410653 and SRR7410652 for adult and metacercaria M. orientalis, respectively. The assembled cDNA sequences were deposited in the Transcriptome Shotgun Assembly (TSA) database of GenBank (accession no. GGVK00000000).

Unigene expression was calculated based on the Fragments Per kb per Million reads (FPKM) method (Mortazavi et al., 2008). The FPKM values were used directly used to compare the differences in gene expression levels between the two developmental stages. The Benjamini-Hochberg procedure is used to perform multiple corrections to p-values and generate false discovery rate (FDR) values. Differentially expressed genes (DEGs) were identified with an adjusted FDR < 0.005 found by DESeq2 version 1.34.0 (Love et al., 2014).

Partial of total RNA same as transcriptome sequencing were used for quantitative real-time PCR (qRT-PCR) validation. Primers designed according to the Illumina sequencing data are listed in Table S1 . CDNA was synthesized from total RNA using the reverse transcription kit (Takara, Dalian, China) following the manufacturer’s instructions. Thermocycling conditions were: 40 cycles each with 95°C for 10 s for denaturation, 60°C for 20 s for annealing, and 72°C for 30 s for extension, performed in StepOne Plus Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) with SYBR Green Pre-mix Ex Taq (Takara, Dalian, China) in triplicate. Relative gene expression was calculated using the 2^-ΔΔCt method with β-actin (GenBank no. EU109284) as the internal control. The correlation coefficients between the transcriptome and qRT-PCR values were calculated.

Each M. orientalis development stage of total proteins (pool of adults comprising of n = 50 and pool of metacercariae comprising n = 2000) were extracted using protein lysis buffer (7 M urea, pH 8.0) in three biological replicates, and were lysed by sonication on ice (2/3 s, 5 min) using a high-intensity ultrasonic processor (Scientz Biotechnology Co. LTD, Ningbo, China).

The lysate was centrifuged at 20,000 × g for 20 min at 4°C in order to remove debris. After centrifugation, the supernatant was treated with 10 mM dithiothreitol for 60 min at 37°C. Then, the samples were alkylated with 55 mM iodoacetamide, protected from light for 45 min at room temperature. The concentration of protein was quantified using the BCA protein assay kit (Pierce, Rockford, IL, USA). For each sample approximately 10 μg of protein was subjected to 12.5% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) to assess protein integrity.

The protein samples for 50 μg were diluted with 30 mM HEPES until the concentration of urea becomes < 2 M. Trypsin was added into each sample at an enzyme to protein ratio of 1:50 and the samples were further digested overnight at 37°C. Enzymatic digestion was terminated by adding 0.5% (v/v) formic acid. Finally, peptides of each sample were desalted and concentrated using Sep-Pack C₁₈ Cartridges (Waters, Worcester, MA, USA).

All samples filtering experiment were separated by an HPLC system (Easy-nanoLC, Thermo Scientific, Chelmsford, MA, USA) connected to an orbitrap fusion mass spectrometer (Thermo Scientific, Chelmsford, MA, USA). Peptides were resuspended with phase A [2% acetonitrile (ACN), 0.1% formic acid (FA)] and centrifuged at 12,000× g for 10 min. The supernatant was loaded on the trap column to be enriched and desalted. Then, the peptides were separated at a flow rate of 300 nL/min on a 15 cm analytical column (Beijing Qinglian Biotech, China, 150 μm ×150 mm, 100A˚, 1.9 μm) connected to the trap column. The linear gradient of LC was set at 3% buffer B (95% ACN, 0.1% FA) (from 0 to 5 min), 8–28% buffer B (from 5 to 107 min), and finally a hold at 28–80% buffer B (from 107 to 120 min).

Peptides were ionized by a nano-electrospray ion source and then identified by the orbitrap fusion mass spectrometer (Thermal Scientific, Chelmsford, MA, USA) in the mode of DDA (data-dependent acquisition). The scan of first-grade MS ranged from 350 to 1550 m/z at a resolution of 120,000 and an automatic gain control (AGC) target of 2×10⁶. The scan of second-grade MS was initiated as 100 m/z at a resolution of 30,000 with a dynamic exclusive time of 30 ms and an AGC target of 5 × 10⁴. The mode of second-grade MS spectra was high-energy collisional dissociation, and inject ions for all available parallelizable time.

Raw mass spectra were searched against M. orientalis transcriptome database in present study, and protein identification was performed using MaxQuant version 1.5.3.30 (Wiśniewski et al., 2009). The search parameters were set as follows: first and main search peptide tolerances of 20 ppm and 6 ppm, respectively, precursor ion mass tolerances of 20 ppm, a maximum of two missed trypsin cleavage sites, fixed cysteine carbamidomethylation, and variable methionine oxidation. Then, the acetylated sites were filtrated at the level of site decoy fraction ≤ 1% to obtain the significant modification. The p-value for identification and quantification of proteins was set as p ≤ 0.05 and acetylated proteins with a fold-change of two were deemed as differentially expressed proteins (DEPs). Protein annotation used BLASTP on the UniProt database (Wang et al., 2011) with default parameters. GO enrichment analysis was used to determine whether the identified proteins were enriched in certain functional groups, as compared with the uniprot Clonorchis sinensis dataset (Wang et al., 2011), and Fisher’s exact test was used for the analysis. Next, the identified proteins were blasted against the kyoto encyclopedia of genes and genomes (KEGG) (Wang et al., 2011) for orthology identification of the corresponding genes, and subsequently mapped to metabolic and regulatory pathways in KEGG. The proteomic data were deposited in the iProX platform (http://www.iprox.org) with the project no. IPX0003502000.

Partial of total proteins same as proteome sequencing were used for further targeted quantification by Parallel Reaction Monitoring (PRM) by Beijing QLBio Biotechnology Co. (China). Briefly, An AQUA stable-isotope peptide as an internal standard reference was spiked in with each sample. Digested peptides were desalted on C₁₈ stage tips prior to reversed phase chromatography on an Easy-nanoLC system (Thermal Scientific, Chelmsford, MA, USA). One hour liquid chromatography at a flow rate of 300 nL/min was used with the following gradients: 3 to 28% buffer B in 107min and 28 to 80% buffer B in 3 min. PRM analyses was performed using an orbitrap fusion mass spectrometer (Thermal Scientific, Chelmsford, MA, USA). Optimal collision energy, charge state, and retention time for the most significantly regulated peptides were generated experimentally using unique peptides of high intensity and confidence for each target protein. The mass spectrometer was operated in position ion mode with the following parameters: the full scan was collected with a resolution of 120,000 at 200m/z, the AGC target was 2×10⁶ and the maximum injection time was at 100 ms. All PRM data analyses and data integration were performed using Skyline version 3.5.0 (MacLean et al., 2010). Three replicates were included for each sample in the PRM-MS analyses. Relative peptide quantification was calculated by dividing the peptide peak area. A two-tailed Student’s t-test was used to estimate the significance of the difference in relative peptide abundance between M. orientalis adults and metacercariae.

In order to examine the detail post-translational regulation between transcriptome and proteome in adult and metacercaria M. orientalis, the fold changes of mRNA and protein were compared. In briefly, the fold changes were got in mRNA and protein level separately between adult stage and metacercariae-stage M. orientalis. And then, fold change ratio was calculated (As the following formula).

We propose for most genes, the fold changes of mRNA and protein are similar. In order to get those genes with significant different fold change, significance A were calculate using MaxQuant version 1.5.3.30 (Wiśniewski et al., 2009), and those genes with the p value less than 0.05 and fold change ratio greater than 2 or less than 2 were consider as the significance up or down genes.

The statistical analyses involved in this study was implemented on R 4.0.3 platform (https://www.r-project.org/). Visualization of graphics were built in GraphPad Prism 8.0 (GraphPad Software, La Jolla, CA, USA) and R environment using the ggplot2 version 3.3.5 and ggpubr version 0.4.0 by online website https://github.com/kassambara/ggpubr/.

A total of 47,396,124 and 51,014,100 clean reads were obtained from M. orientalis adults and metacercariae, respectively. The average ratio of clean reads to raw reads was 94.43% ( Table S2 ). A total of 254,543 (>200 nt) transcripts were produced by the Trinity program for all samples. The N50 size was 1,002 nucleotide base pairs, which is shorter than other helminth transcriptomes (Liu et al., 2016; Zhang et al., 2017). Removing the redundancy resulted in 100,234 unigenes based on at least 0.3 FPKM in all samples. Functional annotation of the 100,234 unigenes data set was carried out using seven in public databases (Nr, Nt, Pfam, COG, Swiss-Prot, KEGG, and GO), resulting in 58.27% (n = 58,402) of the data set being annotated in at least one database ( Table 1 ). Using a cut-off FDR of < 0.005 and a twofold change identified 13,823 DEGs in adult versus metacercariae, of which 4,773 were upregulated and 9,050 were downregulated ( Table S3 ; Figure 1A ).

In total, 21,604 peptides, 10,704 unique peptides, and 3,530 proteins were determined via proteomic analyses ( Table S4 ). In terms of protein mass distribution, most proteins (73.5%) had molecular weights ranging from 10 to 70 ku, 12.1% of proteins had a molecular weight > 100 ku. Using KEGG and GO database annotations, 45.98% of proteins (1,623/3,530) were annotated to 1,036 GO terms and 67.56% proteins (2,385/3,530) were annotated to 1,140 KEGG pathways ( Table 1 ). Using a twofold change and p-value < 0.05 as a threshold resulted in 1,445 DEPs detected by the proteomic analyses (adult versus metacercariae), among which 519 proteins were up-regulated and 926 proteins were down-regulated ( Table S5 ; Figure 1B ).

By comparing the RNA-seq data with the proteomic data, all proteins were matched to the corresponding unigenes. Of these, 2,225 genes, including 780 DEGs and 1445 DEPs displayed differential expression at either the mRNA or protein levels between adult stage and metacercariae-stage M. orientalis. Of these, 496 genes were consistently upregulated between the adult stage and metacercariae stage, and 74 genes were consistently downregulated between the adult stage and metacercariae stage ( Table 2 ; Figure 1C ). However, 1,655 genes showed inconsistent expression at the mRNA and protein levels between the adult stage and metacercariae stage, which might result from post-translational regulation or modifications.

Analyses of gene sharing between the adult stage and metacercariae-stage M. orientalis revealed 34,228 unigenes to be co-expressed by adults and metacercariae, whereas most (n = 66,006) exhibited stage-specific expression (adult = 25,426; metacercariae = 40,850) ( Figure 2A ). This suggested that the transcriptome level during the metacercariae stage is high, which might be related to the infection ability of metacercariae. There were 1,272 terms common to both adults and metacercariae. Only 81 GO terms were unique to the adults, compared with 998 GO terms in the metacercariae. The unique GO terms in the adults included GO:0044163, GO:0075521, and GO:0071479, involving host cytoskeleton, microtubule-dependent intracellular transport of viral material towards nucleus, and cellular response to ionizing radiation. Specific GO entries in the metacercariae included GO:0004664, GO:0004298, GO:0005839, and GO:0051603. These functions involved prephenate dehydratase activity, threonine-type endopeptidase activity, proteasome core complex, and proteolysis involved in cellular protein catabolism.

GO analyses were performed on DEGs and all DEGs were assigned to 4,389 GO terms, of which 446 GO terms were significantly enriched (p < 0.05). We focused on these enriched functions. In terms of biological processes, 267 GO terms were associated with these DEGs, most of which were involved in cellular protein metabolic process, cellular protein modification process, protein modification process. In terms of the cellular component, 77 GO terms were associated with these DEGs and were mainly involved in organelle, organelle part, and intracellular organelle part. By contrast, 102 molecular function GO terms were associated with these DEGs, including protein binding, binding, nucleoside binding, GTPase activity, and purine ribonucleoside triphosphate binding ( Figure 2B ; Table S6 ). In addition, nucleocytoplasmic transport (GO:0006913) was the only function that was significantly enriched among the genes unique to adults. By contrast, 44 GO terms were significantly enriched that involved genes unique to metacercariae ( Figure 2C ).

The up- and down-regulated genes were further analyzed to better understand the DEGs between adult and metacercariae-stage M. orientalis. All the upregulated genes were assigned to 139 GO terms ( Table S7 ). Among these, up-regulated genes were uniquely enriched in 64 GO functions, mainly related to microtubules, and were significantly enriched in functions such as microtubule-based processes, microtubule-based movement, microtubule motor activity, and tubulin binding ( Figure 2D ). In addition, numerous up-regulated genes were enriched in GO functions of catalytic activity, hydrolase activity, and catalytic complexes. By contrast, down-regulated genes were assigned to 234 GO terms, of which 159 were uniquely enriched, including transmembrane receptor protein serine/threonine kinase activity, transmembrane receptor protein serine/threonine kinase signaling pathway, DNA integration, transmembrane receptor protein kinase activity, and enzyme linked receptor protein signaling pathway ( Figure 2E ). In addition, numerous down-regulated genes were concentrated in primary metabolic processes, cellular metabolic processes, macromolecule metabolic processes, and cellular macromolecule metabolic processes. Interestingly, a few genes participated in a variety of negatively regulated functions, such as neuropeptide hormone activity, negative regulation of proteasomal protein catabolic processes, and negative regulation of mitotic nuclear division.

KEGG pathway enrichment analyses of DEGs showed that 215 up-regulated DEGs were significantly enriched in 35 pathways (p <0.05) ( Figure 3A ; Table S8 ). Consistent with the GO analyses, up-regulated genes were enriched in a variety of pathways related to reproduction, substance metabolism, and biosynthesis, such as oocyte meiosis (ko04114), glycolysis/gluconeogenesis (ko00010), pyrimidine metabolism (ko00240), nitrogen metabolism (ko00910), N-glycan biosynthesis (ko00510), purine metabolism (ko00230), starch and sucrose metabolism (ko00500) and arginine biosynthesis (ko00220). Similarly, various up-regulated genes participated in pathways related to genetic information processing, suggesting vigorous reproductive behavior during the adult stage. Of note, up-regulated genes were significantly enriched in the ko05130 pathway, which is mainly involved in the encoding of ACTB_G1 (Actin beta/Gamma 1), a protein with a key role in adult motion. In addition, these genes were enriched in the p53 signaling pathway (ko04115), which is highly associated with cancers and mainly regulates the apoptosis and senescence of cells. Interestingly, genes that were highly expressed in adults were also involved in butirosin and neomycin biosynthesis (ko00524), thyroid hormone synthesis (ko04918), inflammatory mediator regulation of TRP channels (ko04750), and drug metabolism - cytochrome P450 (ko00982). These findings could provide clues to the mechanisms by which the parasites fight bacteria and immune evasion.

Of the down-regulated DEGs, 261 were significantly enriched in 26 pathways (p < 0.05) ( Figure 3B ). Most of these genes expressed in metacercariae were involved in various biological cycles and some metabolic processes, such as the citrate cycle (TCA cycle) (ko00020), oxidative phosphorylation (ko00190), pyruvate metabolism (ko00620), and one-carbon pool by folate (ko00670). In addition, many hyperexpressed genes in metacercariae were enriched in the AMPK signaling pathway (ko04152), a fuel sensor and regulator that promotes ATP-producing and inhibits ATP-consuming pathways in various tissues. This suggests that the parasite regulates its energy use during the metacercariae stage when food intake is not possible. Interestingly, the overexpressed metacercariae genes were also involved in a variety of disease-related pathways, such as herpes simplex infection (ko05168), epithelial cell signaling in Helicobacter pylori infection (ko05120), bacterial invasion of epithelial cells (ko05100), non-alcoholic fatty liver disease (NAFLD) (ko04932), viral carcinogenesis (ko05203), and salmonella infection (ko05132). These findings could provide clues to a potential association between the M. orientalis metacercariae and these diseases.

To validate the transcriptome data, eight genes (three up-regulated and five down-regulated in adult stage versus metacercariae stage) were selected randomly among the DEGs, and their expression levels were verified by qRT-PCR ( Figure S1 ). The results showed that DEGs of phosphoglycerate mutase, translation initiation factor 3 subunit B, and cathepsin F precursor were up-regulated in the adult stage versus metacercariae stage, whereas cytoplasmic 1, serpin, glutamine synthetase, elongation factor 1-gamma, and lactate dehydrogenase were up-regulated in metacercariae stage versus adult M. orientalis, which showed a similar expression trend to the transcriptome analysis, providing evidence on the reliability of the transcriptome sequencing results.

GO analyses were performed on DEPs, and all DEPs were assigned to 630 GO functions. Among these, there were 218, 101, and 311 GO terms for biological processes, cellular components and molecular functions, respectively ( Figure 4A ). Enrichment analyses showed that only protein folding (GO:0006457) was significantly enriched (FDR < 0.05) ( Figure 4B ; Table S9 ). The DEPs up- and down-regulated between the adult and metacercariae M. orientalis were then analyzed. In terms of biological processes, up-regulated proteins were significantly enriched in protein folding, microtubule-based processes and translation, whereas down-regulated proteins were significantly enriched in carbohydrate metabolic processes, TCA cycle, and the pentose-phosphate shunt ( Figure 4C ). In terms of cellular components, up-regulated proteins were significantly enriched in microtubules and ribosomes, whereas down-regulated proteins were significantly enriched in extracellular space, integral membrane components, and synaptic vesicles. In terms of molecular functions, up-regulated proteins were significantly concentrated in the unfold protein binding, structural dynamics of ribosomes, protein kinase activity, and structural constituents of cytoskeleton categories, whereas down-regulated proteins were assigned to 205 GO terms ( Figure 4C ), however, no protein were enriched for any function.

KEGG pathway enrichment analyses were performed for DEPs, and all DEPs were assigned to 335 KEGG pathways ( Table S10 ). These proteins were significantly enriched in glutathione metabolism (ko00480) and glycerolipid metabolism (ko00561). Up-regulated and down-regulated proteins enrichment analyses showed that a total of 388 up-regulated DEPs were involved in 294 KEGG pathways, and 265 down-regulated DEPs were involved in 304 pathways. The up-regulated proteins were significantly enriched in protein processing in endoplasmic reticulum, ribosomes, gap junctions, antigen processing and presentation, and drug metabolism ( Figure 4D ). By contrast, the down-regulated proteins were significantly enriched in 29 pathways, among which numerous proteins were enriched in the biosynthesis of antibiotics, biosynthesis of amino acids, valine, leucine and isoleucine degradation pathways ( Figure 4D ). In addition, a few down-regulated proteins were enriched in inflammatory mediator regulation of TRP channels, PPAR signaling pathways, and cAMP signaling pathways.

PRM is an ion monitoring technology based on high resolution and high precision mass spectrometry, and can selectively detect target proteins and peptides, such as post-translational modifications peptides. To confirm DEPs in label-free analyses between adult stage and metacercariae-stage M. orientalis, eight significant DEPs were selected randomly for PRM analyses. The results showed the similar expression trends for the label-free proteomics data and PRM data ( Table 3 ), verifying the accuracy and reliability of the proteome analyses.

In order to examine the detail post-translational regulation between transcriptome and proteome in adult and metacercaria M. orientalis, the significance up or down genes were calculated by significance A method ( Table S11 ). Further analyses of the significance up or down genes were conducted through GO and KEGG pathways. Analyses of potential translation regulation were performed to obtain possible up-regulation and down-regulation results. Among 669 up-regulated genes, 230, 213, and 364 genes were assigned to biological processes, cellular components and molecular functions, respectively ( Figure 5A ; Table S12 ). These up-regulated genes were associated with actin filament capping, spectrin, glucose metabolic processes, arginine metabolic processes, structural constituents of cytoskeleton, myosin complexes and ATP binding ( Figure 5B ). In addition, some up-regulated genes were enriched in urea cycle, glucose metabolic processes, aromatic amino acid family metabolic processes, TCA cycle, and glycolytic processes ( Table S13 ). This suggested that energy metabolism activities increased during metacercariae development. In addition, these up-regulated genes were also involved in the phagosome, spliceosome, and PPAR signaling and Hippo signaling pathways. Interestingly, many of the up-regulated genes were also involved in a variety of disease processes, such as hypertrophic cardiomyopathy (HCM), dilated cardiomyopathy (DCM), Vibrio cholerae infection, and viral myocarditis. Among 235 down-regulated genes, 53, 55, and 83 genes were assigned to biological processes, cellular components, and molecular functions, respectively. These decreased genes were associated with cilium assembly and movement, motile cilia, microtubule pyridoxal phosphate binding, and Rab GTPase binding ( Figure 5C ). Some down-regulated genes were also enriched in gap junctions, glucagon signaling pathway, aldosterone synthesis and secretion, insulin secretion, and other pathways ( Figure 5D ). These results suggest that although gene and protein expression levels were not always consistent, there was a high degree of consistency between the functions of DEGs and DEPs in the different stages of M. orientalis.