Fasciola gigantica–Derived Excretory-Secretory Products Alter the Expression of mRNAs, miRNAs, lncRNAs, and circRNAs Involved in the Immune Response and Metabolism in Goat Peripheral Blood Mononuclear Cells

Fasciola gigantica produces excretory-secretory products (ESPs) with immune-modulating effects to promote its own survival. In this study, we performed RNA-seq to gain a comprehensive global understanding of changes in the expression of mRNAs, miRNAs, lncRNAs, and circRNAs in goat peripheral blood mononuclear cells (PBMCs) treated with F. gigantica ESPs. A total of 1,544 differently expressed mRNAs (790 upregulated and 754 downregulated genes), 30 differently expressed miRNAs (24 upregulated and 6 downregulated genes), 136 differently expressed circRNAs (83 upregulated and 53 downregulated genes), and 1,194 differently expressed lncRNAs (215 upregulated and 979 downregulated genes) were identified. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses revealed that F. gigantica ESPs altered the expression of genes associated with the host immune response, receptor signaling, disease and metabolism. Results from RNA-seq were validated by qRT-PCR. These findings provide an important resource for future investigation of the role of mRNAs and non-coding RNAs in mediating the immune-modulating effects of F. gigantica ESPs.

Fasciola gigantica produces excretory-secretory products (ESPs) with immunemodulating effects to promote its own survival. In this study, we performed RNA-seq to gain a comprehensive global understanding of changes in the expression of mRNAs, miRNAs, lncRNAs, and circRNAs in goat peripheral blood mononuclear cells (PBMCs) treated with F. gigantica ESPs. A total of 1,544 differently expressed mRNAs (790 upregulated and 754 downregulated genes), 30 differently expressed miRNAs (24 upregulated and 6 downregulated genes), 136 differently expressed circRNAs (83 upregulated and 53 downregulated genes), and 1,194 differently expressed lncRNAs (215 upregulated and 979 downregulated genes) were identified. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses revealed that F. gigantica ESPs altered the expression of genes associated with the host immune response, receptor signaling, disease and metabolism. Results from RNA-seq were validated by qRT-PCR. These findings provide an important resource for future INTRODUCTION Fasciolosis, caused by Fasciola gigantica and F. hepatica, is a significant parasitic disease, causing huge economic losses to livestock and serious adverse impact on public health (1)(2)(3). Fasciolosis is prevalent in developing countries, particularly in Asia, Africa, Latin America, and the Caribbean (2,4,5). The annual global economic losses to livestock production caused by liver fluke infection exceed $3 billion (6). Despite the significant impact of fasciolosis, our knowledge of host defense mechanisms against F. gigantica infection remains limited. Exploring the molecular mechanisms of interaction between the host and liver flukes would enhance our understanding of the host immune response and assist in identifying potential targets for the development of therapeutic interventions against fasciolosis.
Liver flukes employ a raft of strategies to subvert the host's immune response, including the secretion of excretory-secretory products (ESPs) to modulate the host's immune responses, to facilitate their own survival and to establish long-term infection (7,8). In a previous study, by using co-immunoprecipitation coupled with tandem mass spectrometry, we identified 14, 16, 9, and 9 proteins in F. hepatica ESPs (FhESPs) that can bind to cytokines IL2, IL17, IFN-g, and goat peripheral blood mononuclear cells (PBMCs) (9). In addition, a proteomic analysis of FgESPs interacting with buffalo serum at different infection periods (42,70, and 98 dpi) identified 18 proteins after F. gigantica infection, including 13 cathepsin proteins, 4 glutathione S-transferase proteins, and 1 calcium-binding protein (10).
Although ESPs play a significant role in supporting liver flukes' evasion of the host immune response, the molecular mechanisms that are mediated by ESPs remain poorly defined. To reveal new aspects of the immuno-modulatory effects of F. gigantica ESPs, we investigated the expression profiles of messenger RNAs (mRNAs), microRNAs (miRNAs), circular RNAs (circRNAs), and long non-coding RNA (lncRNAs) in goat PBMCs following incubation with ESPs of F. gigantica. The results showed that ESPs can regulate the proliferation, immune response, and metabolism of treated PBMCs. Our data provide new information for further in-depth analysis of the role of ESPs in the immune-evasive strategy of F. gigantica.

Source and Culture Conditions of the PBMCs
Goat PBMCs were isolated as previously described (11). Briefly, healthy goat venous blood samples mixed with ethylene diamine tetraacetic acid (EDTA) were collected into vacutainer tubes. Goat PBMCs were isolated using a commercial kit (TBD, Tianjin, China) according to the manufacturer's instructions. The isolated PBMCs were cultured in RPMI 1640 medium containing 10% fetal bovine serum (FBS) and 1% penicillinstreptomycin (Gibco, California, USA) and were incubated in a humidified atmosphere of 5% CO 2 at 37°C for 4 h.

Preparation of ESPs
ESPs of F. gigantica were prepared according to a previously described method (15,(17)(18)(19). Briefly, adult liver flukes of F. gigantica were obtained from the bile ducts of slaughtered buffaloes naturally infected by F. gigantica and washed three times with phosphate-buffered saline (PBS) at 37°C. The washed adult flukes were then transferred to PBSG buffer (PBS containing glucose, supplemented with 10,000 UI/ml penicillin G, 10 mg/ml streptomycin sulfate, and 25 mg/ml amphotericin B) and incubated in plastic petri dishes (1 worm/2 mL) for 2 h at 37°C. The supernatant was collected three times (i.e., every 2 h over a 6-h period). After incubation, the supernatant was collected and centrifuged at 10,000g for 30 min at 4°C, then filter-sterilized with a 20-mm nylon filter. After filter-sterilization, ESPs were concentrated using Amicon ® Ultra-15 centrifugal filter devices (30 K), and the concentration of ESPs was measured by the Bicinchoninic Acid Protein Assay (BCA-1, Sigma-Aldrich, Corp., St. Louis, MO, USA) and stored at −80°C until analysis.

RNA Extraction and RNA-seq
PBMCs (1 × 10 6 ) were incubated with 80 µg of F. gigantica ESPs or an equal volume of PBS (controls) at 37°C with 5% CO 2 for 24 h. This experiment was performed in triplicate. Total RNA of goat PBMCs was extracted using the Ribo-Zero Kit (Ambion, Austin, Tex, USA) following the manufacturer's protocol. RNA integrity was evaluated using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Samples with a RNA Integrity Number (RIN) ≥ 7 were included in the analysis. The sequencing libraries were constructed using TruSeq Stranded Total RNA with Ribo-Zero Globin, according to the manufacturer's instructions. These libraries were sequenced on an Illumina sequencing platform (HiSeq 2500).

Data Processing and Differential Expression Analysis
Raw reads generated from high-throughput sequencing were fastq format sequences. In order to obtain high-quality reads, trimmomatic software was used to remove low-quality bases, Nbases, and low-quality reads (20). HISTA2 was used to align clean reads to the goat reference genome (21). Also, Q20, Q30, and GC-content of the raw data were calculated (22). The "estimateSizeFactors" function of the DESeq R package was used to normalize the read counts (23), and the "nbinomTest" function was used to calculate P values and fold change values. The transcript with P-values < 0.05 and |log 2 Fold change| > 1 was deemed to be a differentially expressed (DE) transcript.

miRNA Prediction
The small RNA tags were mapped to a reference sequence by Bowtie (24) to analyze their expression and distribution on the reference genome. The data of miRBase20.0 were used as miRNA reference, and modified software miRDeep2 (25) was used to identify the known miRNAs. The software miREvo (26) and miRDeep2 (25) were combined to predict novel miRNA by exploring the secondary structure, the Dicer cleavage site and the minimum free energy of the small RNA tags unannotated in the former steps.

CircRNA Prediction
In order to generate a SAM file, BWA software was used to align the sequencing reads of each sample to the reference genome (27). Then, the CIRI software was used to scan paired chiastic clipping (PCC) signals, and circRNA sequences were predicted based on junction reads and GT-AG cleavage signals (28).

lncRNA Prediction
Stringtie software was used to assemble the known and new transcripts (29). The candidate lncRNA transcripts were screened by comparing with the gene annotation information of the reference sequence produced by Cuffcompare software (30). The transcripts with coding potential were then screened out using CPC (31), CNCI, Pfam (32), and PLEK (27) in order to obtain lncRNA predicted sequences.

GO and KEGG Enrichment Analysis
Gene Ontology (GO) enrichment was used for functional annotation analysis. GOseq based Wallenius non-central hyper-geometric distribution (33), which adjusts for gene length bias, was implemented for GO enrichment analysis. The Kyoto Encyclopedia of Genes and Genomes (KEGG) (34) database (http://www.genome.jp/kegg/) was utilized to test the statistical enrichment of the target genes in KEGG pathways using KOBAS software (35).

Quantitative Real-Time PCR (qRT-PCR)
qRT-PCR was performed to validate the transcriptome data using the SYBR assay according to the manufacturer's instructions. The qRT-PCR reactions were carried out using TB Green ® Premix Ex Taq ™ II (Tli RNaseH Plus) (Takara, Shuzo, Kyoto, Japan) with the following reaction conditions: 95°C for 30 s, followed by 40 cycles at 95°C for 5 s and 60°C for 30 s. The GADPH gene was used for mRNA, circRNA, and lncRNA normalization, while the U6 gene was used for miRNA's normalization. Three independent replicates were examined for each gene and the relative expression of the target mRNAs, mircRNAs, circRNAs, and lncRNAs was calculated using the 2 −DDCT method. All primers used in this study are listed in Table 1.

Cell Proliferation Assay by CCK-8 Kit
The effect of F. gigantica ESPs on goat PBMC's proliferation was investigated using the CCK-8 kit (Beyotime, Haimen, Jiangsu, China) as previously described (13). Briefly, goat PBMCs were incubated with F. gigantica ESPs at concentrations of 5, 10, 20, 40, and 80 µg/ml. Forty-eight hours later, 10 µl of CCK-8 reagent were added to each cell culture and incubated for 4 h. The OD 450 TABLE 1 | qRT-PCR primers used in this study.
of each well was measured using a microplate reader (Bio-Rad, California, USA).

Statistical Analysis
Student's t-test was used to analyze the data from qRT-PCR experiments using GraphPad Prism 5 software (GraphPad Software, Inc., San Diego, CA, USA). All data are presented as means ± standard deviations (SDs).

RESULTS
ESPs of F. gigantica Promoted the Proliferation of PBMCs F. gigantica ESPs promoted the proliferation of goat PBMCs in a concentration-dependent manner ( Figure 1). Significant differences were detected between cells treated with 20 (p < 0.1), 40 (p < 0.0001) and 80 (p < 0.0001) µg/ml ESPs, when compared with untreated (control) cells. However, no significant differences were detected between cells exposed to lower concentrations (5 or 10 µg/ml) and control cells.

DE mRNAs, miRNAs, circRNAs, and lncRNAs
The volcano plots ( Figure 2) and heat maps ( Figure 3) clearly showed the number and hierarchical clustering of differentiated mRNAs (Figures 2A and 3A), miRNAs ( Figures 2B and 3B), circRNAs ( Figures 2C and 3C), and lncRNAs ( Figures 2D and  3D) in ESP-treated and control (untreated) groups. Our analysis identified 1,544 differently expressed mRNAs (including 790 upregulated and 754 downregulated, Figure  2A and Supplementary Table S1), 30 differently expressed miRNAs (including 24 upregulated and 6 downregulated, Figure 2B and Supplementary Table S2), 136 differently expressed circRNAs (including 83 upregulated and 53 downregulated, Figure 2C and Supplementary Table S3), and 1,194 differently expressed lncRNAs (including  10 transformed P values (y-axis) are plotted against the average log 2 fold changes in expression (x-axis). Data points representing RNAs that were not identified as differentially expressed are shown in blue. Transcripts that are differentially expressed after treatment with ESPs (corrected P value < 0.05) with an absolute log 2 fold change (|FC|) greater than or less than 1 are shown as red (upregulated) and green (downregulated) dots. 215 upregulated and 979 downregulated, Figure 2D and Supplementary Table S4).

Validation of DE mRNAs, miRNAs, circRNAs, and lncRNAs by qRT-PCR
To validate the transcriptome data, a total of 12 genes were randomly selected for qRT-PCR verification, including three mRNAs (HDAC1, MCM2, CD47) ( Figure 4A Figure 4D). The qRT-PCR analysis showed that the trend of the data of qRT-PCR expression analysis was consistent with that of the data obtained by transcriptome analysis, supporting the accuracy of the transcriptomic data.  Figure 5A). Among these, the upregulated mRNAs were also significantly enriched in cellular component including the intracellular, intracellular part and organelles (Supplementary Figure S1A), and the downregulated mRNAs were significantly enriched in biological processes, such as the immune system process, regulation of localization, and regulation of the immune system process (Supplementary Figure S1B). The KEGG pathway enrichment was used to identify the significantly enriched pathways of the DE mRNAs. The top 20 pathways are shown in Figure 5B. Systemic lupus erythematosus, alcoholism, DNA replication and cell cycle were the most enriched pathways. Among them, the upregulated mRNAs were also significantly enriched in Systemic lupus erythematosus, alcoholism, DNA replication and cell cycle (Supplementary Figure S1C). The downregulated mRNAs were significantly enriched in Staphylococcus aureus infection, lysosome, complement and coagulation cascades, and phagosome (Supplementary Figure S1D).

GO and KEGG Analysis of the Target Genes of miRNAs
GO and KEGG pathway analysis was used to probe the main functions and the involved pathways of the DE miRNAs. As shown in Figure 6, the DE miRNAs were significantly enriched in biological processes, such as cell differentiation, cellular developmental process, and regulation of nervous system development. The most enriched pathway of the DE miRNA was the TGF-beta signaling pathway.

GO and KEGG Analysis of the Target Genes of circRNAs
The DE circRNAs were significantly enriched in biological processes, such as positive regulation of cellular process, positive regulation of biological process, and regulation of cellular metabolic process ( Figure 7A). The most enriched pathways of the DE circRNAs involved included lysine degradation, the TNF signaling pathway, the B cell receptor signaling pathway, and hepatitis C ( Figure 7B).

GO and KEGG Analysis of the Target Genes of lncRNAs
The biological functions of the DE lncRNAs were mainly involved in biological processes, such as negative regulation of biological process, organelle organization and macromolecular complex subunit organization ( Figure 8A). Among them, the upregulated lncRNAs were also significantly enriched in biological process and cellular component, such as macromolecular complex subunit organization, chromosome and chromosome organization (Supplementary Figure S2A). The downregulated lncRNAs were significantly enriched in the categories biological process and molecular function, such as enzyme binding, kinase binding and single-organism cellular localization (Supplementary Figure S2B). The DE lncRNAs were significantly enriched in Type I diabetes mellitus, Influenza A and Allograft rejection ( Figure 8B). Among them, the upregulated DE lncRNA were also significantly enriched in Influenza A, Type I diabetes mellitus and Graft-versus-host disease (Supplementary Figure S2C). The downregulated DE lncRNAs were significantly enriched in measles, rheumatoid arthritis, influenza A, and the chemokine signaling pathway (Supplementary Figure S2D).

DISCUSSION
In this study, we performed global RNA-seq to profile the expression of mRNAs, miRNAs, lncRNAs, and circRNAs in goat PBMCs exposed to F. gigantica ESPs. GO and KEGG enrichment analysis revealed that F. gigantica ESPs can influence the expression of transcripts associated with the host  immune response, receptor signaling, cell proliferation, disease and metabolism. In previous studies, we have shown that the F. gigantica proteins, rFgRab10, rFgCatB and rFg14-3-34, influence cell proliferation and induce apoptosis in goat PBMCs (11)(12)(13).
In the present study, we found that ESPs of F. gigantica promote cell proliferation as revealed by the CCK-8 assay ( Figure 1) and upregulate cell proliferation-related genes ( Table 2) of goat PBMCs. A previous study revealed that miR-211 suppresses epithelial ovarian cancer proliferation and cell-cycle progression by targeting Cyclin D1 and CDK6 (36), and miR-204-5p suppresses proliferation and invasion, and promotes chemotherapeutic sensitivity of colorectal cancer cells by downregulating RAB22A (37). The reduction in the expression of miR-211 and miR-204-5p in the present study implies a possible involvement of F. gigantica ESPs in promoting the proliferation of PBMCs. TLR3 and TLR4 play a critical role in innate immunity and induction of inflammatory response (38). Previous studies have shown that parasite infections, such as F. hepatica, F. gigantica, Leishmania, and Trypanosoma (39)(40)(41)(42), can modulate the expression of TLRs in order to evade the host's innate immune response. For example, stimulation of human monocyte-derived DCs with live Brugia malayi microfilariae significantly downregulates mRNA expression of TLR3, TLR4, TLR5 and TLR7 (43). In cattle, F. hepatica ESP decreases TLR2/TLR4mediated immune activation in order to limit excessive inflammation and attenuate the pathological reactions that accompany F. hepatica infection (44). Another study suggested that TLR4 is involved in F. hepatica fatty acid-binding protein (FhFABP) and induces alternative activation of human macrophages (45). FhFABP was suggested to suppress TLR4 activation in macrophages and inhibit the inflammatory cytokines induced following incubation with LPS (16). In the present study, PBMCs incubated with F. gigantica ESPs significantly downregulated mRNA expression of TLR3 and TLR4 ( Table 2). TLR3 suppression promotes a Th2 immune response to Schistosoma mansoni infection. In S. mansoni infected TLR3 gene-deficient mice, Th1-associated cytokines (e.g., IL-12) are significantly decreased and Th2-associated cytokines and chemokines are significantly increased (46). rFgCatB protein (component of F. gigantica ESPs) can induce the expression of Th2-associated cytokines (e.g., IL-4, IL-10, and TGF-b) (13). These data suggest that reduced expression of TLR3 skewed the host immune response towards a Th2 response during F. gigantica infection.
Interferon regulatory factor 3 (IRF3) is an important transcription factor involved in the downstream regulation of My88-independent TRIF-dependent TLR signaling (46), which regulates the production of type I IFNs and the induction of inducible nitric oxide synthase (iNOS) expression (47,48). iNOS catalyzes the conversion of arginine into citrulline and nitric oxide (NO), which protect the host against pathogen invasion (49). In this study, the expression of IRF3 was significantly decreased in PBMCs treated with F. gigantica ESPs. The reduced expression of IRF3 supports the link between inhibition of TRIF-dependent TLR3 signaling and decreased expression of NOS2 during F. hepatica infection (50,51), which may be one mechanism employed by the liver fluke to avoid damage caused by host-generated NO. These findings suggest a possible mechanism for F. gigantica ESPs in suppressing TLR expression or its signaling cascade.
Notably, in the present study, the NF-kB signaling pathway was suppressed. Therefore, the upregulation of miR-148a-5p detected in our study may suggest a role for F. gigantica ESPs in shifting Th1/Th2 towards the Th2 immune response.
Apoptosis of immune cells is an immunosuppressive strategy employed by F. hepatica, wherein the ESPs induce apoptosis in eosinophils via a caspase-dependent mechanism, resulting ROS- mediated mitochondrial-membrane potential loss (54). F. gigantica infection can induce apoptosis in sheep PMBCs via two distinct pathways: the extrinsic death receptor pathway, which is induced by the upregulation of TNF and TNFR1/ TNFR2 (40,50), and the intrinsic mitochondrial pathway (55). In our study, most death receptors and FAS receptor were downregulated, whereas the p53 signaling pathway-related genes were upregulated ( Table 2), suggesting that F. gigantica ESPs may induce apoptosis through the p53 signaling pathway, rather than by the Fas/Fasl model. Further investigation of the role of p53 signaling pathway in F. gigantica-induced apoptosis of immune cells is warranted. Fasciola infection can disrupt liver functions, resulting in disruption of the overall metabolism of infected livestock and losses in productivity (56). KEGG analysis revealed that exposure of goat PBMCs to F. gigantica ESPs upregulated genes related to lipid, nucleotide, amino acid, carbohydrate, and glycan metabolism. Alterations in metabolism during infection with F. gigantica and F. hepatica have been reported in previous studies (57)(58)(59). However, the mechanisms by which F. gigantica ESPs alter the host metabolic pathways remain to be elucidated.
Noncoding RNAs (miRNAs, circRNAs, and lncRNAs) participate in the regulation of different biological processes. In this study, a large number of circRNAs and lncRNAs were identified. However, the exact functions of these ncRNAs remain unclear. Pathways analysis indicated that the target genes of lncRNAs were highly enriched in disease-related pathways, such as toxoplasmosis, influenza A, and hepatitis B, and immune-related pathways, such as antigen processing and presentation, the Jak-STAT signaling pathway, and cytokine-cytokine receptor interaction. The target genes of circRNAs were highly enriched in the TNF signaling pathway and Hepatitis C. However, the target genes of miRNAs were only enriched in the TGF-beta signaling pathway. The immune-regulatory roles of these ncRNAs during F. gigantica infection need to be further investigated.
In conclusion, RNA-seq profiling of goat PBMCs treated with F. gigantica ESPs revealed significant changes in the expression of mRNAs, miRNAs, lncRNAs and circRNAs. The abovementioned F. gigantica ESP-induced changes appear to mediate a number of immune-regulatory mechanisms that enable F. gigantica to establish a persistent infection. These data should form the basis for further in-depth investigations of ESP-mediated immune-evasion strategies used by F. gigantica to subvert host immune responses.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The data is deposited in the NCBI repository, accession number is PRJNA686179.

ETHICS STATEMENT
This study was reviewed and approved by the Animal Ethics Committee of Lanzhou Veterinary Research Institute, Chinese Academy of Agricultural Sciences (permit no. 2018-012).

AUTHOR CONTRIBUTIONS
X-QZ, DC, S-SW, and HE conceived and designed the experiments. S-SW and DC performed the experiments. J-JH, W-BZ, A-LT, and G-HZ contributed reagents/materials/analysis tools. S-SW and DC analyzed the data and wrote the paper. HE, J-JH, W-BZ, and X-QZ critically revised the manuscript. All authors read and approved the final version of the manuscript.

ACKNOWLEDGMENTS
The authors would like to thank Novogene Bioinformatics Technology Co., Ltd (Beijing, China) for performing the sequencing and preliminary data analysis. The funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.