Identification and Expression Profiling of Circulating MicroRNAs in Serum of Cysticercus pisiformis-Infected Rabbits

Cysticercus pisiformis (C. pisiformis), the larval form of Taenia pisiformis, parasitize mainly the liver, omentum and mesentery of rabbits and cause huge economic losses in the rabbit breeding industry. MicroRNA (miRNA), a short non-coding RNA, is widely and stably distributed in the plasma and serum. Numerous data demonstrates that, after parasitic infection, miRNAs become the key regulatory factor for controlling host biological processes. However, the roles of serum miRNAs in C. pisiformis-infected rabbits have not been elucidated. In this study, we compared miRNA expression profiles between the C. pisiformis-infected and healthy rabbit serum using RNA-seq. A total of 192 miRNAs were differentially expressed (fold change ≥ 2 and p < 0.05), including 79 up- and 113 downregulated miRNAs. These data were verified by qRT-PCR (real time quantitative polymerase chain reaction) analysis. Additionally, GO analysis showed that the target genes of these dysregulated miRNAs were most enriched in cellular, single-organism and metabolic processes. KEGG pathway analysis showed that these miRNAs target genes were involved in PI3K-Akt, viral carcinogenesis and B cell receptor signaling pathways. Interestingly, after aligning clean reads to the T. pisiformis genome, four (miR-124-3p_3, miR-124-3p_4, miR-124a and novel-miR1) T. pisiformis-derived miRNAs were found. Of these, novel-miR1was upregulated in different periods after C. pisiformis infection, which was verified qRT-PCR, and pre- novel-miR-1 was amplified from the cysticerci by RT-PCR, implying novel-miR-1 was derived from C. pisiformis and has great potential for the diagnosis of Cysticercosis pisiformis infection. This is the first investigation of miRNA expression profile and function in the serum of rabbits infected by C. pisiformis, providing fundamental data for developing diagnostic targets for Cysticercosis pisiformis.


Introduction
Cysticercus pisiformis (C. pisiformis), the larval form of Taenia pisiformis (T. pisiformis), is a heteroxenous parasite and belongs to Taenidae [1]. T. pisiformis parasitizes mainly the small intestines of canines, while C. pisiformis parasitizes mainly the omentum, liver and mesentery of lagomorphs and rodents [1][2][3]. T. pisiformis and C. pisiformis are distributed all over the world. In Cairo, Egypt, the prevalence of T. pisiformis infection in dogs was 70.8% [4]; in northern Italy, the cysticercosis of brown hares has been reported to be 3.28-14.8% [2]. However, in the state of Morelos, Mexico, the infection rate of hare cysticercosis is as high as 67.7% [5]. In China, T. pisiformis has also been reported [6]. Once infected with C. pisiformis, rabbits will show decreased prolificacy, weight loss and weakened immunological resistance, causing serious economic losses to the rabbit breeding industry [7,8]. The Dot-ELISA diagnostic method has been reported for detection of C. pisiformis, but has high false positive rate limit applicability [9]. Multiplex PCR has been studied as a diagnostic method for T. pisiformis in dog feces, but cannot be used for the diagnosis of rabbit cysticercosis [10]. Thus, there is an urgent need to develop a new method for early diagnosis for rabbit in order to reduce the losses due to cysticercosis.
MiRNAs are a class of endogenous non-coding RNA ranging in length from 18-24 nt [11]. Owing to the short length, miRNAs exist stably in various biosamples. miRNAs negatively regulating the expression of coding RNA at the post-transcriptional level have demonstrated that miRNAs play important roles in many diseases [12][13][14]. Some differentially expressed miRNAs have been proven as diagnostic targets in various diseases [15,16]. Numerous studies have shown that complex immune responses against parasites is activated in the host [17][18][19]. Of them, miRNAs have important functions in host-parasite interactions. For example, in Cryptosporidium parvum infection mouse model, the hostderived miR-27 suppressed the expression of KH-type splicing regulatory protein (KSRP) to promote iNOS mRNA stability and activate the TLR4/NF-κB signaling indirectly, which is involved in the epithelial antimicrobial defense process [17]. miR-101b-3p, expressed in mouse brain tissue infected with Angiostrongylus cantonensis (A. cantonensis), can suppress extracellular superoxide dismutase 3 (Acsod3) expression, thus increasing miR-101b-3p expression in mice that can ameliorate the oxidative damages of larvae [20]. miR-155-5p was shown to be involved in the central nervous system (CNS) inflammation induced by A. cantonensis infection [21]. Additionally, 58 host-origin serum miRNAs were significantly differentially expressed in mice infected by E. multilocularis [22]. Taken together, the above studies reveal the key role of miRNAs in the host response against parasitic infection and imply that dysregulated miRNAs could be involved in the biological processes of parasitic diseases.
Parasite-derived circulating miRNAs have shown great potential as diagnostic targets for some diseases. For instance, sja-miR-2b-5p and sja-miR-2c-5p were found to be significantly differentially expressed in patients infected by Schistosoma japonicum (S. japonicum) and can be used as biomarkers for the detection of human S. japonicum infection [23]. The differential expression of egr-miR-71 and egr-let-7 in the serum of cystic echinococcosis patients can also be promising biomarkers for early detection of hydatid cyst infection and post-surgical follow-up [24].
However, to the best of our knowledge, there is no data regarding the serum miRNA expression profile in rabbits infected with C. pisiformis. In the present study, we analyzed the miRNA expression profiles in C. pisiformis-infected and uninfected rabbit serum. Some dysregulated host and C. pisiformis-derived miRNAs were found and analyzed, which will improve understanding of the interactions between rabbits and C. pisiformis and provide fundamental miRNA data for screening serological diagnostic targets of cysticercosis pisiformis.

Ethics Statement
This study was approved by the Animal Ethics Committee of Lanzhou Veterinary Research Institute, Chinese Academy of Agricultural Sciences (Permit No. LVRIAEC-2016-006), and the protocols used were strictly in accordance with the guidelines of animal welfare. The minimum amount of harm to experimental rabbits was required in the present study.

Parasites and Animals
C. pisiformis was isolated from naturally infected rabbits at a slaughterhouse in Jiaozuo city, Henan province, China. To obtain abundant C. pisiformis, each of the three 2-month-old pathogen-free dogs was orally challenged with 20 mature C. pisiformis. At 60 days post infection (PI), the eggs of T. pisiformis were collected from the feces of the C. pisiformisinfected dogs, diluted with PBS to 1000 eggs/ml. Subsequently, six 1.5-month-old pathogenfree New Zealand white rabbits were randomly assigned into infection groups (Cpi groups, which including Cpi-1, Cpi-2 and Cpi-3) and control groups (NC groups, which included NC-1, NC-2 and NC-3). Each rabbit in the Cpi group was orally infected with 1 mL (1000 eggs) mature T. pisiformis eggs while the NC group was orally treated with 1 mL PBS. The rabbits in both groups were reared separately with forage and adequate lukewarm boiled water. The successful infection was determined by autopsy and examination of the number of C. pisiformis.

Blood Samples Collection and RNA Isolation
Fresh blood was collected from the rabbits in the Cpi and NC groups at 60 days PI and the total RNA was isolated from the separated serum samples using Trizol Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. Furthermore, the six rabbits were necropsied to confirm infection by C. pisiformis. The results of the examination are shown in Table 1.

Small RNA Library Generation and Sequencing
Small RNA fragments were separated from the six total RNA samples by polyacrylamide gel electrophoresis (PAGE), and RNA with the lengths of 18-30 nt was collected. Subsequently, 5-adenylated and 3-blocked adaptors were ligated to the 3 end of the small RNA fragments, then unique molecular identifiers (UMI) [25] labeled primers were added to hybridize with the 3 adaptor. Unlinked adaptors were enzyme digested after 5 end adaptor ligation. cDNA first strand was synthesized using RT primers for the UMI. The cDNA library containing both 5 and 3 adaptors were amplified by PCR and products 100-120 bp in length were selected to construct small RNA libraries. After library quantitation, pooling cyclization and quality control, the qualified libraries were prepared for RNA sequencing. The qualified sequences were sequenced on a BGISEQ-500 (Beijing Genomics Institute, China) and the sequences (raw reads) were produced for subsequent analysis.

Reads Analysis and MiRNA Prediction
The clean reads were obtained after removing the low-quality reads (those with 5 primer contaminants, without 3 primer and/or without insertion) from the raw reads. The clean reads were mapped against the Oryctolagus cuniculus genome (https://www. uniprot.org/uniprot/?query=taxonomy:9986 (accessed on 19 December 2019)) and T. pisiformis genome (unreleased) using Bowtie2 software. The mapped reads were classified and annotated into different kinds of small RNAs (miRNAs, snRNA, snoRNA, Rfam other sncRNA, rRNA, tRNA). Of them, the miRNA maturity and precursor were selected by comparing against the known mature miRNA in the miRbase database Unknown tags were used to predict novel-miRNAs through miRDeep2 by exploring the characteristic hairpin structure of the miRNA precursor and the minimum free energy (MFE) of the candidate miRNA secondary structure [26,27].

Verification of Differential Expression MiRNAs by Using qRT-PCR
The miRNAs from the fresh blood samples were separated using an EasyPure miRNA Kit (TransGen, Beijing, China) according to the manufacturer's instructions. Subsequently, the miRNAs were reverse transcribed into cDNA using the Mir-XTM miRNA First-Strand Synthesis Kit (Clontech, Shanghai, China). qRT-PCR was performed in a 25 µL reac-tion volume, including 12.5 µL SYBR Advantage Premix (2×) (Takara, Shiga-ken, Japan), 2 µL cDNA, 0.5 µL ROX Dye (50×), 0.5 µL miRNA-specific forward primer, 0.5 µL mRQ 3 primer and 9 µL RNase free water. qRT-PCR was performed with the following conditions: 95 • C for 10 min and 40 cycles of 95 • C for 15 s, 60 • C for 1 min. The expression of 8 target miRNAs in rabbit serum at two months PI with T. pisiformis were normalized to U6 small nuclear RNA (snRNA), with the primers (Sangon Biotech, Shanghai, China) listed in Table 2. The 2 −∆∆Ct method was used to calculate the relative expression of candidate miRNAs [28].

Function Prediction of Dysregulated MiRNAs
To explore the potential functions of differentially expressed miRNAs (DEMs), target genes were predicted using RNAhybrid [29] and miRanda [30] software. The target genes were functionally annotated by assigning Gene Ontology terms (GO-terms, http: //www.geneontology.org/ (accessed on 29 November 2019)) [31] and KEGG pathway analysis [32]. By calculating the number of genes per-term and applying the hypergeometric test, the DEM target gene functions were subjected to GO function classification. Similarly, KEGG was used to perform pathway enrichment analysis.

Statistical Analysis
Student's t-test was used to analyze the significance of miRNA expression and Graphpad 7 was used to create the relative expression diagram. p < 0.05 was considered statistically significant.

Characterization of the sRNAs Libraries
Among the six small RNA (sRNAs) libraries, the average raw tags were 28,000,000 per library. After removing low-quality tags in each library, about 26,000,000 clean tags were obtained, which accounted for greater than 90% of the raw tags. Clean tags were aligned to the reference genome of Oryctolagus cuniculus (https://www.uniprot.org/uniprot/?query= taxonomy:9986 (accessed on 19 December 2019)) and approximately 23,000,000 tags were mapped for mapping rates between 85.11% to 94.11%. In order to determine whether or not there were sRNAs from C. pisiformis, we mapped the clean tags to the T. pisiformis genome (unreleased). The result showed that 1,337,063 tags were aligned to the T. pisiformis genome for mapping rates from 2.6% to 5.27%. The specific data of each group are shown in Table 1.
The tags which mapped to the Oryctolagus cuniculus genome were annotated into different kinds of sRNAs (miRNAs, snRNA, snoRNA, Rfam other sncRNA, rRNA, tRNA), genes (exon, intron, intergenic), unmapped and repeat sequences. Of them, 3.15% were classified as miRNAs (mature, hairpin and precursor) (Figure 1a). When mapped to the T. pisiformis genome, the proportion of miRNAs was 6.17% ( Figure 1b). Additionally, sRNA populations of between 17 and 26 nt in length (with the most abundant length being 22 nt and 23 nt) were observed in both species (Figure 2), which is consistent with miRNA common length.
The tags which mapped to the Oryctolagus cuniculus genome were annotated into different kinds of sRNAs (miRNAs, snRNA, snoRNA, Rfam other sncRNA, rRNA, tRNA), genes (exon, intron, intergenic), unmapped and repeat sequences. Of them, 3.15% were classified as miRNAs (mature, hairpin and precursor) (Figure 1a). When mapped to the T. pisiformis genome, the proportion of miRNAs was 6.17% ( Figure 1b). Additionally, sRNA populations of between 17 and 26 nt in length (with the most abundant length being 22 nt and 23 nt) were observed in both species (Figure 2), which is consistent with miRNA common length.  It is worth mentioning that in the present study, the results of Cpi-3 and NC-3 samples were quite different from the other samples. To obtain accurate results, we analyzed only two sets of data of Cpi-1/ Cpi-2 and NC-1/ NC-2 in the next study.

Identification of Differentially Expressed miRNAs
By comparison, a total of 530 rabbit-derived miRNAs were found to be differentially expressed in the Cpi-1 and NC-1 groups (Supplementary Table S1), and 548 in the Cpi-2 and NC-2 groups (Supplementary Table S2). Among them, 192 duplicated miRNAs (71 upregulated and 121 downregulated;|log2(Cpi/NC)| > 1 and FDR < 0.001; Supplementary  Table S3) were significantly differentially expressed in the two infected groups. Of these 192 miRNAs, 20 known miRNAs and 172 novel predicted miRNAs were included. The sequences and reads of the 20 known miRNAs are shown in Figure 3. It is worth mentioning that in the present study, the results of Cpi-3 and NC-3 samples were quite different from the other samples. To obtain accurate results, we analyzed only two sets of data of Cpi-1/ Cpi-2 and NC-1/ NC-2 in the next study.

Identification of Differentially Expressed miRNAs
By comparison, a total of 530 rabbit-derived miRNAs were found to be differentially expressed in the Cpi-1 and NC-1 groups (Supplementary Table S1), and 548 in the Cpi-2 and NC-2 groups (Supplementary Table S2). Among them, 192 duplicated miRNAs (71 upregulated and 121 downregulated;|log2(Cpi/NC)| > 1 and FDR < 0.001; Supplementary Table S3) were significantly differentially expressed in the two infected groups. Of these 192 miRNAs, 20 known miRNAs and 172 novel predicted miRNAs were included. The sequences and reads of the 20 known miRNAs are shown in Figure 3.  Interestingly, four deregulated T. pisiformis-derived miRNAs (miR-124-3p_3, miR-124-3p_4, miR-124a and novel-miR1) were found in the two Cpi groups (Supplementary  Table S4 and S5). However, novel-miR1 was the only upregulated miRNA, which was verified by qRT-PCR at different time points after infection. The stem-loop structure for the putative pre-miRNAs of novel-miR1 was predicted by miRdeep2 and the sequence of pre-novel-miR1 was specifically amplified from the C. pisiformis (shown in Figure 4). Further analysis found that miR-124-3p_3, miR-124-3p_4, miR-124a were also present in the rabbit samples. The miRNA distribution of each group is shown in Figure 5. Interestingly, four deregulated T. pisiformis-derived miRNAs (miR-124-3p_3, miR-124-3p_4, miR-124a and novel-miR1) were found in the two Cpi groups (Supplementary  Table S4 and S5). However, novel-miR1 was the only upregulated miRNA, which was verified by qRT-PCR at different time points after infection. The stem-loop structure for the putative pre-miRNAs of novel-miR1 was predicted by miRdeep2 and the sequence of pre-novel-miR1 was specifically amplified from the C. pisiformis (shown in Figure 4). Further analysis found that miR-124-3p_3, miR-124-3p_4, miR-124a were also present in the rabbit samples. The miRNA distribution of each group is shown in Figure 5.

qRT-PCR Verification of miRNA Expression
Eight differentially expressed miRNAs from the rabbits were randomly selected for qRT-PCR to confirm sequencing data. The results showed that the expression trend of all candidate miRNAs, seven upregulated (novel-mir128, novel-mir303, novel-mir307, novel-mir318, novel-mir476, novel-mir857 and novel-mir150) and one downregulated miRNA (novel-mir57) were consistent with the sequencing results, which support the accuracy of miRNA sequencing ( Figure 6).

qRT-PCR Verification of miRNA Expression
Eight differentially expressed miRNAs from the rabbits were randomly selected for qRT-PCR to confirm sequencing data. The results showed that the expression trend of all candidate miRNAs, seven upregulated (novel-mir128, novel-mir303, novel-mir307, novel-

Prediction of miRNA target genes
In order to further identify the importance of miRNA/mRNA-regulated axes in T. pisiformis infection, RNAhybrid [29] and miRanda [30] were used to predict the DEM target genes. A total 2,202,871 genes were predicted as target genes of 886 differentially expressed rabbit miRNAs (Supplementary Table S6). A total of 140 target genes were predicted from four T. pisiformis-derived miRNAs (Supplementary Table S7). In order to confirm these findings more definitively, 10 rabbit DEMs (5 up-and 5 downregulated miR-NAs) were selected to construct an interaction network using Cytoscape_3.7.2, which included 348 nodes and 412 edges (Figure 7a). Additionally, four T. pisiformis-derived miR-NAs were used to construct an interaction network that included 140 nodes (Figure 7b). From the network diagram, it was clear that the target XM.008265418.2 (tyrosine kinase with immunoglobulin-like and EGF-like domains 1, TIE1) was regulated by the novel-miR388, novel-miR62 and novel-miR635. Additionally, the target genes of miR-124a such as Mark00009255, Mark00000443, Mark00007007 and Mark00009550 can be clearly seen through the network diagram.

Prediction of miRNA Target Genes
In order to further identify the importance of miRNA/mRNA-regulated axes in T. pisiformis infection, RNAhybrid [29] and miRanda [30] were used to predict the DEM target genes. A total 2,202,871 genes were predicted as target genes of 886 differentially expressed rabbit miRNAs (Supplementary Table S6). A total of 140 target genes were predicted from four T. pisiformis-derived miRNAs (Supplementary Table S7). In order to confirm these findings more definitively, 10 rabbit DEMs (5 up-and 5 downregulated miRNAs) were selected to construct an interaction network using Cytoscape_3.7.2, which included 348 nodes and 412 edges (Figure 7a). Additionally, four T. pisiformis-derived miRNAs were used to construct an interaction network that included 140 nodes (Figure 7b). From the network diagram, it was clear that the target XM.008265418.2 (tyrosine kinase with immunoglobulin-like and EGF-like domains 1, TIE1) was regulated by the novel-miR388, novel-miR62 and novel-miR635. Additionally, the target genes of miR-124a such as Mark00009255, Mark00000443, Mark00007007 and Mark00009550 can be clearly seen through the network diagram.

Functional Analysis
To determine the roles of the rabbit DEMs, the top 10 terms (biological processes, cellular component and molecular function) in the Cpi and NC groups were obtained by gene ontology term analysis (Figure 8a). Among them, the most abundant GO terms under the biological process (BP) category were cellular process, single-organism process and metabolic process. The most enriched cellular component (CC) terms were organelle, cell part and membrane. Finally, molecular transducer activity, catalytic activity and binding were the most enriched molecular function (MF) terms. The KEGG pathway analysis of potential target genes were involved in some of the signaling pathways responsible for tight junction, PI3K-Akt signaling pathway, microRNAs in cancer, cell adhesion molecules and apoptosis (Figure 8b).
The target genes of T. pisiformis miRNAs enriched by GO terms of BP, MF and CC are shown in Figure 9a. The notable GO terms were cellular process, metabolic process and catalytic activity. KEGG pathway analysis showed that the majority of the target genes of the T. pisiformis-derived miRNA participated in viral carcinogenesis, osteoclast differentiation, herpes simplex infection and B cell receptor signaling pathway (Figure 9b). The target genes of T. pisiformis miRNAs enriched by GO terms of BP, MF and CC are shown in Figure 9a. The notable GO terms were cellular process, metabolic process and catalytic activity. KEGG pathway analysis showed that the majority of the target genes of the T. pisiformis-derived miRNA participated in viral carcinogenesis, osteoclast differentiation, herpes simplex infection and B cell receptor signaling pathway ( Figure  9b).
In addition, the miRNAs were relatively stable in body fluids and tissues [39]. Many miRNAs have been proven as biomarkers with great potential in the diagnosis and treatment of infectious diseases and cancer [40][41][42]. Of course, there are many studies about miRNAs as biomarkers to diagnose or detect parasitic infection. For instance, miR-277, miR-3479-3p and bantam, in human serum from Schistosome endemic areas, were reported to detect infected S. mansoni individuals from low and high infection intensity sites with specificity/sensitivity values of 89%/80% and 80%/90%, respectively [43]. Additionally, parasite-derived miRNAs sja-miR-2b-5p, sja-miR-2c-5p, sja-miR-277 and sja-miR-3479-3p in mouse serum were shown to have potential as diagnostic markers for Schistosomiasis japonica [23,44]. Previous research has shown that hydatidosis can be diagnosed early and monitored through the biomarker of serum-derived egr-miR-71 [24]. The expression of hsa-miR-125b-5p was stably upregulated in the plasma and liver tissue samples from patients with alveolar echinococcosis (AE), which suggested that hsa-miR-125b-5p may be a promising biomarker for early non-invasive diagnosis of AE [45]. Serum aca-miR-146a was found to be expressed significantly higher in mice infected by Angiostrongylus cantonensis and its receiver operating characteristic (ROC) curve analysis showed that aca-miR-146a was an effective biomarker for discriminating the infected from uninfected cases with an area under the ROC curve (AUC) of 0.90. Its diagnostic accuracy was assessed on patients (n = 30) and healthy controls (n = 30), and the sensitivity and specificity reached 83% and 86.7%, respectively, indicating that the aca-miR-146a in serum is an effective biomarker to track infection of A. cantonensis [46].
Furthermore, in the present study, we amplified the precursor sequence of novel-miR1 from C. pisiformis and found that novel-miR1 is unique to C. pisiformis. qRT-PCR analysis revealed that novel-miR-1 was upregulated in the serum of rabbits infected with C. pisiformis 2, 3 and 7 months after infection ( Figure 4). As the candidate diagnosis target of Cysticercosis pisiform, novel-miR1 has a strong specificity and should be further investigated.

Conclusions
To the best of our knowledge, this study was the first to determine miRNA expression pattern in rabbits infected with C. pisiformis. A total of 192 DEMs were identified, including 20 known and 172 novel-miRNAs. Additionally, four T. pisiformis miRNAs were found to be differentially expressed, but only novel-miR1 was derived strictly from C. pisiformis and upregulated in different periods in rabbits infected with T. pisiformis. Further miRNAs data analysis showed that these deregulated miRNAs were most enriched in cellular processes, single-organism processes, metabolic processes and molecular transducer activity. In addition, the signaling pathways of those miRNA were involved in PI3K-Akt, viral carcinogenesis and B cell receptor signaling pathways. Most importantly, novel-miR1 shows promising diagnostic potential for cysticercosis pisiformis.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/genes12101591/s1, Table S1: Rabbit miRNAs expressed differentially between the Cpi-1 and the NC-1 groups; Table S2: Rabbit miRNAs expressed differentially between the Cpi-2 and the NC-2 groups; Table S3: Duplicate rabbit miRNAs differentially expressed between the Cpi1 vs. NC-1 groups and Cpi2 vs. NC-2 groups; Table S4: T. pisiformis miRNAs differentially expressed between the Cpi-1 and the NC-1 groups; Table S5: T. pisiformis miRNAs differentially expressed between the Cpi-2 and the NC-2 groups; Table S6: Total target genes were predicted for differentially expressed rabbit miRNAs; Table S7: Total target genes were predicted for differentially expressed T. pisiformis miRNAs.