Review of: "Timing of transcriptomic peripheral blood mononuclear cell responses of sheep to Fasciola hepatica infection differs from those of cattle, reflecting different disease phenotypes" peripheral blood mononuclear cell responses of sheep to Fasciola hepatica infection from

Timing of transcriptomic differs those of cattle, reflecting different disease phenotypes comments This manuscript is about Timing of transcriptomic peripheral blood mononuclear cell responses of sheep to Fasciola hepatica infection differs from those of cattle, reflecting different disease phenotypes. The aim of this paper was to further investigate the transcriptomic response of ovine peripheral blood mononuclear cells (PBMC) to F. hepatica infection, and to elucidate the differences between ovine and bovine PBMC responses. Authors have employed variety bioinformatics tools for the analysis of transcriptome data and identified different immunological pathways. One could also suggest that future studies must involve the use of flowcytometry for characterization of immune cells and real time PCR for specific selected immune genes. A disclaimer is that as a reviewer I did not go into the details of the statistical methods employed due to limited expertise towards those fields of study. Above all else, I truly, recommend that the article be accepted provided some corrections and suggestions are made. I have made some suggested corrections, typographic errors, etc, are indicated in the text of the manuscript to improve the article quality. Abstract Well written and structured, for specific comments see the text. Abstract Infection with the zoonotic trematode Fasciola hepatica, common in many regions with a temperate climate, leads to delayed growth and loss of productivity in cattle, while infection in sheep can have more severe effects, potentially leading to death. Previous transcriptomic analyses revealed upregulation of TGFB1, cell death and Toll-like receptor signalling, T-cell activation, and inhibition of nitric oxide production in macrophages in response to infection. However, the differences between ovine and bovine responses have not yet been explored. The objective of this study was to further investigate the transcriptomic response of ovine peripheral blood mononuclear cells (PBMC) to F. hepatica infection, and to elucidate the differences between ovine and bovine PBMC responses. data including IL-10 signalling, Th2[TN1] pathway, and Th1 and Th2 activation were upregulated only in the chronic phase in sheep. We propose that the earlier activation of anti-inflammatory responses in cattle, as compared with sheep, may be related to the general absence of acute clinical signs in cattle. These findings offer scope for “smart vaccination[TN2] ” strategies for this important livestock parasite.


Introduction
Fasciola hepatica is a zoonotic trematode that affects livestock (sheep, cattle, goats) as well as humans.
Human cases of fasciolosis occur mainly in Latin America, Africa and Asia, with approximately 2.5 million people affected worldwide (Mas-Coma et al., 1999, WHO, 2006. Sporadic human cases also occur in Europe as a result of the consumption of raw vegetables such as parsley or watercress, grown on irrigated or flooded land (Rondelaud et al., 2000).
F. hepatica infection is an important animal health issue in ruminants in endemic areas (Torgerson and Claxton, 1999). In Ireland, F. hepatica seroprevalence is reported as 75.4% to 82% of dairy herds, and within-herd prevalence is estimated at ≤ 50% (Bloemhoff et al., 2015, Selemetas et al., 2015. In Northern Ireland, 63.15% cattle herds were positive for F. hepatica, with an overall animal level prevalence of 23.68% (Byrne et al., 2016). Horses can also become infected, with a 9.5% prevalence of F. hepatica antibodies reported in a post-mortem survey of horses in Ireland (Quigley et al., 2017). The cost of livestock infection includes productivity losses (e.g., meat, milk, wool) as well as reduced fertility and condemnation of livers as unfit for human consumption (Hurtrez-Boussès et al., 2001, Schweizer et al., 2005. While F. hepatica infection in cattle leads to delayed growth and loss of productivity, acute infection in sheep can have more severe effects. The severe disease presentation in sheep is linked to a rapid acute response and can lead to sudden death (Behm and Sangster, 1999). Anthelmintic resistance in F. hepatica populations is an increasingly important issue, with resistance to triclabendazole of particular concern (Martínez- Valladares et al., 2014, Kelley et al., 2016, as this is the only drug effective against early immature flukes. Hence, there has been significant attention given to developing a vaccine effective in protecting ruminants against fasciolosis. While promising vaccine trials have been conducted, repeatability and consistency has been an issue [TN3] (Molina-Hernández et al., 2015, McManus, 2020. The ruminant immune response to F. hepatica infection involves an increase in peripheral eosinophil numbers, alternative activation of macrophages, and suppression of Th1 responsiveness to the parasites themselves and to bystander antigens (Flynn et al., 2010, O'Neill et al., 2001. The transcriptomic response of peripheral blood mononuclear cells (PBMC) to F. hepatica in cattle has been characterised by apoptosis of antigen-presenting cells (APCs), liver fibrosis and a Th2 response, with TNF, IL1B, and DUSP1, APP, STAT3 and mir-155 as important upstream regulator genes leading to hepatic fibrosis and apoptosis of APCs or migration and chemotaxis of leukocytes (Garcia-Campos et al., 2019). A balanced Treg-Th2 response was present in the acute phase with increased polarization towards a Th2 response during the chronic phase of infection.
The transcriptomic response of ovine PBMC has also been characterised previously. The response at 2 and 8 weeks post-infection (wpi) was associated with upregulation of the TGF-β signalling pathway, the Qeios, CC-BY 4.0 · Review, July 30, 2021 Qeios ID: 3NP14M · https://doi.org/10.32388/3NP14M 3/27 complement system, chemokine signalling and T-cell activation. Groups of genes were identified which were important for the immune response to the parasite and included those coding for lectins, proinflammatory molecules from the S100 family and transmembrane glycoproteins from the CD300 family.
Early infection was characterised by positive activation of T-cell migration and leukocyte activation, while the in the later stage[TN4] lipoxin metabolic processes and Fc gamma R-mediated phagocytosis were among the upregulated pathways (Alvarez Rojas et al., 2016). Another study which looked at the ovine acute and chronic response at 1 and 14 wpi identified upregulation of TGFβ signalling, production of nitric oxide in macrophages, death receptor signalling and IL-17 signalling at 2 wpi, as well as Toll-like receptor (TLR) and p53 signalling in both the acute and chronic phase (Fu et al., 2016). Overall, these studies suggest a drive towards an anti-inflammatory response and tissue repair in F. hepatica infection in both sheep and cattle.
The objective of this study was to further investigate the pathways involved in the response of ovine PBMC to F. hepatica infection, and to compare those with bovine responses to aid in developing effective vaccine strategies.

Animal trial outline
The animal trial has been described previously by Naranjo-Lucena et al. (2021). Briefly, sixteen male Merino sheep were obtained from a liver fluke-free farm at 9 months of age. All animals were tested monthly for parasite eggs by faecal sedimentation. Prior to challenge, all animals were tested for serum IgG specific for FhCL1 by an enzyme-linked immunosorbent assay (ELISA). Animals were housed indoors (100 m2 covered and 100 m2 uncovered facilities) and fed with hay and pellets and given water ad libitum. Sheep were randomly allocated to the infected and control group (n = 8 in both groups). At week 0, eight sheep were orally infected with one dose of 120 F. hepatica metacercariae of the Ridgeway (South Gloucester, UK) strain. Animals were monitored for 16 weeks, with faecal egg counts measured weekly starting from 7 wpi, and haematology measurements taken at 0, 2, 7 and 16 wpi. PMBC were purified from blood samples taken at 0, 2 and 16 wpi [TN5] and used for the transcriptomic analysis (Fig. 1A

Overview of the animal infection protocol (A) and bioinformatics analysis (B). Sixteen male Merino-breed
sheep were randomly assigned to an infected or control group and orally infected with a single dose of liver fluke (Fasciola hepatica). Blood was then collected for isolation of peripheral blood mononuclear cells (PBMC). These cells were used to generate RNA sequencing data, at 0, 2 and 16 wpi. These data were then analysed for DE genes using two approaches: analysis 1 (infected vs control) -infected animals were compared to control animals at each time point; analysis 2 (longitudinal) -the infected and control animals were compared separately to their pre-infection time points. The DE gene lists of both groups were then compared, and the DE genes which occurred in the infected group but not in the control group at each time point were considered to be the potential DE genes in response to F. hepatica.
wpi -weeks post-infection. Image created using biorender.com

RNA extraction from peripheral blood mononuclear cells (PBMC)
Blood samples were collected in 10 ml EDTA vacutainers (BD) from the jugular vein. Under a flow hood, 30 [TN8] ml of whole blood per animal was added into Falcon tubes and centrifuged at 1900 × g for 30 min without brake at 4°C to separate the buffy coat, red blood cells and plasma. The buffy coat was aspirated, added to a 15 ml Falcon tube and diluted 1:1 with PBS. After resuspension, the diluted buffy coat was gently layered on the top of 4 ml Ficoll Histopaque (Sigma-Aldrich) and centrifuged at 1900 × g for 30 min without brake at 4°C. Then, the buffy coat (PBMC) formed in the interphase between Histopaque and medium was aspirated. PBMC were washed once with 10 ml sterile PBS-EDTA and centrifuged at 1100 × g for 10 min at 4°C. After discarding the supernatant, the cells were resuspended in 5 ml of sterile ammonium chloride potassium, ACK buffer (155 mM NH4Cl, 10 mM KHCO3 and 0.1 mM ethylenediaminetetraacetic acid, EDTA) for 2 min at room temperature to lyse the erythrocytes. Then, sterile PBS-EDTA was added up to 14 ml and centrifuged at 1100 × g for 10 min at 4°C. The supernatant was discarded, and 10 ml of sterile PBS were added. At this step, PBMC were counted using Trypan blue for viability assessment. After counting, tubes were centrifuged at 1100 × g for 10 min at 4°C and PBMC were resuspended in a minimum volume of sterile PBS and divided into cryotubes centrifuged at 900 × g for 10 min at 4°C. After the PBS was removed, 1 ml of RNAlater (Sigma-Aldrich) was added per cryotube. The tubes were kept at 4°C overnight and stored at −80 °C until further assays.
Cell pellets were stored in RNAlater at −80°C until RNA extraction. Prior to extracting RNA, 950 μl of sample in RNAlater was mixed with 950 μl PBS and centrifuged at 3000 × g for 5 min. RLT buffer (Qiagen) was then added to the pellet and RNA was extracted using an RNeasy Plus Mini kit (Qiagen) as per the manufacturer's instructions. RNA was eluted in 35 μl of RNAse-free water (Sigma).

Quality control, assembly and DE gene analysis
All the bioinformatics pipeline bash and R scripts used for computational analyses were deposited in a GitHub repository at https://github.com/DagmaraNiedziela/RNAseq_Fhepatica_ovine_PBMC. RStudio v1.3.1093 and R v4.0.3 were used for the analyses. Sequence quality was assessed using FastQC v0.11.8 (Andrews, 2018). Trimming was performed using fastp v0.19.7 (Chen et al., 2018) with default settings as well as enabling base correction for paired end data, overrepresented sequence analysis, paired-end adapter removal, a minimum trimmed read length of 30 bases and a minimum phred score of 20. On average, 1% of the reads were trimmed per sample. Trimmed fastq files were re-checked with FastQC to confirm sequence read quality. Trimmed sequences were mapped to the ovine genome (Rambouillet v1.0) using the STAR aligner v2.7.3a (Dobin et al., 2013), with gene counts generated simultaneously with the STAR software using the Ensembl annotation v101 of the Ovis aries reference genome and annotation.
Genes with read counts < 10 across all samples as well as non-protein coding genes were removed in a filtering step (Smedley et al., 2009) Animal as batch effect was not included in this model, as animals from each group were considered as their own individual controls. A difference between the gene lists from infected and control group DE genes at each time point was examined as a potential response to F. hepatica. The longitudinal analysis was performed for comparison purposes, as several previous transcriptomic studies have used this approach.

Pathway and gene correlation analysis
Significant (Padj. < 0.1, primary analysis, Padj. < 0.05, longitudinal analysis) DE genes were converted to their one-to-one human orthologs using AnnotationDbi in R, ordered by adjusted P value and analysed for enrichment of Gene Ontology (GO) terms, Kyoto Encyclopedia of Genes and Genomes (KEGG) and REACTOME pathways using gProfiler R package v0.7.0 (Reimand et al., 2007). When a group of DE genes Qeios, CC-BY 4.0 · Review, July 30, 2021 Qeios ID: 3NP14M · https://doi.org/10.32388/3NP14M 6/27 matched multiple pathways, the pathway that was represented by the highest number of DE genes was selected in a filtering step. Ingenuity Pathway Analysis (IPA; Qiagen) was also used to examine overrepresented pathways and predicted upstream and downstream transcriptional regulators.

Cell composition analysis
Transcripts per million (TPM) values were generated from raw gene counts using a custom R script. Ovine Gene IDs converted to one-to-one human orthologs were used for the analysis. To estimate the proportion of the different immune cell populations in each somatic cell sample, a quanTIseq algorithm which is based on expression of cell specific surface marker genes was used (Finotello et al., 2019).

Parasitological and Clinical Examinations
No animals showed evidence of F. hepatica  Table 1).

Data exploration
Raw read counts were used for DE gene analysis using DESeq2. After filtering for low counts (sum of reads from all samples < 10) and non-protein coding genes, 15,363 genes were retained for further analyses.
The expression of the 500 most variable genes was used for generating PCA plots within the DESeq2 software package. Samples were observed to separate by time rather than by group, with 2 wpi and 16 wpi forming distinct clusters ( Fig. 2A). The largest sample variation was observed at time 0, with samples dispersed across the PCA plot. This variation may reflect inter-animal differences being more pronounced at time 0 due to an acclimatisation period. For subsequent individual time points, the clearest separation Qeios, CC-BY 4.0 · Review, July 30, 2021 Qeios ID: 3NP14M · https://doi.org/10.32388/3NP14M 7/27 between the infected and control groups was observed at 2wpi (Fig. 2B). PCA plots generated by time point, with infected animals in red and control animals in blue.

Cell composition
Normalised Undefined cells did not differ between groups or time points and represented 55.3 ± 0.8% in the control group and 54.8 ± 3.1% in the infected group. These "other" undefined cells most likely represent CD4 and CD8 T [TN11] lymphocytes. Both CD4 and CD8A were present in the data set and not differentially expressed, however the CD8B gene was not found in the ovine genome annotation, which was likely the reason for CD8 T-cells not being identified by QuanTISeq. [TN12] Numbers of T regulatory cells appeared higher in the infected group than in the control group at 2 wpi; however, the difference was not significant.
Overall, lack of major differences in cell composition indicates that the response seen in this study is likely to be due to changes in gene expression by particular PBMC cell types, rather than by changes in predominant cell type.

Differentially expressed genes[TN13]
Analysis of infected animals versus control animals (Analysis 1) yielded 59 DE genes at 0 wpi, 453 DE genes at 2 wpi, and 2 DE genes at 16 wpi (Padj. < 0.1; Fig. 3A). Only 7 of the 0 wpi DE genes overlapped with the 2 wpi DE genes, and these were removed from further pathway analysis (Fig. 3B). The single DE gene which was common between 2 wpi and 16 wpi was HSPB8, which encodes the heat shock protein family B (small) member 8, and was downregulated at both time points.   signalling, death receptor signalling and protein kinase A signalling (Table 3). There was a concordance between IPA and KEGG[TN19] , because TLR signalling, role of RLRs in antiviral innate immunity, and protein ubiquitination were also among the significant IPA canonical pathways. Furthermore, examination of GO biological processes confirmed the findings of KEGG, REACTOME and IPA, notably since pathways related to stress, protein ubiquitination, apoptosis and cell death were overrepresented (Table 2).
Selected pathways overrepresented in response to F. hepatica at 2 wpi (infected vs control analysis), as identified by gProfileR.  The key receptor in this pathway, RIG-I, is encoded by the DDX58 gene, which was upregulated in this study at 2 wpi (Fig. 5C). RIG-I is responsible for the innate recognition of uncapped 5'-phosphate double-stranded or single-stranded RNA found in numerous viruses including influenza and Ebola (Loo and Gale Jr, 2011).
Another key gene that was upregulated in our study and encodes an important transcription factor is the TNF receptor associated factor 3 gene (TRAF3). Since TRAF3 is key to intracellular responses, molecules secreted by the parasite may undergo endocytosis and be detected while inside the cell. Therefore, the RLR signalling pathway could be involved in recognition of and response to extracellular vesicles of F. hepatica, which can contain RNA and miRNA and are secreted throughout the life cycle of the parasite, including newly excysted juveniles and adults (Sánchez-López et al., 2020, Fromm et al., 2017. TRAF3 activity in RLR signalling is further regulated by the action of Polo-like kinases (Vitour et al., 2009), which were also upregulated at 2 wpi (Fig. 5A). Overall, the upregulation of these key genes in the RLR signalling pathway would ultimately lead to an increase in interferon production. While no interferon genes were detected as upregulated in this study, the activation may occur at a later stage in the infection.
[TN22] 3.6.1.2 Toll-like receptor signalling TLR signalling was identified using both the KEGG and REACTOME pathway databases. The only TLR gene that was upregulated in this pathway was TLR9 (Fig. 5B), which is activated by unmethylated CpG sequences in DNA molecules due to cancer, infection, or tissue damage. TLR9 is the receptor involved in the response to oncogenic viruses, and in autoimmune disorders. The overexpression of TLR9 coinciding with overexpression of DDX58 (RIG-I) further suggests that a response to nucleic acids secreted by F.
hepatica is taking place. TLR7 and TLR9 are selectively expressed by plasmacytoid dendritic cells (pDCs; also known as IFN-producing cells), which are a subset of DCs with a plasmacytoid morphology and unique in their capacity to rapidly secrete large quantities of type I IFN in response to viral infection, including response to influenza (Kawai and Akira, 2008). Glycans produced by F. hepatica were previously reported to induce a semi-mature state of DCs in mice, with reduced IFN-γ production (Rodríguez et al., 2016).
Excretory-secretory products (ESP) from F. hepatica have also been found to inhibit the production of effector molecules such as IL-6 and IL-12p70 in DCs following stimulation of TLR9 with CpG (Falcón et al., 2010). The DC receptor genes CD80 and CD83 were among the most upregulated DE genes in this study at 2 wpi (Table 1, Fig. 4A) -which may further suggest initial activation of DCs in response to F. hepatica, with further downstream response being inhibited by the parasite.
The lack of differential expression of other TLR genes in this study is in contrast to a previous study of transcriptomic response of ovine PBMC, where TLR1, TLR5, TLR6, TLR7 and TLR10 were downregulated in PBMC at the acute stage of infection (Fu et al., 2016). However, in a separate study of ovine PBMC, no TLR genes were DE at 2 wpi (Alvarez Rojas et al., 2016). TLR signalling also involves upregulation of TRAF3, which undergoes endocytosis together with TLR4 in response to lipopolysaccharide (LPS). LPS has been identified as a predicted upstream regulator by IPA at 2 wpi, which may indicate a potential influence of Ultimately, the TLR signalling pathway leads to upregulation of NF-kB and MAPK related pathways. In this study, the NF-kB inhibitor delta gene NFKBID was upregulated, and MAP kinases downregulated (Fig. 5B), both of which would lead to a decreased cytokine response. In fact, we did not detect any interleukin genes to be DE apart from IL16, which points to an inhibition of cytokine response taking place.
Caspase 8 is a common protein in both RIG-I and TLR signalling pathways and is known to be a key driver of cellular apoptosis. This caspase was among the most downregulated DE genes at 2 wpi, and suggests that an inhibition of cellular apoptosis was taking place at the acute stage of infection. Caspase 8 also may be key to the inflammatory response as casp8-/-mice exhibited a reduction in IL1B and IL12B expression, reduced recruitment of monocytes and an increased susceptibility to Toxoplasma gondii (DeLaney et al., 3.6.2 Several pathways related to apoptosis and cell death were identified 3.6.2.1 Transcriptional regulation by TP53 was overrepresented Transcriptional Regulation by TP53 is a REACTOME pathway that groups the upstream and downstream effects of the transcriptional regulator p53, most known for being involved in cancer regulation. The transcription factor recognizes specific responsive DNA elements and regulates the transcription of genes involved in cellular metabolism, survival, senescence, apoptosis and the DNA damage response Prives, 2009, Kruiswijk et al., 2015). PLK2, PLK3 and BTG2 were the most upregulated genes in the pathway (Fig. 5A), encoding the polo-like kinase (PLK) 2 and 3 and BTG anti-proliferation factor 2 proteins, respectively. PLKs play pivotal roles in cell cycle progression. PLK2 and 3 are activated by spindle checkpoints and DNA damage, and are required for entry into S phase (Pellegrino et al., 2010). The fact that the PLK2 and PLK3 genes are highly upregulated here would indicate a progression of the cell cycle and activation of cell proliferation. BTG2 belongs to the anti-proliferative (APRO) family of genes that regulate cell cycle progression by enhancing or inhibiting the activity of transcription factors. The BTG2 gene has a potential role in muscle fibre size, intramuscular fat deposition and weight loss in sheep and can lead to a decrease in cell proliferation or an increase in energy expenditure (Kashani et al., 2015).
Therefore, upregulation of this gene may be related to increased cellular metabolism.
Other DE genes at 2 wpi that may indicate perturbation of the cell cycle include the cyclin D3, L1, Q and T2 genes (CCND3, CCNL1, CCNQ, CCNT2) and cyclin-dependent kinase genes CDK17 and CDK5R1, with the majority of these being upregulated (Supplementary Table S2). The product of the CCND3 gene, downregulated in this study, is responsible for the progression through the G1 phase (Sarruf et al., 2005).
Cyclin types L and T, which were upregulated in this study, are involved in transcription (Chen et al., 2006, Herrmann et al., 2007. It is therefore likely that transcription is activated at 2 wpi rather than PBMC proliferation.

The role of PKR in interferon induction and anti-viral response pathway was inhibited at 2 wpi
Protein kinase receptor (PKR) belongs to a family of pattern-recognition receptors (PRRs especially the production of type I IFNs and increased expression of pro-apoptotic factors (Kang and Tang, 2012). The interferon alpha and beta receptor subunit 2 (IFNAR2) and interferon regulatory factor 9 (IRF9) genes were downregulated in this pathway (Supplementary Figure 2). IRF9 is a transcription factor which mediates signalling by type I IFNs (IFN-α and IFN-β). Following type I IFN binding to cell surface receptors such as IFNAR, STAT1 and STAT2 are phosphorylated and IRF9 associates with the phosphorylated STAT1:STAT2 dimer to form a complex termed ISGF3 transcription factor, which then enters the nucleus and activates the transcription of interferon stimulated genes thereby driving the antiviral response of the cell (Hernandez et al., 2018). Downregulation of these genes at 2 wpi, indicates that the interferon response is inhibited. The NFKB inhibitor delta gene (NFKBID) is upregulated, which suggests inhibition of the PKR signal transduction that occurs through the NF-κB pathway. MAPK genes are downregulated, which also suggests inhibition of the PKR signal transduction which occurs through the MAPK pathway. Overall, the downregulation of this pathway suggests that an inhibition of innate immune response and cellular apoptosis are taking place, therefore disabling key mechanisms that may facilitate host defence against F.
hepatica. Notably, antiviral response activation is suggested by upregulated genes in the TLR and RIG-I signalling pathways. This leads to the hypothesis that downstream interferon response is actively inhibited, even though foreign nucleic acids are being detected inside the cell.

Death receptor signalling was inhibited at 2 wpi
Death receptor signalling and apoptosis signalling were observed to be upregulated in the acute stage of F. hepatica infection in sheep previously (Fu et al., 2016). In this study, the IPA death receptor signalling pathway was inhibited[TN23] . Caspase-8 is a major initiator in the death signalling pathway (Ho and Hawkins, 2005). Previously, it has been shown that F. hepatica ESPs induced apoptosis of eosinophils, and caspase 3, 8 and 9 were activated in the process (Flynn et al., 2010). CASP8 and CASP9 were both downregulated in this study, and CASP8 downregulation has also been described previously in sheep (Fu et al., 2016). Overall, the death receptor signalling pathway was inhibited; however, the key genes involved in the pathway such as TNF and CASP3 were not DE. This pathway should therefore be studied in more detail in order to elucidate the mechanisms of cell survival or apoptosis during infection with F. hepatica.
Overall, these pathways suggest that immune response to exogenous nucleic acids as well as changes in energy metabolism and cell viability occurred in the acute phase of ovine F. hepatica infection.  Table S5). The number of DE genes at the acute stage of infection in cattle was much lower than during chronic infection, while in sheep more than 1000 DE genes were detected at both the acute and chronic phases. This suggests that although in cattle penetration of immature flukes from the gut to the peritoneum elicits a relatively muted host response, the ovine response at the acute phase is more marked.
DE gene lists from the current ovine study and the bovine study were converted to human orthologs and gene lists at the acute and chronic stage were compared (Fig. 6A). There were three common DE genes between cattle and sheep at the acute stage of infection: IL1RL1, CHD1 and RASSF1, and 232 common DE genes at the chronic stage of infection (Supplementary Table S6). When all time points in both species were compared, one DE gene was common to all time points in both sheep and cattle -the interleukin 1 receptor like 1 gene (IL1RL1). Common chronic DE genes included those involved in the action of TGFβ, such as the TGFB[TN24] induced factor homeobox 1 (TGIF1) and transforming growth factor beta receptor 3 (TGFBR3) genes. TGIF1 was downregulated in both species, while TGFBR3 was upregulated in cattle and downregulated in sheep. Pro-inflammatory cytokine and chemokine genes were also among the common chronic phase DE genes, including CCL18, CCL2, and IL1B, as well as the immune receptor genes CCRL2, IL10RA, IL2RG and TLR4. Several genes were common between the species, but exhibited opposing direction of expression (Fig. 6B). Pathway analysis of the longitudinal DE gene lists at acute and chronic phases was performed using IPA for both species. Three IPA canonical pathways were enriched in acute phases in both sheep and cattle, and eight canonical pathways were enriched in chronic phrases for both species. The common acute pathways between sheep and cattle included Toll-like receptor signalling, PPAR activation, and RAR activation (Fig.   7). In sheep two of these common pathways were overrepresented in both acute and chronic phases of infection. TREM1 signalling, STAT3 pathway and role of IL-17F in allergic inflammatory airway diseases were upregulated in both species at the chronic stage. CD28 signalling in T Helper cells, Phospholipase C signalling and NF-kB signalling were upregulated in the chronic phase in sheep, but downregulated in Qeios, CC-BY 4.0 · Review, July 30, 2021 cattle, while the activation or inhibition of the hepatic fibrosis pathway could not be determined in either species. A total of seven pathways were significant in the chronic phase in sheep and in the acute phase in cattle, and included IL-10 signalling, IL-6 signalling, Th2 pathway and Th1 and Th2 activation pathway, which were all upregulated in sheep. 3.7.2 TLR signalling The TLR pathway was upregulated in sheep in the acute phase and downregulated in the chronic phase of infection, and overrepresented (with no consensus on activation or inhibition) only at the acute stage in cattle. In cattle, the TLR2 and TLR4 genes were upregulated in the chronic phase, and TLR10 was The IPA hepatic fibrosis pathway was enriched at the chronic stages in both sheep and cattle. TGF-β1 is an essential cytokine promoting collagen production and fibrosis, which then leads to encapsulation of flukes and limiting migration of the parasites through the liver parenchyma (Cutroneo, 2007, Behm andSangster, 1999). Liver fibrosis in ruminants has previously been associated with expression of IL10 and TGFB, with an increased expression of these genes potentially leading to increased fibrosis and control of fluke burdens (Haçariz et al., 2009). In cattle, liver fibrosis has been proposed as a downstream effect of upregulated TNF and IL1B, which would activate the genes that were also upregulated in the bovine study, i.e., IL6, PLAU, SERPINE1, TNFRSF1A, SOCS1, and CTSB. In sheep, however, the pathway was observed to exhibit a more recognisable outcome, with the involvement of IFNG, the collagen genes COL11A1 and COL11A2, and matrix metalloprotease MMP9 and SMAD2 genes. Even though typically collagen type I and III are involved in chronic liver fibrosis, the upregulation of genes encoding alpha and beta subunits of collagen type XI as well as a matrix metalloprotease gene, which is involved in reorganisation of the extracellular matrix, supports a potential role of these proteins in fibrosis in chronic liver fluke infection in sheep. Interferons alpha and gamma, typically secreted by NK cells, are also thought to inhibit liver fibrosis (Altamirano- Barrera et al., 2017, Tiggelman et al., 1995. IFNG was downregulated in the chronic phase in sheep, which further supports activation of fibrosis. Surprisingly, SMAD2 was downregulated, which could suggest inhibition of fibrosis, even though other genes in the pathway were activated. IL1B and TLR4 were involved in the pathway in both sheep and cattle, so there could also be an role in liver fibrosis activation for IL1B. TNF [TN28] was not DE in sheep. The SERPINE1 gene, which was detected at chronic stages in the bovine study and in a previous study of ovine PBMC response following infection with F. hepatica (Fu et al., 2016), was also not DE in this study.
The low number of DE genes in the acute response of cattle, coupled with specific anti-inflammatory pathways such as Th2 activation or IL-10 signalling, which were upregulated during the chronic phase in sheep and that were upregulated in the acute phase in cattle, underscore the significance of the acute response in sheep and its potential detrimental effect on disease phenotype. It is important to note that in which then causes a rapid migration of the flukes through the gut and maturation in the liver causing a large inflammatory response, coupled with an attempt to repair tissue damage and slow inflammation through Th2 responses. The conclusion, therefore, would be that a severe acute inflammatory response to F. hepatica, as seen in sheep, has major pathological consequences.
3.7.4 Comparison with a previous ovine PBMC study [TN29] Garcia-Campos et al. (2019) discussed their findings in comparison with a previous transcriptomic study of ovine PBMC by Fu et al. (2016) and concluded that in sheep, the majority of DE genes were present in the acute phase when compared to the chronic phase, while in cattle only 5% of the total DE genes were expressed in the acute phase of infection. While the proportion of acute DE genes in the bovine response is indeed very low, our study shows that using the same longitudinal analysis as in the cattle transcriptomic study, with the inclusion of control animals in the analysis, leads to a result where 63% of all DE genes were allocated to the chronic phase. Even when examining the infected animal group only, as in the study by Fu et al. (2016), the majority of DE genes in our study were expressed during the chronic phase of infection, with 72% of the DE genes detected at this stage. It is important to note, however, that in the previous ovine transcriptomic study by Fu et al. (2016) the DE genes at the chronic phase were compared to the acute phase (i.e., T3 vs T2), not to the baseline of week 0. In the bovine DE gene analysis, the number of DE genes from the chronic phase compared to acute phase (WK14 vs WK1) was three-fold lower than when the chronic phase was compared to baseline. Therefore, the choice of DE gene analysis could explain the skew towards the acute phase in the previous study of the transcriptomic response in sheep.
Breed differences, as well as a differences in the F. hepatica strain used could also partially account for the different results when comparing our study to Fu et al. (2016).
In both our study and the bovine[TN30] study by Garcia-Campos et al. (2019) the log2 fold changes of DE genes were relatively low in comparison to the previous ovine transcriptomic study: in our study the longitudinal DE gene log2 fold changes ranged from −3.1 to 3.3 in the acute phase, and from −3.2 to 4.3 in the chronic phase. While these log2 fold changes are still higher than in the bovine study, the differences are relatively minor, especially given the different DE gene analysis packages used (edgeR in the bovine study vs DEseq2 in the current ovine study), which could lead to slightly different results (Love et al., 2014). The previous ovine study (Fu et al., 2016) observed much larger log2 fold changes, ranging from −27 to 12. The study used a Limma-Voom package for DE gene analysis, which could at least partially account for the differences in log2 fold changes. The DE gene analysis in Limma is based on a linear model, with gene counts typically log-transformed, while DESeq2 and edgeR use a negative binomial generalised linear model, with raw counts used for the analysis. A difference in normalisation of raw counts prior to DE gene analysis could also cause differences in log2 fold changes.

General Discussion
F. hepatica is an important pathogen that causes zoonotic disease and economic losses in animal agriculture. Due to lack of acquired immunity following natural infection (Piedrafita et al., 2004, Piedrafita et al., 2007 and anthelmintic resistance concerns, there is an urgent need for the development of novel control methods such as vaccination. Over the last few decades, progress has been made in the isolation, characterisation and testing of a number of native and recombinant molecules as vaccines against liver fluke disease in ruminant hosts (Toet et al., 2014). However, data from vaccine trials have been inconsistent (Molina-Hernández et al., 2015), and some recent trials with recombinant antigens failed to induce protection against fluke infection[TN31] (Fu et al., 2016) or achieved only a minor reduction of fluke burdens (Zafra et al., 2021). In spite of evidence that effective vaccination may be possible, design of smart vaccines will require more knowledge about the immune response to F. hepatica, which can be provided by big-data approaches such as transcriptomics. This study provides an insight into the response of ovine PBMC to F. hepatica at both the acute and chronic stages of infection, and is the first long-term ovine PBMC transcriptomic study that includes uninfected control animals in the experimental design.

Comparison with previous ovine transcriptomics studies[TN32]
Our study confirmed the need to use control animals in the infection. A greater number of DE genes were identified in sheep in response to F. hepatica by Fu et al. (2016). However, this previous study did not use uninfected control animals in its design, and DE genes were compared to pre-infection. In fact, when the DE gene analysis was performed in this study in comparison to time 0 (longitudinal analysis), a large number of DE genes was identified in both infected and control animals, a phenomenon also found in a previous bovine study ( Garcia-Campos et al., 2019). This suggests that presence of DE genes related to developmental changes that occur in 9 month sheep could obscure the response to infection. Therefore, a study on older animals may be warranted, where growth effects may not be so dominant. Another possible explanation may be the seasonality of the infection and the onset of low temperatures. A previous study by Liu et al. (2014) identified DE genes in sheep due to wool production between the months of August and December. The resulting DE genes included IL1B, IL5, IL8 and CXCL1, which suggests that immune proteins may be produced or inhibited due to seasonal changes. Our study was conducted between September and January, which could suggest a similar seasonal effect. In addition, immune response in sheep could be breed-dependent, as demonstrated by relative resistance of the Indonesian thin tail sheep breed against F. gigantica when compared to Merino sheep (Hansen et al., 1999, Pleasance et al., 2011.
Our study was conducted on Merino sheep, while the previous ovine study used the Suffolk breed.
Longitudinal studies of gene expression in sheep with seasonal changes, potentially in conjunction with measuring hormones and other physiological variables could assist future transcriptomic studies of host response.
In a previous ovine PBMC study by Alvarez Rojas et al. (2016), genes encoding for galectin 9, SMAD2 and TGF-β were all upregulated at 2 wpi. The genes and pathways described in Alvarez Rojas et al. (2016) differed from genes and pathways found in this study at 2 wpi.  Rojas et al. (2016) was performed using NOISeq, which has been shown to produce a comparably large number of false positive results when less than six biological replicates are used (Schurch et al., 2016).
Another ovine transcriptomic study was also performed on the response to F. hepatica in the HLN at 16 wpi (Naranjo-Lucena et al., 2021), in the same sheep as the ones used in our study. There were common pathways found in PBMC in the infected vs control analysis in this study at 2 wpi and in the HLN at 16 wpi.
Specifically, B cell receptor signalling was upregulated in both tissues, and apoptosis and NK cell signalling were downregulated. Crosstalk between DCs and NK cells was downregulated in the hepatic lymph node, while in this study the direction of pathway activation could not be determined, mainly due to a low number of DE genes involved (Supplementary Table 4). In our study, the CD-NK pathway was mainly significant due to the activation of DC receptor genes CD69, CD80 and CD83. It is important to note that the comparison is made between the acute phase in PBMC and chronic phase in HLN, and therefore the DC-NK cross-talk pathway findings in both studies could be an indication of disease progression from activation of DCs in the acute phase, potentially caused by recognition of parasitic RNAs by TLR9, followed by an inhibition of DC activity in the chronic phase. The involvement of DCs and NK cells in response to F.
hepatica is a novel finding in both our study and the study on response to F. hepatica in HLN. DCs and NK cells are present in the liver as a part of the resident immune cell population (Robinson et al., 2016).
Therefore, the involvement of these cell types in the response of PBMC and lymph node lymphocytes to F.
hepatica could indicate that these are important cell types responding to the parasite in the liver, and that these cells could be targeted for stimulation of an adaptive immune response.
The aforementioned ovine HLN study used an infected vs control analysis of DE genes, and over 5000 DE genes were detected at 16 wpi. However, in our PBMC study only 2 DE genes were found in the infected vs control analysis at 16 wpi. Therefore, a comparison of the chronic response in PBMC and HLN in the present study and the study by Naranjo-Lucena et al. (2021) was not possible. In the HLN study, fibrosis was not found in the lymph node tissue of the infected sheep, and the genes TGFB1 and COL1A1 typically associated with liver fibrosis found in the previous ovine studies (Fu et al., 2016, Alvarez Rojas et al., 2015 were not DE in the HLN. In our PBMC study, these genes were also not DE, even though their expression would be expected to increase due to the "echoing" capacity of PBMC. Fibrosis of the liver was not extensive in the sheep from our study, which is confirmed by the lack of significant overexpression of fibrosis-related genes.

Limitations of longitudinal analysis
Larger numbers of DE genes were found in the longitudinal analysis than when the infected and control groups were compared at each time point. However, it is important to note that the P values and foldchanges of these longitudinal DE genes relate to a comparison between the infected animals at week 2 or 16 versus week 0. The statistical difference between the infected and control group is not accounted for in Qeios,.0 · Review, July 30, 2021 Qeios ID: 3NP14M · https://doi.org/10.32388/3NP14M 19/27 this case and rather the change in response over time in each group is assessed. This analysis was performed to highlight the potential of developmental or seasonal changes in both groups of sheep, which may lead to a small number of DE genes when infected and control animals are compared. It also highlights the necessity of a suitable control in transcriptomics experiments, where changes in gene expression may be caused by factors other than the infection studied.
The infected vs control analysis identified genes where expression difference were high enough enough to overcome individual animal differences between the infected and control animal groups. The low number of DE genes found in the main infected vs control analysis also indicates that the systemic response to F. hepatica was muted, since larger numbers of DE genes were found in the lymph nodes of the same animals at 16 wpi (Naranjo-Lucena et al., 2021).

Concluding remarks
This is the first study which describes the transcriptomic response of ovine PBMC to F. hepatica with infection duration up to 16 wpi and with control animals included in the experimental design. The study was systematically compared to a similar study in cattle, both of which represent major livestock species that are the targets for development of vaccines as a control measure for fasciolosis. Overall, the infection vs control analysis in our study revealed downregulation of apoptosis and increased cellular metabolism in the acute phase, with the upregulation of the TP53 pathway, and a response to intracellular foreign nucleic acids through TLR9 and RIG-I in the acute phase of infection. The intracellular immune response to RNA molecules could indicate involvement of F. hepatica extracellular vesicles in triggering the acute phase response. Based on the most significant DE genes including CD83, expression of TLR9 and a significant DC-NK cross-talk pathway, the involvement of DCs may be key to the acute phase response in sheep. A further longitudinal analysis revealed that an anti-inflammatory response tends to occur in sheep only in the chronic phase of infection, while in cattle the same response is activated as early as 1 wpi. This reflects the different phenotype of disease in the two host species, and may point to new control strategies for sheep, such as a potential benefit of using adjuvants in ovine vaccines to trigger an early anti-inflammatory response, or triggering dendritic cell activation. Overall, this study leads to better understanding of hostparasite [TN33] relationships.
Breed differences and the infecting strain may have affected the DE genes found in this study when compared to other studies. Furthermore, it is apparent that differences in experimental design and bioinformatic approaches used, including selection of controls, use of software packages, or inclusion of batch effects can significantly affect the resulting DE genes. It is recommended that uninfected control animals are used in the experimental design to account for potential seasonal or developmental differences, especially in long term studies. [TN34] 5 Data Availability Statement The datasets generated for this study can be obtained from the European Nucleotide Archive (www.ebi.ac.uk/ena) with accession PRJEB45790. Qeios, CC-BY 4.0 · Review, July 30, 2021 This experiment was approved by the Bioethics Committee of the University of Córdoba (UCO, Spain) (code No. 1118) and conducted in accordance with European (2010/63/UE) and Spanish (RD 1201(RD /2005 Directives on animal experimentation. Written informed consent was obtained from the owners for the participation of their animals in this study.

Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. 10 Acknowledgments We would like to thank the staff of the experimental facility from the Universidad de Córdoba for the care of the animals and help provided during the trial.