Transcriptome Analysis for Genes Associated with Small Ruminant Lentiviruses Infection in Goats of Carpathian Breed

Small ruminant lentiviruses (SRLV) are economically important viral pathogens of sheep and goats. SRLV infection may interfere in the innate and adaptive immunity of the host, and genes associated with resistance or susceptibility to infection with SRLV have not been fully recognized. The presence of animals with relatively high and low proviral load suggests that some host factors are involved in the control of virus replication. To better understand the role of the genes involved in the host response to SRLV infection, RNA sequencing (RNA-seq) method was used to compare whole gene expression profiles in goats carrying both a high (HPL) and low (LPL) proviral load of SRLV and uninfected animals. Data enabled the identification of 1130 significant differentially expressed genes (DEGs) between control and LPL groups: 411 between control and HPL groups and 1434 DEGs between HPL and LPL groups. DEGs detected between the control group and groups with a proviral load were found to be significantly enriched in several gene ontology (GO) terms, including an integral component of membrane, extracellular region, response to growth factor, inflammatory and innate immune response, transmembrane signaling receptor activity, myeloid differentiation primary response gene 88 (MyD88)-dependent toll-like receptor signaling pathway as well as regulation of cytokine secretion. Our results also demonstrated significant deregulation of selected pathways in response to viral infection. The presence of SRLV proviral load in blood resulted in the modification of gene expression belonging to the toll-like receptor signaling pathway, the tumor necrosis factor (TNF) signaling pathway, the cytokine-cytokine receptor interaction, the phagosome, the Ras signaling pathway, the phosphatidylinositol 3-kinase (PI3K)/protein kinase B (AKT) (PI3K-Akt) signaling pathway and rheumatoid arthritis. It is worth mentioning that the most predominant in all pathways were genes represented by toll-like receptors, tubulins, growth factors as well as interferon gamma receptors. DEGs detected between LPL and HPL groups were found to have significantly enriched regulation of signaling receptor activity, the response to toxic substances, nicotinamide adenine dinucleotide (NADH) dehydrogenase complex assembly, cytokine production, vesicle, and vacuole organization. In turn, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway tool classified DEGs that enrich molecular processes such as B and T-cell receptor signaling pathways, natural killer cell-mediated cytotoxicity, Fc gamma R-mediated phagocytosis, toll-like receptor signaling pathways, TNF, mammalian target of rapamycin (mTOR) signaling and forkhead box O (Foxo) signaling pathways, etc. Our data indicate that changes in SRLV proviral load induced altered expression of genes related to different biological processes such as immune response, inflammation, cell locomotion, and cytokine production. These findings provide significant insights into defense mechanisms against SRLV infection. Furthermore, these data can be useful to develop strategies against SRLV infection by selection of animals with reduced SRLV proviral concentration that may lead to a reduction in the spread of the virus.


Introduction
The maedi visna virus (MVV) and the caprine arthritis encephalitis virus (CAEV) belong to the group called small ruminant lentiviruses (SRLV) within the Retroviridae family. Molecular epidemiology studies revealed that both viruses represent a broad spectrum of genetic variants that can infect sheep and goats. To date, four main genotypes have been described, but molecular information on new subtypes is continuously updated, showing a high genetic and antigenic heterogeneity [1]. Infections with SRLV, which are spread worldwide, cause multi-organ failure usually over a long period of time and can lead to severe diseases such as pneumonia, mastitis, arthritis, wasting, and encephalitis [2]. Moreover, they contribute to economic losses in small ruminant production and affect animal welfare deterioration [3,4].
There is no effective vaccine or treatment preventing animals from SRLV infection. Several practices for controlling or preventing SRLV infection have been developed, such as serological testing with culling or segregation of infected animals, replacement of infected animals with offspring from seronegative mothers, or artificial rearing of newborn animals separated from the infected mothers immediately after birth [5,6]. These practices can be effective when carefully designed and applied continuously to eradicate the progression of the infection [6][7][8][9][10]. However, such an approach is often costly and time-consuming. The high genetic variability of SRLV and the absence of sensitive diagnostic tests that are able to detect all strains are additional challenges reducing the effective implementation of eradication programs.
In many persistent viral infections, viral load is reported to estimate the likelihood of pathogenesis and disease progression. For retroviruses, including lentiviruses, in which genomes are integrated with the host genome, proviral load (PL) is a risk factor determining disease prediction [34]. It was shown that animals with high PL showed more tissue lesion severity, indicating that proviral concentration is positively correlated with the presence and severity of clinical disease symptoms [35,36]. Elimination of animals predisposed to high PL can limit the outcome of clinical signs and spread of the virus since these animals are also highly efficient in shedding the virus [37]. On the contrary, some studies indicated potential restriction in low PL carriers and referred to them as long-term non-progressors. These animals showed competent antibody response in the absence of productive virus replication leading to minimizing the spread of the virus within the flock [38]. Consistent with this, approaching genetics factors associated with low PL in animals infected with SRLV could be used to control SRLV infection, especially in flocks with a high level of seroprevalence.
In this study, RNA sequencing (RNA-seq) was used to identify genes associated with high (HPL) and low (LPL) proviral load in goats of Carpathian breed naturally infected with SRLV. The results provide unique insights for further exploration and understanding patterns of the host responses to SRLV infection in goats. A deep understanding of Viruses 2021, 13, 2054 3 of 20 transcriptome profile changes may provide additional information on the contribution of specific genes responsible for the course of infection with SRLV. Additionally, these data can be useful to develop strategies against SRLV infection by elimination and/or selection of animals with reduced SRLV provirus concentration, which may lead to limitation of viral spread and can improve the welfare of animals and prolong their life.

Animals and Blood Sample Collection
Whole blood was taken from 27 adult goats by jugular venipuncture and stabilized in ethylenediaminetetraacetic acid (EDTA) and Tempus blood RNA tubes. Goats represented the Carpathian breed, and they were owned by the National Research Institute of Animal Production in Krakow. All goats were healthy and maintained in one flock in the same environment. This involved being housed indoors in sheds, aside from during the grazing season (April-November), where they spent daily hours in pastures outdoors. It also involved the feeding conditions (summer feeding based on pasture or green fodder and winter feeding consisting mainly of hay and oats). Serological status of animals for SRLV infection was confirmed by enzyme-linked immunosorbent assay (ELISA) (ID Screen MVV/CAEV Indirect Screening test, IDVet, Grabels, France) according to the manufacturer's recommendations. Blood samples were collected from all animals on the same day. At the time when blood was taken, none of the goats exhibited any clinical signs of the disease. All procedures associated with animal handling and treatments were approved (no 37/2016) by the Local Ethical Committee on Animal Testing at the University of Life Sciences in Lublin (Poland).

Proviral Load Quantification
DNA was extracted from peripheral blood leukocytes (PBLs), and proviral DNA was quantified by the real-time polymerase chain reaction (PCR) using Rotor-Gene Q Series ver. 2.0.3 (Qiagen, Hilden, Germany) with primers and probe specifically designed for SRLV A5 subtype, which circulation in this flock was previously confirmed [39]. Sequences of forward and reverse primers and probe were CA5F (5 TGGGAGTAGGA-CAAACAAATCA 3 ), CA5R (5 TGACATAT GCCTTACTGCTCTC 3 ) and CA5P (5 6-FAM-TCACCCATTGTAGGCATAGCTGCC-BHQ-1 3 ), respectively. A reference plasmid encompassing the target gag region was generated by the cloning of a 625 bp fragment into pDrive plasmid used to generate a standard curve based on 10-fold serial dilutions of plasmid DNA from 10 8 to 1 0 . Amplification was performed in a total volume of 20 µL, according to the following cycling conditions: initial incubation and polymerase activation at 95 • C for 15 min and followed by 45 cycles of 94 • C for 60 s and 60 • C for 60 s. The reaction mixture for each PCR test contained 10 µL 2× QuantiTect Multiplex NoROX PCR buffer (Qiagen, Hilden, Germany), 400 nM of each of the primers, 200 nM of the specific probe, 5 µL of extracted genomic DNA. A non-template control (diethylpyrocarbonate (DEPC) H 2 O) was included in each run. All samples were tested in duplicate, and the results were expressed as a mean copy number of provirus per 500 ng of genomic DNA of each goat. Then, these data were used for further statistical analysis, thereby allowing the identification of goats with a high (HPL) and low (LPL) proviral load.

Transcriptome Sequencing and Data Analysis
The total RNA was isolated from the whole blood of goats using MagMAX™ for Stabilized Blood Tubes RNA Isolation Kit (Thermo Fisher Scientific, Waltham, MA, USA ), according to the protocol. The possible RNA contamination with DNA was removed using TURBO DNase™ (Thermo Fisher Scientific, Waltham, MA, USA). The quality and quantity of obtained RNA were checked using the Nanodrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA) and TapeStation 2200 System (Agilent, Santa Clara, CA, USA) using Agilent RNA ScreenTape (Agilent, Santa Clara, CA, USA). The samples with an RIN (RNA integrity number) value above 8 were used for further analysis. The cDNA Viruses 2021, 13, 2054 4 of 20 libraries were prepared from 300 ng of total RNA with the use of the TruSeq RNA Kit v2 kit protocol (Illumina, San Diego, CA, USA). Each library was ligated with adaptors with different index barcodes. The quality and quantity of libraries were assessed using Qubit 2.0 fluorometer (Invitrogen, Carlsbad, CA, USA ) and TapeStation 2200 (Agilent, Santa Clara, CA, USA) with D1000 ScreenTape (Agilent, Santa Clara, CA, USA). Libraries were pooled and sequenced by synthesis, using HiSeq High-Output v4-SR (Illumina, San Diego, CA, USA) into 50 single-end cycles, according to the protocol. The quality of the reads was assessed with FastQC software [40]. Then, we used Flexbar software [41] to remove adapters, reads shorter than 35 base pairs, and those with a phred quality score lower than 30. Processed reads were mapped to the goat reference genome Capra hircus ARS1 (GCA_001704415.1) with Tophat software [42] on default parameters. Next, the mapped reads were counted into Ensembl GTF version 97 annotation intervals using HTSeq-count software [43]. Differential expressed genes (DEG) were estimated using DESeq2 software [44] with default parameters. Genes with p-adjusted < 0.05 (Benjamini-Hochberg p-value adjustment) and fold change >1.3 were regarded as differentially expressed and included in further annotation analysis.

GO Enrichment and Pathways Analysis
The gene ontology analysis (GO) and pathways analyses were performed on all significant differentially expressed genes (DEGs) sets. The gene set enrichment analysis (GSEA) was performed with the use of WebGestalt software with Fisher's exact test (http://webgestalt.org/ accessed on 31 July 2021). For pathway functional analysis, David software v6.8 (Fisher Exact test) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) with KEGG mapper search pathway tools were applied. The significantly enriched pathways were identified based on the p-values obtained from a Fisher exact test [45]. The latest version of the Capra hircus ARS1 (GCF_001704415.1) reference genome was used.

Quantitative Polymerase Chain Reaction (qPCR) Analysis
Validation of the RNA-seq results was carried out using the real-time PCR method for nine DEGs (Supplementary Materials Table S1), selected based on their important function in viral infection and/or previous reports confirming their association with proviral concentration. The qPCR was performed for the validation of RNA-seq results of all samples analyzed using the RNA-seq method (for validation), as well as for all samples tested by qPCR, which was tasked with the estimation of the correlation between PL (SRLV copy number) and gene expression levels (for a correlation with SRLV copy number). The cDNA was prepared from 250 ng of total RNA using a high-capacity RNA-to-cDNA Kit (Thermo Fisher Scientific, Waltham, MA, USA) according to protocol. The transcript level of selected genes was estimated on QuantStudio 7 Flex (Applied Biosystems, Thermo Fisher Scientific, Waltham, MA, USA), and for each gene, the reaction was carried out in three replications using Sensitive RT HS-PCR Mix EvaGreen (A&A Biotechnology, Gdynia, Poland). The expression was calculated using the delta-delta CT method according to Pfaff [46] and based on HPRT1 and ACTB reference controls [47]. The comparison between next-generation sequencing (NGS) data (RNA-seq) and relative quantity obtained by real-time PCR method was performed using Spearman correlation with the use of R software [48].

Statistical Analysis
To classify the animals into HPL and LPL groups, a copy number calculated by qPCR per each animal was used as a potential cut-off value, and a Box-Cox transformation was used to achieve normal distribution. The t-Student test was employed to determine the significance of the differences between the two potential groups of animals. Finally, the cut-off value was chosen based on the lowest p-value to distinguish between the HPL and LPL groups, which was additionally confirmed by the Welch test (p < 0.0001).
The phenotypic and physiological variables between HPL and LPL goats were analyzed using TIBCO Software Inc. Statistica (Data Analysis Software System, Palo Alto,  13 (2017). The association between physiological differences (stage of lactation, mastitis and other diseases, abortion) between HPL and LPL goats were tested by a chi-square test with appropriate correction. The variables included age, parity number, body weight, and milk yield (kg) and were analyzed with Pearson correlation. All variables were also tested by logistic regression models. The normality of the distribution and the homogeneity of variance was tested by Shapiro-Wilk and Brown-Forsyth tests, respectively. Differences were considered significant when p > 0.05.

Classification of Goats on HPL and LPL
A total of 24 animals from the flock tested in this study were seropositive by ELISA and positive to quantitative polymerase chain reaction (qPCR), whereby confirming the infection of SRLV. Three goats were negative in both ELISA and qPCR. The average number of proviral copies varied from 1 to 106 per 500 ng of genomic DNA, and these values showed skewed distribution with a relatively limited number of animals with a high concentration of provirus. Animals were classified into HPL (mean ± SD; 82.39 ± 13.14) and LPL group (mean ± SD; 14.31 ± 14.23) ( Figure 1). significance of the differences between the two potential groups of animals. Finally, th cut-off value was chosen based on the lowest p-value to distinguish between the HPL an LPL groups, which was additionally confirmed by the Welch test (p < 0.0001).
The phenotypic and physiological variables between HPL and LPL goats were ana lyzed using TIBCO Software Inc. Statistica (Data Analysis Software System, Palo Alto, CA USA), version 13 (2017). The association between physiological differences (stage of lacta tion, mastitis and other diseases, abortion) between HPL and LPL goats were tested by chi-square test with appropriate correction. The variables included age, parity numbe body weight, and milk yield (kg) and were analyzed with Pearson correlation. All varia bles were also tested by logistic regression models. The normality of the distribution an the homogeneity of variance was tested by Shapiro-Wilk and Brown-Forsyth tests, re spectively. Differences were considered significant when p > 0.05.

Classification of Goats on HPL and LPL
A total of 24 animals from the flock tested in this study were seropositive by ELISA and positive to quantitative polymerase chain reaction (qPCR), whereby confirming th infection of SRLV. Three goats were negative in both ELISA and qPCR. The average num ber of proviral copies varied from 1 to 106 per 500 ng of genomic DNA, and these value showed skewed distribution with a relatively limited number of animals with a high con centration of provirus. Animals were classified into HPL (mean ± SD; 82.39 ± 13.14) an LPL group (mean ± SD; 14.31 ± 14.23) ( Figure 1). In addition, the statistical analyses of the phenotypic and physiological difference between HPL and LPL goats were performed. Variables analyzed included age, parit number, body weight, milk yield, stage of lactation, disease occurrence (including mast tis), and abortion. A moderate correlation between proviral load and age (r = 0.35, p 0.035), as well as parity number (r = 0.38, p = 0.024), was observed; however, no significan In addition, the statistical analyses of the phenotypic and physiological differences between HPL and LPL goats were performed. Variables analyzed included age, parity number, body weight, milk yield, stage of lactation, disease occurrence (including mastitis), and abortion. A moderate correlation between proviral load and age (r = 0.35, p = 0.035), as well as parity number (r = 0.38, p = 0.024), was observed; however, no significant differences between HPL and LPL goats were noted when other variables were analyzed. Finally, eight goats (four goats with HPL and four goats with LPL) that were phenotypically and physiologically the most homogeneous were then carefully selected for whole blood transcriptome analysis. All these animals were female, unrelated within biological groups, clinically healthy, and multiparous after parturition. Their average age was 7.25 ± 1.98 years, the average body weight was 45.3 ± 5.15 kg, and they produced on average 266.6 ± 96.0 kg of milk per year.
Moreover, partial (gag and LTR) sequences of the SRLV genome were analyzed, and no significant mutations/differences between sequences obtained from HPL and LPL goats that could alter proviral load were observed.

Transcriptome Quantification
Transcriptome analysis was performed using whole blood taken from 11 selected goats. These animals included four goats with high (HPL group) and four goats with low (LPL group) SRLV proviral load, as well as three uninfected goats (control group). After next-generation sequencing (NGS) and data filtration, the average number of reads obtained per sample was about 34.3 mln. On average, 83.8% of reads were mapped to the reference Capra hircus genome (GCA_001704415.1) (Supplementary Materials Table S2). Furthermore, the principal component analysis (PCA) showed that the analyzed groups of animals formed distinct clusters, which confirmed the presence of two groups of goats with high and low proviral concentrations (Supplementary Materials Figure S1).

DEGs Analysis
The whole blood transcriptome sequencing using the NGS approach allowed the identification of 1130 DEGs between control and LPL groups, 411 between control and HPL groups, and 1434 significant DEGs between HPL and LPL groups ( Figure 2). no significant mutations/differences between sequences obtained from HPL and goats that could alter proviral load were observed.

Transcriptome Quantification
Transcriptome analysis was performed using whole blood taken from 11 selec goats. These animals included four goats with high (HPL group) and four goats with (LPL group) SRLV proviral load, as well as three uninfected goats (control group). A next-generation sequencing (NGS) and data filtration, the average number of reads tained per sample was about 34.3 mln. On average, 83.8% of reads were mapped to reference Capra hircus genome (GCA_001704415.1) (Supplementary Materials Table  Furthermore, the principal component analysis (PCA) showed that the analyzed gro of animals formed distinct clusters, which confirmed the presence of two groups of g with high and low proviral concentrations (Supplementary Materials Figure S1).

DEGs Analysis
The whole blood transcriptome sequencing using the NGS approach allowed identification of 1130 DEGs between control and LPL groups, 411 between control HPL groups, and 1434 significant DEGs between HPL and LPL groups ( Figure 2). When the set of DEGs between uninfected and both HPL and LPL goats was c pared, the 1046 genes were detected. Among this DEGs panel, 408 genes were identi as downregulated in both HPL and LPL groups, while 638 genes were upregulated. T When the set of DEGs between uninfected and both HPL and LPL goats was compared, the 1046 genes were detected. Among this DEGs panel, 408 genes were identified as downregulated in both HPL and LPL groups, while 638 genes were upregulated. This gene set was used for subsequent analysis as being potentially associated with immune response to SRLV infection. Among all 1434 DEGs identified between LPL and HPL goats, 571 were upregulated and 863 downregulated.
When the groups of LPL vs. HPL were compared, the highest differences in gene expressions were detected for downregulated genes for which the decrease in transcriptlevel abundance reached up to 419-fold change (FC) for solute carrier family 22 member 1 (SLC22A1) gene. The most significant deregulated genes showing the highest FC in gene expression between analyzed groups are presented in Figure 3 and Supplementary Materials Table S3. The group of genes for which expression was deregulated in response to the provirus concentration were: the zinc-finger gene family (ZNF; 23 genes); a transmembrane protein (TMEM; 16 genes); the solute carrier family genes (SLC; 22 genes); the NADH: ubiquinone oxidoreductase supernumerary subunits (NDUF) genes (20 genes); ATP synthase genes (15 genes); interleukin and interleukin receptors (12 genes); the sorting nexin family (SNX; 5 genes) and the translocase of inner mitochondrial membrane family (TIMM, 5 genes).
When the groups of LPL vs. HPL were compared, the highest differences in gene expressions were detected for downregulated genes for which the decrease in transcriptlevel abundance reached up to 419-fold change (FC) for solute carrier family 22 member 1 (SLC22A1) gene. The most significant deregulated genes showing the highest FC in gene expression between analyzed groups are presented in Figure 3 and Supplementary Materials Table S3. The group of genes for which expression was deregulated in response to the provirus concentration were: the zinc-finger gene family (ZNF; 23 genes); a transmembrane protein (TMEM; 16 genes); the solute carrier family genes (SLC; 22 genes); the NADH: ubiquinone oxidoreductase supernumerary subunits (NDUF) genes (20 genes); ATP synthase genes (15 genes); interleukin and interleukin receptors (12 genes); the sorting nexin family (SNX; 5 genes) and the translocase of inner mitochondrial membrane family (TIMM, 5 genes).

DEGs Detected between Control Group and Groups with Proviral Load
The gene ontology analysis showed the significant enrichment of several GO terms ( Table 1). The most overrepresented were integral components of the membrane, which involved 13 DEGs, 6 extracellular region genes, and 5 response to growth factor genes. Other detected GO terms related to inflammation and innate immunity were represented by four genes belonging to the toll-like receptors family, TLR2, TLR4, TLR8, and TLR7, which were all significantly upregulated in the infected goats (Figure 4).
Our results also showed significant deregulation of selected pathways in response to viral infection ( Table 2). The presence of SRLV proviral load in blood cells resulted in modification of expression in genes belonging to toll-like receptor signaling pathway (FDR > 0.05); TNF signaling pathway (FDR > 0.01); Cytokine-cytokine receptor interaction (FDR > 0.04) and phagosome (FDR > 0.02) ( Figure 5). It is worth mentioning that the most predominant genes in all pathways were the genes represented by toll-like receptors, tubulins, growth factors, as well as interferon gamma receptors. The highest number of downregulated genes were detected within the Ras signaling pathway. These pathways allowed the identification of PLA2G1B (phospholipase A2 group IB) and KITLG (KIT ligand) DEGs, and both were considered as strongly related to viral infection.  EER REVIEW 9 of 20 Figure 4. The interaction between differentially expressed genes involved in toll-like receptor signaling and cytokine-cytokine receptor interaction pathways (red-genes belonged to NF-kappa-B signaling pathway, and TIR domain; dark blue-I-kappa-B kinase/NF-kappa-B signaling, and interleukin-1 receptor binding; yellow-TNFR1-induced NF-kappa-B signaling pathway, and TICAM1 deficiency-HSE; purple-innate immunity; green-cytokine; light blue-inflammatory responses; String software; detected genes showed no more than five interactions).
Our results also showed significant deregulation of selected pathways in response to  False discovery rate (p-value adjusted).

Figure 5.
DEGs for which expression has been modified through SRLV infection involved in phagosome pathways (KEGG chx04145). The genes identified as differentially expressed (adjusted p-value < 0.05) between uninfected and infected goats were highlighted red.

DEGs Detected Between LPL and HPL Groups
GO enrichment analysis allowed the detection of the most represented GO term, 10 up-and 10 downregulated ( Figure 6).

Gene
Protein Name FC Adjpval Protein Function

TLR2
Toll-like receptor 2 1.60 0.03 Related to mediating the innate immune response to bacterial lipoproteins or lipopeptides, related to cytokine secretion and the inflammatory response.

TLR6
Toll-like receptor 6 1.66 0.04 Acts via MYD88 and TRAF6, leading to NF-kappa-B activation, cytokine secretion, and the inflammatory response.

CHUK
Inhibitor of nuclear factor kappa-B kinase subunit alpha 2.44 0.03 Plays an essential role in the NF-kappa-B signaling pathway activated by multiple stimuli also by viral products.

CSF1R
Macrophage colony-stimulating factor 1 receptor 1.80 0.03 Controlling the proliferation and differentiation of hematopoietic precursor cells, especially mononuclear phagocytes, such as macrophages and monocytes.

IRF1
Interferon regulatory factor 1 1.41 0.02 Regulation of IFN and IFN-inducible genes, host response to viral and bacterial infections.

NRLP3
NACHT, LRR, and PYD domain-containing protein 3 Plays a crucial role in innate immunity and inflammation.

IFIH1
Interferon-induced helicase C domain-containing protein 1 2.14 0.05 Plays a major role in sensing viral infection and in the activation of a cascade of antiviral responses, including the induction of type I interferons and proinflammatory cytokines.

TBK1
Serine/threonine-protein kinase TBK1 1.80 0.02 Regulation of transcriptional activation of proinflammatory and antiviral genes including IFNA and IFNB.

MYD88
Myeloid differentiation primary response protein MyD88 1.82 0.01 Acts via toll-like receptor and IL-1 receptor signaling pathway in the innate immune response.
Identified DEGs were also analyzed for their involvement in biological pathways. Thus, genes have been assigned to pathways involved in the acquired or antigen-specific immune response (B-cell receptor and T-cell receptor signaling pathways; natural killer cell-mediated cytotoxicity and Fc gamma R-mediated phagocytosis). Moreover, DEGs belonged to the pathways responsible for recognition of pathogen, signal transduction, and early immune responses: toll-like receptor signaling pathway; tumor necrosis factor (TNF) signaling pathway; mammalian target of rapamycin (mTOR) signaling; and forkhead box O (Foxo) signaling pathway. The comparison of the whole blood transcriptome of goats with different provirus copy numbers allowed the detection of the Ras signaling pathway (17 DEGs), inflammatory mediator regulation of transient receptor potential (TRP) channels (12 DEGs), and hypoxia-inducible factor 1 (HIF-1) signaling pathway (12 DEGs), which are considered as critical to control cytokine production, cell differentiation, function, and cytotoxicity. The panel of genes was identified (paxillin (PXN); profilin 1(PFN1); actinrelated protein 2/3 complex subunit 5 (ARPC5); cytoplasmic FMR1 interacting protein 1 (CYFIP1); IQ motif containing GTPase activating protein 1 (IQGAP1)), which also represent regulation of the actin cytoskeleton. The genes belonging to the selected pathways are presented in Table 6.

qPCR Results
DEGs from different functional groups, including the following genes: CCL2, CXCL5, IL15, C-X-C motif chemokine receptor 3 (CXCR3), MIF, NOD2, CCR, B-cell lymphoma 2 (BCL2), and IL-2-inducible T-cell kinase (ITK), were selected for further validation by qRT-PCR. This analysis revealed an agreement with the RNA-seq results: a high and significant correlation was detected for IL15, CXCR3, and NOD2 genes ( Table 7). For other genes, the correlation was not significant, which may be related to genome annotation still being under development and continued limited knowledge of all spliced variants of the studied genes.
The correlation between provirus copy number and gene expression levels carried out using samples from all animals tested from the flock showed that selected DEGs as CCL2 and CXCL5 (p-value < 0.001), CCR and BCL2 (p-value < 0.01), and ITK and NOD2 genes (p-value < 0.05) were significantly associated with SLRV copy number.

Discussion
To better understand the role of genes involved in the host response to SRLV infection, the RNA-seq method was applied to compare the whole gene expression profile in uninfected goats with those carrying relatively high and low SRLV proviral loads. Data obtained in this study enabled us to identify 1130 DEGs between control and LPL groups, 411 between control and HPL groups, and 1434 significant genes showing changed expression levels depending on provirus copy number. Out of the panel of 1434 DEGs differentiated HPL and LPL goats, only 10 DEGs were shared between both the control vs. LPL and control vs. HPL groups. This indicates that proviral load might be the main driver and risk factor determining disease prediction in goats infected with SRLV [35]. Here, we focused on the analysis of some DEGs only being involved in immunological processes since both innate and adaptive immune responses are known to play a crucial role in controlling the course of SRLV infection.
It was shown that SRLV infection influences the expression of a cytokine network that plays a pivotal role in the activation of the immune system and SRLV-related pathogenesis [20]. Our findings indicated that 95 of the DEGs that were involved in multiple biological processes of cytokine production were overexpressed in HPL goats. These animals showed upregulated expression of interleukin 15 (IL-15) and interleukin 1 alpha (IL-1α) and receptors for IL-10 (IL10Rβ), IL-13 (IL13Rα1), IL-15 (IL15Rα), IL-2 (IL2Rα) and IL-4 (IL4R). This observation partly confirmed results obtained by Ravazzolo et al. [35], who did not find prominent differences in the expression of several interleukins in goats with different SRLV proviral loads. The level of IL-15 seems to be associated with the proviral load, as was also seen in patients infected with human immunodeficiency virus (HIV) with high viral load [49]. There is limited knowledge about the expression of IL-1α in the course of SRLV infection. Jarczak et al. [50] observed down-regulation of IL-1α mRNA in the blood of infected goats, suggesting that lentivirus infection may inhibit the expression of this gene. However, this fact was not confirmed in our study where IL-1α was upregulated in HPL goats. IL-1α is a proinflammatory cytokine that induces the expression of a variety of genes and synthesis of several proteins, which, in turn, induce acute and chronic inflammatory changes [51,52]. However, animals tested in this study did not show any clinical signs of infections, but we cannot exclude the presence of inflammatory processes, especially in goats with HPL, as the association between virus load and presence of inflammatory lesions was clearly evidenced [35,53].
Among the most interesting DEGs detected in this study and engaged in numerous biological processes of cytokine production were also toll-like receptor 2 (TLR2), toll-like receptor 4 (TLR4), toll-like receptor 6 (TLR6), cluster of differentiation 14 (CD14), and myeloid differentiation primary response gene 88 (MyD88), which are involved in TLR signaling. Toll-like receptors, a family of pattern recognition receptors (PRRs), are key elements of native immunity. While the role of SRLV-induced TLR signaling has not been widely studied in sheep and goats, it was shown that mutations in TLR7 and TLR8 may play an important role in susceptibility and/or resistance to SRLV infection [54,55]. Our results indicated that four genes belonging to the toll-like receptors family, TLR2, TLR4, TLR8, TLR7, were significantly upregulated in the infected goats compared to uninfected goats. However, we did not observe different expressions of TLR7 and TLR8 genes between goats with HPL and LPL. Only the genes encoding TLR2, TLR4, and TLR6 were differentially expressed and were found to be upregulated in HPL groups. TLR2, TLR4, and heterodimers TLR2-TLR6 are TLR family members that have been involved in the recognition of viral structural and nonstructural proteins leading to inflammatory cytokine production [56,57]. The expression of CD14 (a co-receptor for the TLR4 and TLR2 response [58] and MyD88), an adaptor molecule that is critical for the signaling responses initiated through most TLRs, was also upregulated in HPL goats. We can conclude that immune response against SRLV is at least partially dependent upon TLR2 and TLR4 and correlated with the concentration of proviral DNA, as was shown in the HIV model based on the expression of TLR2 and TLR4 in monocytes [59].
Interferons (IFNs α, β, and γ) response is a highly robust and effective first line of defense against a wide variety of viral infections; however, results on IFNs expression during SRLV infection are contradictory [20,35,60,61]. When the IFN is synthesized, it binds to the interferon alpha receptor (IFNAR), the specific receptor for IFN-I on the cell membrane, formed by two subunits: IFNAR-1 and IFNAR-2. This binding activates the tyrosine kinases TYK-2 and JAK-1, leading to the activation of the JAK-STAT pathway, which is important in cytokine-mediated immune responses [62]. In our study, no differences in IFNs expression were observed between infected and uninfected goats, as well as between HPL and LPL goats. However, genes involved in IFN signaling, interferon regulatory factor 1 (IRF1), interferon receptor (IFNAR1), interferon-induced transmembrane protein 1 (IFITM1), interferon-inducible protein 1 (IFIH1), genes-encoded proteins of JAK/STAT family (STAT2, STAT3, JAK2, JAK3, TYK2) and interferon-induced protein with tetratricopeptide repeats (IFIT1, IFIT2, and IFIT3), were found to be one of the most upregulated genes in goats with HPL. Overexpression of these genes may result in higher activation of factors involved in antiviral responses. Our results may indicate that gene expression of INFs did not necessarily correspond with the protein concentration, which was also suggested by Jarczak et al. [50].
The zinc-finger (ZNF) proteins provide a particular interest in this analysis because 23 genes encoding these proteins were differentially expressed in HPL and LPL goats. Zincfinger proteins have nucleic acid-binding domains that can serve to regulate multiple gene transcription. It has been established that a deletion variant near ZNF389 influenced SRLV proviral concentration in multiple sheep flocks [63]. The functional importance of ZNF during regulation of SRLV infection in goats is currently unknown, but our results strongly suggest that expression of these genes is dominant in response to the SRLV infection and can be associated with SRLV proviral concentration as was observed for HIV [64].
Another group of genes in which expression was dysregulated in response to the infection with SRLV and was associated with proviral concentration was the genes-encoded transmembrane proteins (TMEM; 16 genes). Recently, TMEM154 and TMEM38A genes were identified as suitable candidates for SRLV resistance in sheep [15,16]. An amino acid substitution (E/K) at position 35 of the TMEM154 was associated with the lower concentration of SRLV provirus in sheep [18,19]. In the present study, we did not find any correlation between the expression of TMEM154 and proviral load. However, our results provided evidence that other genes, such as TMEM238, TMEM223, TMEM151, TMEM147, TMEM53, were upregulated and may play an important role in the course of SRLV infection and provirus concentration.

Conclusions
In this study, we have demonstrated the changes in the transcriptome profile of goats showing high and low proviral load following infection with SRLV. A total of 1434 differentially expressed genes were involved in a variety of molecular and cellular defense mechanisms of immune response, cell cycle regulation, and cellular metabolism. Numerous genes have not been previously associated with lentiviral infection and may extend structural and/or regulatory networks implicated in the course of infection with SRLV (Supplementary Materials Figure S2). The knowledge about these genes provides the basis for further work to identify genetic markers associated with SRLV infection and provirus concentration. Such markers may be used to eliminate animals predisposed to high proviral load and limit the outcome of clinical signs and spreading of the virus.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/v13102054/s1, Figure S1: Principal component analysis (PCA) plot presented samples clustering based on similarities of their gene expression variation, Figure S2: The Summary of study designed and obtained results based on whole transcriptome profiling under proviral load factor, Table S1: The primers used for real-time PCR method, Table S2: Summary of reads quality and mapping results of RNA-seq,

Informed Consent Statement: Not applicable.
Data Availability Statement: The data sets used and/or analyzed during the current study are available from the corresponding author on reasonable request.