Fasciola hepatica Infection in Cattle: Analyzing Responses of Peripheral Blood Mononuclear Cells (PBMC) Using a Transcriptomics Approach

The parasitic helminth Fasciola hepatica (liver fluke) causes economic loss to the livestock industry globally and also causes zoonotic disease. New control strategies such as vaccines are urgently needed, due to the rise of drug resistance in parasite populations. Vaccine development requires a comprehensive understanding of the immunological events during infection. Previous in vivo studies by our group have investigated global differentially expressed genes (DEGs) in ovine peripheral blood mononuclear cells (PBMC) in response to both acute and chronic F. hepatica infection. This work demonstrated that pathways involved in the pathogenesis of ovine fasciolosis included fibrosis, inhibition of macrophage nitric oxide production, and antibody isotype switching, among others. Transcriptomic changes in PBMC populations following F. hepatica infection in cattle, in which the disease phenotype is quite different, have not yet been examined. Using RNA sequencing we investigated gene expression changes in PBMC isolated from 9 non-infected and 11 F. hepatica-experimentally-infected calves immediately before infection, at 1 and at 14 weeks post-infection. Longitudinal time-course comparisons between groups revealed 21 and 1,624 DEGs driven exclusively by F. hepatica infection in cattle at acute and chronic stages, respectively. These results show that fewer DEGs at the acute stage of infection can be identified in cattle, as compared with sheep. In addition, the log2 fold-changes of these DEGs were relatively low (−1 to 3) reflecting the different clinical presentation of F. hepatica infection in cattle. Gene pathways for hepatic fibrosis and hepatic cholestasis along with apoptosis of antigen-presenting cells were enriched at chronic stages. Our results reflect the major differences in the disease phenotype between cattle and sheep and may indicate pathways to target in vaccine development.

The parasitic helminth Fasciola hepatica (liver fluke) causes economic loss to the livestock industry globally and also causes zoonotic disease. New control strategies such as vaccines are urgently needed, due to the rise of drug resistance in parasite populations. Vaccine development requires a comprehensive understanding of the immunological events during infection. Previous in vivo studies by our group have investigated global differentially expressed genes (DEGs) in ovine peripheral blood mononuclear cells (PBMC) in response to both acute and chronic F. hepatica infection. This work demonstrated that pathways involved in the pathogenesis of ovine fasciolosis included fibrosis, inhibition of macrophage nitric oxide production, and antibody isotype switching, among others. Transcriptomic changes in PBMC populations following F. hepatica infection in cattle, in which the disease phenotype is quite different, have not yet been examined. Using RNA sequencing we investigated gene expression changes in PBMC isolated from 9 non-infected and 11 F. hepatica-experimentally-infected calves immediately before infection, at 1 and at 14 weeks post-infection. Longitudinal time-course comparisons between groups revealed 21 and 1,624 DEGs driven exclusively by F. hepatica infection in cattle at acute and chronic stages, respectively. These results show that fewer DEGs at the acute stage of infection can be identified in cattle, as compared with sheep. In addition, the log 2 fold-changes of these DEGs were relatively low (−1 to 3) reflecting the different clinical presentation of F. hepatica infection in cattle. Gene pathways for hepatic fibrosis and hepatic cholestasis along with apoptosis of antigen-presenting cells were enriched at chronic stages. Our results reflect the major differences in the disease phenotype between cattle and sheep and may indicate pathways to target in vaccine development.

INTRODUCTION
The liver flukes Fasciola hepatica and Fasciola gigantica are parasitic trematodes that affect cattle, sheep and goats, worldwide causing significant economic loses to agriculture (1,2). These parasites also infect people, with most human fasciolosis cases concentrated in South America, Africa and Asia (3), although human fasciolosis has also been reported in other areas such as Turkey (4), Serbia (5), Denmark (6), and Germany (7). In addition, to ruminants and humans a wide range of other mammalian species can be infected.
Fasciola hepatica infection of the definitive host occurs after ingestion of metacercariae, the infective stage, dispersed on pasture. Once the metacercariae are in the small intestine, newly excysted juveniles (NEJs) hatch, and migrate via the peritoneum for 4-6 days (early stages of infection) until they reach the liver capsule. In the liver, juvenile flukes commence migrating and feeding throughout the liver parenchyma. In this phase, haemorrhagic tracts and also the early stages of tissue repair can be seen. Eight to ten weeks later, the flukes reach the bile ducts, where they become mature and commence egg-laying. For several decades, triclabendazole has been the drug of choice for the clearance of F. hepatica infection as it can target early immature, immature and mature flukes, unlike other flukicides. Nevertheless, the use of this drug is restricted for dairy cattle in certain countries. In addition, an increase in triclabendazole-resistant liver fluke populations is evident (8), and presents difficulties for the control of fasciolosis in farmed ruminants. Accordingly, several research projects are pursuing alternative prophylactic and therapeutic strategies including new active compounds (9) and vaccines (10). The latter are especially desirable as vaccines could reduce the use of flukicides, thus slowing the emergence of anthelmintic resistance in the future, and do not leave drug residues that present risks both to the environment and to consumers of animal products.
In order to develop an effective vaccine strategy, it is important to understand the key immunological features that the parasite triggers in the mammalian host. It has been shown that stimulation with liver fluke-derived extracts, such as excretory/secretory (ES) products and tegument, induce unbalanced Th1/Th2 responses (11)(12)(13), alternative activation of macrophages (14) and cell death (15). In addition, glycoconjugates derived from liver fluke extracts induce activation of peritoneal eosinophils (16) and apoptosis of macrophages (17). Transcriptional analysis revealed upregulated TGFB1 following F. hepatica infection in ovine liver (18) and peripheral blood mononuclear cells (PBMC) (19,20). This is consistent with the fibrosis seen in chronic infection in sheep, although there are indications that by 14 weeks postinfection this is limited via SMAD7 down-regulation (20). Sheep transcriptome studies also revealed differentially expressed genes (DEGs) involved in T-cell activation (19), activation of the death receptor signaling pathway and the inhibition of nitric oxide production in macrophages along with Toll-like receptor signaling pathways after F. hepatica infection (20,21).
The use of RNA sequencing (RNA-seq) has also been useful in clarifying the underlying mechanisms involved in the immunology of fasciolosis, and in particular, how immunoregulation induced by the parasite can be overcome. For example, it has been hypothesized that the IL-18 and IL-12 complex, two important up-stream regulators, are inhibited during F. hepatica infection of sheep, and this could be responsible for switching the humoral response toward the production of IgG1 antibodies rather than IgG2 (20), as is seen during infection (22). In line with this observation, it has been shown that the involvement of genes associated with nitric oxide production, IL-12 signaling and IL-8 signaling could confer protective immunity when mice were vaccinated with T-cell epitopes derived from F. hepatica cathepsin B (23).
Cattle, along with sheep, are the primary target for a liver fluke vaccine. The disease phenotype in cattle differs significantly from that in sheep. Whereas, cattle may exhibit delayed growth and loss of productivity due to chronic infection, sheep can suffer from sudden death following heavy infection (24), as well as acute disease. Accordingly, we set out to carry out transcriptomic analysis of PBMC during the acute and chronic stages of F. hepatica infection in cattle, and to compare this with similar studies previously performed in sheep.

Ethics, Experimental Design, and Group Size Calculation
Animals and procedures carried out were approved by the UCD Animal Research Committee (AREC-  and licensed by the Health Products Regulatory Authority (AE18982/P098). Twenty castrated male Holstein-Friesian and crossbred Holstein-Friesian calves, 3-6 months of age, obtained from a small number of local farms with no history of F. hepatica infection were included in the study. Animals were divided into uninfected (n = 9; 4 Holstein-Friesian and 5 crossbred Holstein-Friesian) and infected (n = 11; 6 Holstein-Friesian and 5 crossbred Holstein-Friesian) groups and housed at UCD Lyons Research Farm throughout the experimental period. Calculation of the minimum number of animals required was based on Hart et al. (25) considering the following factors: two-tailed test with 90% power, minimum depth of coverage of 10 million reads per sample, fold change in expression between groups equal to 2, and the coefficient of variance of gene expression as an average of the variance of gene expression seen in the bovine atlas (26). Inclusion of animals in one group or another was based on blocking by age and breed in order to reduce confounding.
Animals were housed in two adjacent pens from 3 weeks prior to infection (W-3) to 15 weeks (W15) post-infection. In order to reduce pen housing as a confounding factor, animals from both uninfected and infected groups were mixed and allocated to each pen, while each also contained 5 Holstein-Friesian and 5 crossbred Holstein-Friesian animals. Throughout the experimental period animals were fed silage and had access to water ad libitum. Health checks were conducted daily during the study.

Experimental Infection
Sedimentation performed on fecal samples for F. hepatica eggs, and ELISA testing for the absence of F. hepatica cathepsin L1 (CL1) antibodies in sera as previously described (22) were used to demonstrate that animals were free from F. hepatica infection at the beginning of the study (W-3). Three weeks later (W0), each animal in the infected group was orally challenged with 200 F. hepatica metacercariae, Italian strain (Ridgeway Research Ltd, UK) suspended in 10 mL of distilled water. Animals in the uninfected group received 10 mL distilled water only. The viability of metacercariae was verified by excystation tests as described elsewhere (27).

Sampling Procedures
Blood samples were collected by jugular venepuncture at the day of infection (W0), but before the animals were infected; then at 1 week (W1) and at 14 weeks (W14) post-infection; using plain, EDTA, and lithium heparin-coated vacutainers. Clotted samples were centrifuged at 2,300 × g for 5 min, and the serum subsequently removed and frozen at −20 • C until ELISA testing. Antibodies specific for F. hepatica CL1 were detected by ELISA as previously described (22). EDTA samples were used to obtain total and differential white cell counts using an ADVIA R 2120 hematology analyser (Siemens). Fecal samples were collected at W14 in order to detect fluke eggs, nematode eggs and lungworm larvae, using Sedimentation, McMaster and modified Baermann testing, respectively.

PBMC Isolation and RNA Extraction
Blood collected in lithium heparin tubes was used for PBMC isolation as previously described (20) with the substitution of 15 mL of Ficoll-Paque PLUS (GE Healthcare) for Histopaque 1.077 g/mL. Following evaluation of cell viability by trypan blue staining and cell counts, aliquots of 2 × 10 7 PBMC per sample were washed and suspended in 600 µL of RLT buffer (Qiagen) supplemented with 10 µL/mL of β-mercaptoethanol. Cells were disrupted immediately by vortexing for 30 s and stored at −80 • C. Total RNA was isolated using the RNeasy Plus Mini Kit (Qiagen) following the manufacturer's instructions. This kit includes a genomic DNA eliminator column. The quantity of RNA extracted was measured using a Nanodrop 1000 spectrophotometer (Thermo Fisher Scientific), and RNA quality was evaluated using an Agilent 2100 Bioanalyzer with an RNA 6000 Nano Lab Chip Kit (Agilent Technologies). All samples exhibited A260/280 ratio >2.0 and RNA integrity number (RIN) >8.

Sacrifice of Animals and Post-mortem Examination
At week 15, animals were sent to the abattoir, and livers were recovered. After gross examination of the liver parenchyma and the bile ducts the liver fluke burden was estimated for each animal as previously described (28). Mean fluke numbers collected from Holstein-Friesian and crossbred Holstein-Friesian calves of the infected group were compared with the Wilcoxon rank-sum test.

RNA Preparation and Sequencing
RNA-seq library preparation and sequencing was performed by the Brigham Young University (BYU) DNA Sequencing Center (https://dnasc.byu.edu, UT, USA). Sixty barcoded paired-end (PE) mRNA-seq libraries were prepared using a KAPA stranded mRNA-seq kit and KAPA Dual-Indexed Adapter Kit (Roche), following the manufacturer's instructions. After quality control, a pool of all 60 barcoded samples was loaded across 8 lanes of v4 chemistry high output flow cells at 2 × 125 nucleotide reads on an Illumina HiSeq 2500 Sequencing System. All RNA-seq data has been made available via the European Nucleotide Archive repository under the accession number PRJEB32022.

Differential Expression Analysis
Bioinformatics analyses were performed using scripts developed in GNU bash version 4.3.48 (http://ftp.gnu.org/gnu/bash), and R version 3.5.1 (29). All scripts can be accessed at a public GitHub repository (https://github.com/carolcorreia/Liver-Fluke-RNAseq). Deconvoluted FASTQ files were evaluated with FastQC v0.11.5 (www.bioinformatics.babraham.ac.uk/projects/ fastqc). After filtering out adapter sequences, removing poor quality reads (those where more than 25% of the nucleotides had a Phred Score (Q) lower than 20) and removing PE reads which lengths were shorter than 100 nucleotides with ngsShoRT v2.2 (30), the clean FASTQ files were re-assessed with FASTQC as before. Filtered PE reads were then aligned and annotated to the bovine genome (Bos taurus_ ARS-UCD1.2) using the alignment software STAR v2.6.0c (31). After the 2-pass alignment and elimination of unmapped and multimapping reads, the gene counts were generated using the package featureCounts from the Subreads v1.5.1 software package (32). In addition, for each sample, filtered reads were subjected to variant calling for the evaluation of SNPs. The identity of samples was re-confirmed by the identity-by-state (IBS) comparison with PLINK v1.9 (33). The Z2 values between two samples were ranked and the first 60 Z2 values were used to verify whether they belonged to the same animal.
Prior to normalization of the libraries, lowly-expressed genes with <1 count per million in 9 or more libraries were discarded from further analysis. Normalization factors were then calculated using the trimmed mean of M-values with the R/Bioconductor EdgeR package v3. 22.3 (34, 35). The overall data structure was reviewed with multidimensional scaling (MDS) at each time point to detect potential outliers. Two approaches were used for identification of DEGs ( Figure 1A) non-paired comparisons between uninfected and infected groups at each time point and ( Figure 1B) longitudinal time-course analysis of DEGs for each group and then selecting those obtained exclusively in the infected group. Model matrices were generated for each comparison. For the first approach, infection status and time point were modeled, including interaction and blocking by breed. For the second approach, separate design matrices were generated for infected vs. non-infected animals, examining interactions between infection status and time, whilst blocking by animal. For both approaches, common and trended dispersions were estimated using the Cox-Reid profile-adjusted likelihood function (36) and tagwise dispersion was estimated with the empirical Bayes method to fit negative binomial generalized linear models. DEGs were determined using the quasi-likelihood F-test (37) choosing a false discovery rate (FDR) cut-off of 5%.

Ingenuity Pathway Analysis of Differentially Expressed Genes
Genes identified as differentially expressed exclusively in the infected group using the second time course approach ( Figure 1B) were analyzed using Ingenuity Pathway Analysis (IPA; Qiagen Inc., www.qiagenbioinformatics.com/products/ ingenuitypathway-analysis). Human orthologs of bovine genes were searched and merged using the HGNC Comparison of Orthology Predictions (HCOP) (38). Where a bovine gene was annotated to multiple human genes just one human ortholog was used for the analysis. Detection of significant functional canonical pathways was established when the enrichment P-value was <0.05.
The IPA "upstream regulators" and "causal network" features were used to analyse upstream and master regulators. The Ingenuity Knowledge Base was used to predict the expected causal effects between upstream regulators and DEG targets. Ingenuity Knowledge Base was also used to predict the expected downstream effects due to the DEG targets by using the "disease and function" and "regulatory effect" functions. Hits were considered significant when the P-value of overlap was lower than 0.05. Z-score values > +2 or below −2 were used to predict the activation or inhibition status of the upstream and master regulator or of the increased or decreased representation of a disease or functions, as appropriate.

PCR Validation
Six genes involved in immunological or signaling pathways (IL17RB, IL1B, TNF, TGFB2, PRG3, and NOS2) were selected for validation using reverse transcription quantitative real-time PCR (RT-qPCR). For all 60 samples, cDNA was synthesized with 1.0 µg of total RNA from each sample used for the RNA-seq library preparation. The High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems) was used for this purpose according to the manufacturer's instructions. After cDNA synthesis the 20.0 µL reaction was diluted 1:20 with distilled H 2 O prior to use. Primers were designed based on the corresponding cDNA sequences (obtained from the Ensembl database), and primer efficiencies were determined using a 1:4 dilution series over 7 points. PCR amplification specificity was assessed using melting curve analysis to confirm the presence of a single sharp peak. RT-qPCR was performed using SYBR Green Master Mix (Applied Biosystems) on a 7500 Real-Time PCR System (Applied Biosystems). All samples were analyzed in duplicate. Each 20 µL PCR reaction contained 5.0 µL diluted cDNA samples (or appropriate controls), 10.0 µL of SYBR Master Mix and 300 nM final concentration of each primer pair. Each RT-qPCR reaction was performed using the following cycling parameters: 2 min at 95 • C (heat-activation step); 40 cycles of 15 s at 95 • C, 1 min at 60 • C. A dissociation curve was also included. The stability of eight potential reference genes (ACTB, GAPDH, H3F3A, PPIA, RPL19, RNF11, SDHA, and YWHAZ) was assessed using the GeNorm algorithm within the qBase+ software package (Biogazelle) (39,40). H3F3A and RPL19 were shown to be the most stably expressed reference genes (M<0.14) and were used for subsequent normalization. Calibrated normalized relative quantities were calculated using the qBase+ software. Two-tailed paired sample t-tests were used to identify statistically significant gene expression changes using the same software.

Confirmation of F. hepatica Infection by ELISA, Sedimentation, and Fluke Burden
Prior to the beginning of the experimental infection, all the animals were healthy and showed no signs of disease. Baermann tests were negative for all animals at W0 and W14 ruling out the presence of lungworm infection during the experimental period. Two animals, one from each group, were observed to have nematode eggs in their feces, with fecal egg counts (FECs) of 50 strongyle eggs per gram at W0. None of the animals showed evidence of anti-F. hepatica CL1 IgG1 antibodies or presence of liver fluke eggs at W0. After challenge with F. hepatica metacercariae, with viability and excystation confirmed by a 65% successful hatching rate in vitro, F. hepatica eggs were detected in nine out of eleven animals belonging to the infected group at W14. All 11 animals in the infected group seroconverted, as determined by ELISA, with O.D. values above the cut-off at W14 (Figure 2A).
When livers were recovered, all the infected animals were shown to harbor adult and immature F. hepatica specimens in bile ducts and liver parenchyma. Individual liver fluke counts are shown in Figure 2B. No significant differences were detected in fluke numbers between Holstein-Friesian and crossbred Holstein-Friesian animals of the infected group (Mann-Whitney U, P = 0.405). Hematological parameters were within the reference range for all the animals at W0. There were no significant differences in the numbers of lymphocytes, monocytes, and eosinophils between uninfected and infected groups at any time point. There was an increase in the eosinophil counts between W1 and W14 in both groups ( Figure 2C), but values remained within the reference range. Lymphocyte numbers in infected crossbred Holstein-Friesian at W1, and eosinophil numbers in both infected and control crossbred Holstein-Friesian animals at W14, were significantly greater than comparator cohorts ( Figure 2C). Infected animals showed thickening of bile ducts, fibrotic areas, and haemorrhagic tracts in the liver (Supplementary Figure S1).

Bioinformatics and Mapping to the Bovine Genome
Sequencing of the pooled RNA-seq libraries across the 8 lanes generated a mean of 29.94 million raw reads per library. From these, 87.72% (∼26.24 million reads) of the reads passed the filtering steps, and a mean of 21.82 million reads uniquely mapped to the bovine genome (Supplementary Table S1). Calculation of the IBS distance between samples confirmed the correct labeling of samples throughout the experiment. Z2 values calculated between samples from different individuals were lower than 0.43 (Supplementary Table S2). Post-alignment quality control showed that libraries from all samples displayed a similar density of gene counts (Supplementary Figure S2).
Following the two approaches used for identification of DEGs ( Figure 1A) non-paired comparisons between uninfected and infected groups at each time point and ( Figure 1B) longitudinal time-course analysis of DEGs for each group and then selecting those obtained exclusively in the infected group, different insights were obtained. MDS showed no separation between uninfected and infected animals after F. hepatica infection at any time point. As shown in the MDS plot (Figure 3), there were clear overlaps between animals from uninfected and infected groups at each of the time-points, suggesting that infection was not the largest source of variation among animals. When the MDS plot was re-analyzed including breed as an additional factor, there was a clear separation between Holstein-Friesian and Holstein-Friesian crossbred cattle on the x-axis (BCV distance 1). Indeed, there were 3,410 DEG identified between Holstein-Friesians and crossbred Holstein-Friesian at W0 (Supplementary Table S3). Because breed was considered in the experimental design, sub-groups of each breed were analyzed separately in the MDS plot. Again, there was no separation between groups at W1 and W14, suggesting again that F. hepatica infection might not have had the largest impact on gene expression within PBMC.

Longitudinal Time-Course Contrasts Conducted Between Time-Points in Uninfected and Infected Groups Showed the Largest Numbers of DEGs at W14
For the first approach (Figure 1A), there was no distinction of DEGs between uninfected and infected animals at any time-point. Using the second approach (Figure 1B), which overcomes individual variability, we identified 21, 3,992 and 2,324 DEGs at the W1 vs. W0, W14 vs. W0, and W14 vs. W1 comparisons, respectively. The list of these DEGs is provided in Supplementary Table S4. After performing the same comparison in the uninfected group and removing the DEGs that overlapped in both groups for each of the time points (Supplementary Figure S3), 21, 1,624, and 551 DEGs were found exclusively in the infected group at W1 vs. W0, W14 vs. W0, and W14 vs. W1, suggesting that these are DEGs driven exclusively by F. hepatica infection (Supplementary Table S5). RT-qPCR results for the six genes selected confirmed the significance and trend of expression (Supplementary Table S6). For all comparisons, the magnitude of change in gene expression was relatively low. The log 2 fold change (log 2 FC) values for the DEGs varied from −1.15 (FDR = 3.5E-04), for the nebulin gene (NEB) at the W14 vs. W0 comparison to +3.2 (FDR = 7.3E-05) for the proteoglycan 3-like gene (LOC781990) at the W14 vs. W1 comparison (Figure 4).
As shown in the Venn diagram ( Figure 5A) there were 4 genes upregulated and 8 genes downregulated common to animals       could not be predicted. There were 8 pathways activated and 25 inhibited at W14 compared to W0. Ten canonical pathways that were also detected as enriched at W1, including IL-6, IL-10, and Toll-like receptor signaling, were also enriched at this stage. However, new enriched canonical pathways such as hepatic cholestasis and hepatic fibrosis/hepatic stellate cell activation were also identified. Within the inhibited canonical pathways were CD40, apoptosis signaling, TNFR1 signaling, hepatic growth factor (HGF) signaling and death receptor signaling. The production of nitric oxide and reactive oxygen species in macrophages had a trend toward inhibition with a z-score value very close to the threshold (−0.83).

Activation of Apoptosis of Antigen-Presenting Cells Is One the Main Regulator Effects Identified in the Chronicity of F. hepatica Infection
When the presence and nature of upstream and master regulators was evaluated, that is, genes that can interfere with the transcription of other genes at different levels, IL-10 was the only upstream regulator whose z-score value (1.715) was predicted in the W1 vs. W0 comparison. It was possible to identify the activation status of 105 and 27 genes at W14 vs. W0 and W14 vs. W1 (Supplementary Table S8). The upstream regulators common to both comparisons were the activation status of IL-1B, APP, TNF, and TLR7 and the inhibition of dual specificity phosphatase 1 (DUSP1) ( Table 3). IL-10 was also predicted in both comparisons with an inhibitory trend but did not reach the z-score cut-off. In terms of prediction of regulatory effects, DEGs of the W14 vs. W1 comparison clustered within 5 primary groups ( Table 4). The cluster with the highest consistency score predicted that the inhibition of the master regulator DUSP1, the inhibition of the upstream regulator mir-155 and the activation of the upstream regulators STAT3 and APP lead to the downstream effects of activation of apoptosis of antigen-presenting cells, migration of phagocytes and mononuclear leukocytes, recruitment of myeloid cells and phagocytes and chemotaxis (Figure 6).

DISCUSSION
Due to the economic importance of F. hepatica infection in animal agriculture, and anthelmintic resistance concerns, there is an urgent need for the development of novel control methods such as vaccination. To date, inconsistent results have been obtained in ruminant vaccine trials with a wide repertoire of recombinant proteins (10). Even though considerable work has been undertaken to understand the immune response to  F. hepatica infection, along with the immunoregulatory effects induced in infected animals, there are still many gaps in our knowledge, which could be important for the design of new vaccine strategies. In this regard, transcriptomics using technologies such as RNA-seq, represents a key tool in bridging these knowledge gaps.
To the best of the authors' knowledge this is the first transcriptomics analysis conducted in cattle (Bos taurus) investigating the immunological changes in PBMC following F. hepatica infection. The immunological effects triggered by the parasite can be detected in PBMC as "echoes" of the events occurring locally in the peritoneal cavity and the liver at the acute and chronic phases of the infection, respectively.
In this study, the phenotypic hallmarks of hepatic fibrosis and cholestasis in the chronic phases of F. hepatica infection in cattle were reflected in the altered gene expression changes seen at W14 in the infected animals. Although no differences in the numbers of eosinophils in peripheral blood were seen between infected and uninfected groups, PRG3, whose transcription and translation leads to eosinophil major basic protein, was the second most highly upregulated gene observed at the W14 compared to W1. Up-regulation of PRG3 was also detected in ovine PBMC at 2 and 8 weeks post F. hepatica-infection (19) indicating a greater degree of eosinophil activation in infected animals.
It was somewhat surprising that the MDS plots in this study did not show a clear separation according to F. hepatica infection status. While sample identity verification tools have been widely used in humans and mice (41), no such tools have been developed for cattle as yet. Therefore, to validate that the ambiguous MDS plot results were not a consequence of sample mislabelling in our study SNPs derived from RNA-seq were used to calculate IBS distances and verify animal identity across time points. Animals were not exposed to F. hepatica previously and they only received one infective dose of metacercariae. It is known that small repeat infections, similar to those which occur naturally (12), produce more severe pathology than a single infective dose with equal numbers of metacercariae in cattle (42); therefore, stronger infection effect might have been observed on the MDS analysis if trickle infections had been used. MDS plots also suggest that there could be intrinsic confounders, such as breed (43), temporal changes affecting animal development, or additional environmental factors that could have had a higher impact on gene expression than infection status. This emphasizes the importance of including an uninfected time course control in these types of studies. Almost 50% of the DEGs identified at the time-course comparisons W14 vs. W0 and W14 vs. W1 were observed in both infected and uninfected animals, all of them with the same direction of expression. We consider that the use of animals aged 3-6 months of age, as a primary target group of any novel anti-fluke vaccine, was appropriate. However, it is possible that greater temporal variation in expression levels of genes related to growth and development in this cohort, as compared with adult animals, may be relevant to the interpretation of the results, and that decreases in γδ T-cell levels, for example, as animals age, may affect the immune response. As no concurrent diseases were detected and animals were housed and managed in identical conditions throughout the experimental period, we were able to identify a subset of DEGs due to F. hepatica infection. Particular attention must be attributed to breed as it is known that differences in animal susceptibility and immunological features vary across breeds. For example, differences in nitric oxide production by macrophages after in vitro stimulation of the Brown Swiss cattle breed compared to Holstein-Friesian have been demonstrated (44) and relative resistance of the Indonesian thin tail sheep breed against F. gigantica (45) has been documented. Changes in global gene expression involved in some biological pathways, including immune responses, from different cattle breeds in leukocytes (46) have also been demonstrated. In our study, 3,410 DEGs were identified between Holstein-Friesian and crossbred Holstein-Friesian cattle, and the MDS plot also confirmed breed as a significant driver of variation. However, this breed effect could be mitigated in the experimental design as equal numbers of animals of each breed per group were included. As there were no significant differences in the liver fluke burden between the two breeds (P = 0.405) no further evaluation of the DEGs between the Holstein-Friesian and crossbred Holstein-Friesian groups in response to infection was conducted. However, further studies should be designed, with more biological replicates per breed, in order to confirm whether breed should be considered as an important factor in investigating response to F. hepatica infection. This requirement is further highlighted by differences between the two breeds in terms of their lymphocyte eosinophil counts at specific timepoints.
The clinical consequences of fasciolosis differ markedly between cattle and sheep, with the latter suffering potential acute consequences, and mortality, whereas in cattle disease is usually chronic (2). In comparing the results of this study with those in a previous study in sheep (20), parallel differences between the transcriptomic responses in the two species are striking. It is not possible to conclude definitively that these differences do not arise, at least in part, from variation in experimental protocol (comparative level of challenge, for example). However, in neither case the level of challenge was sufficient to produce clinical signs, so we suggest that this evidence indicates two clear differences in the response to F. hepatica between cattle and sheep. The first key difference is the proportion of DEGs detected at acute vs. chronic phases with F. hepatica responsible for almost 5% of the total DEGs at the acute phase in cattle whereas in sheep this percentage increases to more than 70% for the same time period (19,20). In cattle, most of the DEGs were observed at the chronic stage of infection. This indicates that when immature flukes penetrate the gut and migrate through the peritoneum of cattle, relatively few changes in gene expression in PBMC are triggered. However, when flukes are migrating through the liver parenchyma and into the bile ducts, they induce a greater number of geneexpression changes. This accurately reflects the different clinical presentations of fasciolosis between the two species. Whereas, sheep are susceptible to acute disease that can be fatal, fasciolosis in cattle has a much more muted clinical presentation, where animals in most cases appear clinically-normal even though they may experience production losses (47) and bystander immunoregulation (12,48). The second key difference is the magnitude of the log 2 FC in response to infection. In cattle, log 2 FC in DEGs at the different time-points varied from −1 to 3 whereas in sheep log 2 FC up to 22 in PBMCs (20) and 14 in liver tissue (18) were observed. These results also reflect differences in the disease phenotype between these two host species. A study in water buffalo found DEGs with log 2 FC varying from −5 to 5 from liver tissue at chronic stages of F. gigantica infections (49), comparable to the magnitude of log 2 FC detected in the present study.
Apart from the two main differences between cattle and sheep outlined above, a combination of information from the DEGs, canonical pathways, upstream regulators, and regulatory effects detected pinpoint immunological events in bovine fasciolosis.

Polarized T-Helper and Immune Responses
This study shows that the Th2 response pathway is enriched after only 1 week of F. hepatica infection. Due to the low numbers of DEGs found at W1, no clear mechanisms related to T-cell responses can be detected at this acute stage of infection, but more information can be deduced from effects after chronic infection. For instance, upregulation of the IL6 and IL4I1 genes and involvement of IL-4 signaling pathways is characteristic of a skewed Th2 response (50) and this response is typically driven by chronic F. hepatica infection (48). The fact that the activation status of IL-10 as an upstream regulator changed from an activated trend at W1 toward an inhibitory trend at W14 suggests that during the acute phase the immune response in cattle is a more balanced Treg-Th2 response, with increased polarization toward a Th2 response during chronic stages.
The expression of the intercellular adhesion molecule 1 gene (ICAM1), which contributes to recruitment of circulating leukocytes to infection sites and stabilizes cell-target interactions (51), was also upregulated at W14. The same effect was observed in buffaloes but at the acute phase of the infection (49). This could reflect that at chronic stages, there is a need to recruit cells to the site of parasite development to clear the infection. The soluble molecules sICAM-1 have been seen to be increased in complicated cases of humans infected with F. gigantica (52).
The class II major histocompatibility complex (MHC-II) transactivator gene (CIITA) is downregulated at the chronic stage of the infection. The CIITA protein is known to be a positive regulator of MHC-II gene transcription (53). Consequently, downregulation of CIITA could explain the fact that gene expression of the MHC-II DO beta gene (BOLA_DOB) is downregulated at the same period (log 2 FC = −0.323). Therefore, antigenic components from F. hepatica could interfere with maturation of antigen presentation cells and in the cross-talk to CD4 T-cells in cattle. Suppression of MHC-II has been described as an effect of fasciolosis in murine peritoneal dendritic cells (54,55) and in water buffalo (49).
We observed upregulation of TLR2, TLR4, and downregulation of TLR10 at the chronic stages of infection. This is consistent with previous findings-for example, increased expression of TLR2 gene by PBMC from F. hepatica-infected cattle for 8 weeks (56), of TLR4 protein by the mucin-derived peptide from F. hepatica (57) and downregulation of TLR10 gene in ovine PBMC (20). Nevertheless, Toll-like receptor signaling had a trend of inhibition in the time period from the acute to chronic fasciolosis, suggesting that F. hepatica could somehow overcome and interfere in the downstream effect triggered by increased TLR2 and TLR4. The suppression of TLR4 and TLR functions by antigenic components such as one of the F. hepatica fatty acid binding proteins (58) and glycoconjugates (59) has previously been reported. We also observed the activation of TLR7 in the time period from the acute to chronic fasciolosis in cattle. Activation of TLR7 gene has been reported in chronic stages of F. gigantica infection in water buffalos (60). The TLR7 receptor protein has been reported to mediate early innate immune responses in malaria (61) and it can contribute to the proliferation and activation of B cells and promotion of Th2 responses.

Liver Fibrosis
Hepatic fibrosis and hepatic cholestasis are clear clinical features seen in chronic fasciolosis. Previous studies have shown that TGFB1 is one of the main molecules orchestrating hepatic fibrosis after F. hepatica infection in sheep (18,20,62) and in cattle (63). Here, we detected the downregulation of TGIF1 as early as 1 week post-infection. This gene encodes a repressor of SMAD2 which leads to a controlled cell proliferation in acute liver injury (64). The plasminogen activator inhibitor-1 (PAI-1), which encodes a molecule with high fibrotic activity (65), was detected in ovine fasciolosis (20), and upregulated after 14 weeks of infection. Upregulation of TGFB1, TGFB2, COL1A1, and COL3A1 were found at W14 post-infection but in both infected and uninfected groups, suggesting that there may be alternative mechanisms for hepatic fibrosis in liver fluke-infected cattle. We suggest that the "fibrosis of liver" downstream effect could have been initiated at the chronic stage by the upregulation of TNF and that TGFB1 would contribute to the fibrotic effect at more advanced stages. Further work comparing TNF and TGFB1 changes from W14 onwards of F. hepatica infection would clarify the fibrosis events in cattle. TNF leads to the upregulation of the TNF gene, and the IL6 (66-68), urokinase-type plasminogen activator (PLAU) (69), SERPINE1, TNFRSF1A, SOCS1, and cathepsin B (CTSB) genes responsible for the fibrosis downstream effect (Figure 7). In addition, TNF activates IL1B, which also contributes to the upregulation of TLR4 and IL6. Upregulation of IL1B in water buffalos after 10 weeks of F. gigantica infection has previously been observed (49). Upregulation of TNF is also consistent with previous findings in sheep PBMC (20) and liver tissue (18) but not with other in vitro studies following incubation of F. hepatica antigens with human cells (70), murine cells in vivo and in vitro (58) and the bovine macrophage cell line BOMA in vitro (71). Discrepancies in TNF production between studies could be due to the differences in cell type populations when evaluated in vitro and in vivo, differences in parasite susceptibility between hosts and the antigen composition employed in each case. The activation of hepatic stellate cells by IL6/STAT3 to promote liver fibrosis has been documented in F. hepatica infection in humans and mice (72) which is in line with our findings in terms of the activation status of the IL6 and STAT3 signaling pathways at W14.

Apoptosis
Fasciola hepatica infection is known to induce apoptosis in peritoneal leukocytes (73), as well as in liver eosinophils (74), and in PMBC (20,21) in sheep. In the case of ovine PBMC, the activation of the extrinsic apoptosis pathway by the upregulation of TNF and TNFR1/TNFR2 genes (20) and the potential association of the intrinsic apoptosis pathway with the BH3-only member genes BIK, BNIP1, and BLC2L11 (21) has been postulated. In this study, apoptosis, and death receptor signaling pathways trended toward inhibition at W14. A close analysis reveals that activation of apoptosis of antigen presenting cells is one of the clearest regulatory effects, seen in combination with the migration of mononuclear leukocytes and migration and recruitment of phagocytes, orchestrated by the inhibitory effect of DUSP1 and mir-155 and by the activation of STAT3 and APP regulators. Although activation of STAT3 is known to reduce apoptosis, it could contribute to the upregulation of IL1B in combination with APP, mir-155, and DUSP1. Although DUSP1 was not detected as a DEG in this study it did trend toward upregulation (Log 2 FC = 0.263, FDR = 2.94E-01). Moreover, DUSP4 and DUSP16, were also observed to be upregulated. Other members of the DUSP gene family were upregulated in ovine PMBC at 2 and 8 weeks post-infection, respectively (19). Upregulation of DUSP13 has been shown in a bovine macrophage cell line stimulated with ES products from two F. hepatica isolates (71). DUSP13 inactivates important cell signaling mediators (JNK, p38, and ERK kinases as well as the transcription factor AP-1), which may lead to reduced inflammation (75).

CONCLUSIONS
In conclusion, this is the first in-depth transcriptomic analysis conducted in cattle for the assessment of the immune response against F. hepatica infection. This study pinpoints significant differences in terms of numbers of DEGs and the magnitude of transcriptional responses between cattle and sheep at the acute and chronic stages of infection. These differences are consistent with the different presentations of infection in these two species. It also predicts the potential role of TNF, IL-6, and IL-1B, in the development of early hepatic fibrosis events where TGF-β is not the key player. Apoptosis of antigen presenting cells is orchestrated by mir-155 and the products of the STAT3, APP, and DUSP1 genes at the chronic stages, demonstrating novel mechanisms triggering apoptosis during F. hepatica infection in cattle. This analysis is of pivotal importance for the better understanding of the immune response against F. hepatica in cattle, and for the development of novel Fasciola vaccines, where overcoming parasite-immunoregulatory strategies will be key to success.

DATA AVAILABILITY
The datasets generated for this study can be found in the European Nucleotide Archive Repository|PRJEB32022.

AUTHOR CONTRIBUTIONS
AG-C, JB, and GM conceived and designed the experiments. AG-C, AN-L, and LG-C collected samples from animals and carried out PBMC isolation. AG-C performed lab work and analyzed data. JB supported lab work. CC, GF, and DM assisted in data analyses and bioinformatics. AG-C and GM contributed to the writing of the manuscript. GM, JB, CC, GF, and DM contributed to the editing of the manuscript critically for important intellectual content. All authors read and approved the final manuscript.