Next Article in Journal
SARS–CoV–2 Spike Impairs DNA Damage Repair and Inhibits V(D)J Recombination In Vitro
Next Article in Special Issue
Evaluation of the Presence of ASFV in Wolf Feces Collected from Areas in Poland with ASFV Persistence
Previous Article in Journal
Molecular and Serological Characterization of Hepatitis B Virus (HBV)-Positive Samples with Very Low or Undetectable Levels of HBV Surface Antigen
Previous Article in Special Issue
Phylogenetic Analysis of European Brown Hare Syndrome Virus Strains from Poland (1992–2004)
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Department of Biochemistry, National Veterinary Research Institute, 24-100 Pulawy, Poland
2
Department of Animal Molecular Biology, National Research Institute of Animal Production, Krakowska 1, 32-083 Kraków, Poland
3
Center for Experimental and Innovative Medicine, University of Agriculture in Krakow, Rędzina 1c, 30-248 Kraków, Poland
*
Author to whom correspondence should be addressed.
Viruses 2021, 13(10), 2054; https://doi.org/10.3390/v13102054
Submission received: 3 September 2021 / Revised: 11 October 2021 / Accepted: 11 October 2021 / Published: 13 October 2021
(This article belongs to the Special Issue State-of-the-Art Animal Virus Research in Poland)

Abstract

:
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.

1. 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.
The dynamics of the host immune response to SRLV infection are still not fully understood. Several attempts to identify host factors associated with resistance to SRLV infections have been made, and some loci were identified [11,12,13,14,15]. In particular TMEM154 gene (transmembrane protein 154), TMEM38A (transmembrane protein 38A), CCR5 (chemokine(C-C motif) receptor type 5), MHC (major histocompatibility complex), ZNF389 (zinc-finger protein 389), TLRs (toll-like receptors), APOBEC3 (apolipoprotein B editing complex 3), TRIM5 (tripartite motif protein 5 alpha), Tetherin/BST-2 (bone marrow stromal cell antigen 2) and other cytokines (interleukin 2 (IL2), interleukin 2 receptor (IL2R), tumor necrosis factor alfa (TNF-α), interleukin 4 (IL4), interleukin 8 (IL8), interleukin 6 (IL6), interleukin 16 (IL-16), interferon gamma (IFN-γ), transforming growth factor beta (TGF-β1), monocyte chemoattractant protein-1 (MCP-1), granulocyte-macrophage colony-stimulating factor (GM-CSF)) and chemokine ligands (C-C motif ligand 2 (CCL2), C-C motif ligand 5 (CCL5), C-C motif ligand 20 (CCL20)) seem to have an important role in the SRLV infection susceptibility/resistance [12,13,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33]. However, additional studies are still necessary to learn more about genes involved in SRLV immunity.
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 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.

2. Methods

2.1. 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).

2.2. 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′ TGGGAGTAGGACAAACAAATCA 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 108 to 10. 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) H2O) 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.

2.3. 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 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.

2.4. 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.

2.5. 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].

2.6. 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, CA, USA), version 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.

3. Results

3.1. 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).
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.

3.2. RNA-seq Data

3.2.1. 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).

3.2.2. 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).
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.
In the panel of 1046 DEGs differentiated uninfected animals from those with SRLV proviral load, the genes with the highest expression changes and downregulated in animals with proviral load were KITLG (KIT ligand—mast cell growth factor) 256-fold change; HHPI (hedgehog interacting protein) 234-fold change; SLC17A6 (vesicular glutamate transporter 2) 222-fold change and P2RX2 (purinergic receptor P2X 2) 186-fold change. In turn, the most upregulated genes in the blood of goats with proviral load were PDGFRB (platelet-derived growth factor receptor beta) 76-fold change; TUBA4A (tubulin alpha 4a) 67-fold change; TNNI3 (troponin I3, cardiac type) 66 FC and TMEM176B (transmembrane protein 176B) 59FC.
Moreover, the following family of genes in which expression was deregulated in response to SRLV infection were identified: the zinc-finger gene family (ZNF; 8 genes); a transmembrane protein (TMEM; 10 genes); the solute carrier family genes (SLC; 14 genes); signaling protein (11 genes), toll-like receptors (TLR; 4 genes), and tubulins (TUBB; 4 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 transcript-level 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).

3.3. Gene Ontology and Pathways Annotation

3.3.1. 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.

3.3.2. 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).
In the group of downregulated genes, most of them were associated with: regulation of signaling receptor activity (p-value < 0.0001) (53 genes, e.g., chemokine (C-C motif), ligand 2 (CCL2), chemokine (C-X-C motif), ligand 5 (CXCL5), TNF superfamily member 11 (TNFSF11), C-C motif chemokine ligand 17 (CCL17), C-X-C motif chemokine ligand 9 (CXCL9), macrophage migration inhibitory factor (MIF)); the response to toxic substances (43 genes, e.g., hemoglobin subunit mu (HBM), cholinergic receptor nicotinic beta 2 subunit (CHRNB2), LY6/PLAUR domain containing 1 (LYPD1), gonadotropin-releasing hormone 1 (GNRH1)), and NADH dehydrogenase complex assembly (21 genes, e.g., NADH: ubiquinone oxidoreductase subunit A13 (NDUFA13), NADH: ubiquinone oxidoreductase subunit S5 (NDUFS5), NADH: ubiquinone oxidoreductase core subunit S7 (NDUFS7), NADH: ubiquinone oxidoreductase subunit B9 (NDUFB9)). The specific genes belonging to detected GO were presented in Table 3 and in Supplementary Materials Table S4. In the HPL group, an increased expression for 95 genes was identified (e.g., ADAM metallopeptidase with thrombospondin type 1 motif 3 (ADAMTS3), interleukin 15 (IL15), chemerin chemokine-like receptor 1 (CMKLR1), nucleotide-binding oligomerization domain containing 2 (NOD2), C-C motif chemokine receptor 2 (CCR2), interleukin 6 receptor (IL6R), interleukin 1 alpha (IL1A)) that represented GO term cytokine production (FDR < 0.0001) (Table 4, Supplementary Materials Table S5). In addition, a high number of upregulated genes related to vesicle organization was detected (38 genes; FDR < 0.030), as well as vacuole organization (34 genes, FDR < 0.035).
To better show the interaction between the broad set of genes represented by the cytokine production GO term, the gene network was prepared according to String software. The analysis indicated that some of the identified DEGs were involved in multiple biological processes related to positive and negative regulation of cytokine production and control in response to stimulus in both adaptive and innate immunity (Figure 6). Such an approach allowed to pinpointed the upregulated genes with the highest number of interactions (Table 5).
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); actin-related 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.

3.4. 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.

4. 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. Zinc-finger 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.

5. 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, Table S3: The most significantly dysregulated genes between HPL and LPL goats, Table S4: The detail fold change values and GO term analysis of down-regulated genes in group HPL goats, Table S5: The detail fold change values and GO term analysis of upregulated genes in group HPL goats.

Author Contributions

M.O. and K.R.-M., participated in experiments, analyzed data, and wrote the manuscript. T.S., made the analysis of the data. K.P., participated in experiments. J.K., designed the study and revised the paper. All authors have read and agreed to the published version of the manuscript.

Funding

Funded by the KNOW (Leading National Research Centre) Scientific Consortium “Healthy Animal—Safe Food”, under Ministry of Science and Higher Education decision No. 05-1/KNOW2/2015.

Institutional Review Board Statement

All procedures associated with animal handling and treatments were approved by the Local Ethical Committee on Animal Testing at the University of Life Sciences in Lublin, Poland (approval number No. 37/2016).

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.

Acknowledgments

The authors are grateful to Łukasz Bocian for help in statistical analyses.

Conflicts of Interest

The authors declare that they have no conflict of interest.

Abbreviations

MVV: Maedi visna virus; CAEV: caprine arthritis encephalitis virus; SRLV: small ruminant lentiviruses; PL: proviral load; HPL: high proviral load; LPL: low proviral load; DEGs: differentially expressed genes; NGS: next-generation sequencing; PCA: principal component analysis; FC: fold change; GO: gene ontology; FDR: false discovery rate; TLR: toll-like receptor; PRRs: pattern recognition receptors; IFN: interferon; TNF: tumor necrosis factor; RNA-seq: RNA sequencing; KEGG: Kyoto Encyclopedia of Genes and Genomes; HIV: human immunodeficiency virus; qPCR: quantitative polymerase chain reaction; qRT-PCR quantitative reverse transcription polymerase chain reaction; PBLs: peripheral blood leukocytes; TMEM154: transmembrane protein 154; CCR5: chemokine(C-C motif) receptor type 5; MHC: major histocompatibility complex; ZNF389: zinc-finger protein 389; APOBEC3: polipoprotein B editing complex 3; TMEM38A: transmembrane protein 38A; TRIM5: tripartite motif protein 5 alpha; BST-2: bone marrow stromal cell; EDTA: ethylenediaminetetraacetic acid; ELISA: enzyme-linked immunosorbent assay; RIN: RNA integrity number.

References

  1. Minguijón, E.; Reina, R.; Pérez, M.; Polledo, L.; Villoria, M.; Ramírez, H.; Leginagoikoa, I.; Badiola, J.J.; García-Marín, J.F.; de Andrés, D.; et al. Small ruminant lentivirus infections and diseases. Vet. Microbiol. 2015, 181, 75–89. [Google Scholar] [CrossRef] [PubMed]
  2. Blacklaws, B.A. Small ruminant lentiviruses: Immunopathogenesis of visna-maedi and caprine arthritis and encephalitis virus. Comp. Immunol. Microbiol. Infect. Dis. 2012, 35, 259–269. [Google Scholar] [CrossRef] [PubMed]
  3. Peterhans, E.; Greenland, T.; Badiola, J.; Harkiss, G.; Bertoni, G.; Amorena, B.; Eliaszewicz, M.; Juste, R.A.; Krassnig, R.; Lafont, J.P.; et al. Routes of transmission and consequences of small ruminant lentiviruses (SRLVs) infection and eradication schemes. Vet. Res. 2004, 35, 257–274. [Google Scholar] [CrossRef] [Green Version]
  4. Martínez-Navalón, B.; Peris, C.; Gómez, E.A.; Peris, B.; Roche, M.L.; Caballero, C.; Goyena, E.; Berriatua, E. Quantitative estimation of the impact of caprine arthritis encephalitis virus infection on milk production by dairy goats. Vet. J. 2013, 197, 311–317. [Google Scholar] [CrossRef]
  5. Scheer-Czechowski, P.; Vogt, H.R.; Tontis, A.; Peterhans, E.; Zanoni, R. Pilotprojekt zur Sanierung der Maedi-Visna bei Walliser Schwarznasenschafen (Pilot project for eradicating maedi-visna in Walliser blacknose sheep). Schweiz. Arch. Tierheilkd. 2000, 142, 155–164. [Google Scholar] [PubMed]
  6. Kalogianni, A.I.; Bossis, I.; Ekateriniadou, L.V.; Gelasakis, A.I. Etiology, Epizootiology and Control of Maedi-Visna in Dairy Sheep: A Review. Animals 2020, 10, 616. [Google Scholar] [CrossRef] [Green Version]
  7. Tavella, A.; Bettini, A.; Ceol, M.; Zambotto, P.; Stifter, E.; Kusstatscher, N.; Lombardi, R.; Nardeli, S.; Beato, M.S.; Capello, K.; et al. Achievements of an eradication programme against caprine arthritis encephalitis virus in South Tyrol, Italy. Vet. Rec. 2018, 182, 51. [Google Scholar] [CrossRef] [Green Version]
  8. Reina, R.; Berriatua, E.; Luján, L.; Juste, R.; Sánchez, A.; de Andrés, D.; Amorena, B. Prevention strategies against small ruminant lentiviruses: An update. Vet. J. 2009, 182, 31–37. [Google Scholar] [CrossRef]
  9. Cirone, F.; Maggiolino, A.; Cirilli, M.; Sposato, A.; De Palo, P.; Ciappetta, G.; Pratelli, A. Small ruminant lentiviruses in goats in southern Italy: Serological evidence, risk factors and implementation of control programs. Vet. Microbiol. 2019, 228, 143–146. [Google Scholar] [CrossRef]
  10. Gjerset, B.; Rimstad, E.; Teige, J.; Soetaert, K.; Jonassen, C.M. Impact of natural sheep-goat transmission on detection and control of small ruminant lentivirus group C infections. Vet. Microbiol. 2009, 135, 231–238. [Google Scholar] [CrossRef]
  11. Herrmann-Hoesing, L.M.; White, S.N.; Mousel, M.R.; Lewis, G.S.; Knowles, D.P. Ovine progressive pneumonia provirus levels associate with breed and Ovar-DRB1. Immunogenetics 2008, 60, 749–758. [Google Scholar] [CrossRef]
  12. Larruskain, A.; Bernales, I.; Luján, L.; de Andrés, D.; Amorena, B.; Jugo, B.M. Expression analysis of 13 ovine immune response candidate genes in Visna/Maedi disease progression. Comp. Immunol. Microbiol. Infect. Dis. 2013, 36, 405–413. [Google Scholar] [CrossRef] [Green Version]
  13. White, S.N.; Knowles, D.P. Expanding possibilities for intervention against small ruminant lentiviruses through genetic marker-assisted selective breeding. Viruses 2013, 5, 1466–1499. [Google Scholar] [CrossRef]
  14. Cecchi, F.; Dadousis, C.; Bozzi, R.; Fratini, F.; Russo, C.; Bandecchi, P.; Cantile, C.; Mazzei, M. Genome scan for the possibility of identifying candidate resistance genes for goat lentiviral infections in the Italian Garfagnina goat breed. Trop. Anim. Health Prod. 2019, 51, 729–733. [Google Scholar] [CrossRef]
  15. White, S.N.; Mousel, M.R.; Herrmann-Hoesing, L.M.; Reynolds, J.O.; Leymaster, K.A.; Neibergs, H.L.; Lewis, G.S.; Knowles, D.P. Genome-wide association identifies multiple genomic regions associated with susceptibility to and control of ovine lentivirus. PLoS ONE 2012, 7, e47829. [Google Scholar] [CrossRef] [Green Version]
  16. Alshanbari, F.A.; Mousel, M.R.; Reynolds, J.O.; Herrmann-Hoesing, L.M.; Highland, M.A.; Lewis, G.S.; White, S.N. Mutations in Ovis aries TMEM154 are associated with lower small ruminant lentivirus proviral concentration in one sheep flock. Anim. Genet. 2014, 45, 565–571. [Google Scholar] [CrossRef]
  17. Colussi, S.; Desiato, R.; Beltramo, C.; Peletto, S.; Modesto, P.; Maniaci, M.G.; Campia, V.; Quasso, A.; Rosati, S.; Bertolotti, L.; et al. A single nucleotide variant in the promoter region of the CCR5 gene increases susceptibility to arthritis encephalitis virus in goats. BMC Vet. Res. 2019, 15, 230. [Google Scholar] [CrossRef] [PubMed]
  18. Molaee, V.; Eltanany, M.; Lühken, G. First survey on association of TMEM154 and CCR5 variants with serological maedi-visna status of sheep in German flocks. Vet. Res. 2018, 49, 36. [Google Scholar] [CrossRef] [Green Version]
  19. Molaee, V.; Otarod, V.; Abdollahi, D.; Lühken, G. Lentivirus Susceptibility in Iranian and German Sheep Assessed by Determination of TMEM154 E35K. Animals 2019, 9, 685. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  20. Jarczak, J.; Kaba, J.; Reczyńska, D.; Bagnicka, E. Impaired Expression of Cytokines as a Result of Viral Infections with an Emphasis on Small Ruminant Lentivirus Infection in Goats. Viruses 2016, 8, 186. [Google Scholar] [CrossRef] [PubMed]
  21. Thompson, J.; Ma, F.; Quinn, M.; Xiang, S.H. Genome-Wide Search for Host Association Factors during Ovine Progressive Pneumonia Virus Infection. PLoS ONE 2016, 11, e0150344. [Google Scholar] [CrossRef] [Green Version]
  22. Ramírez, H.; Echeverría, I.; Benito, A.A.; Glaria, I.; Benavides, J.; Pérez, V.; de Andrés, D.; Reina, R. Accurate Diagnosis of Small Ruminant Lentivirus Infection Is Needed for Selection of Resistant Sheep through TMEM154 E35K Genotyping. Pathogens 2021, 10, 83. [Google Scholar] [CrossRef]
  23. Jáuregui, P.; Crespo, H.; Glaria, I.; Luján, L.; Contreras, A.; Rosati, S.; de Andrés, D.; Amorena, B.; Towers, G.J.; Reina, R. Ovine TRIM5α can restrict visna/maedi virus. J. Virol. 2012, 86, 9504–9509. [Google Scholar] [CrossRef] [Green Version]
  24. White, S.N.; Mousel, M.R.; Reynolds, J.O.; Lewis, G.S.; Herrmann-Hoesing, L.M. Common promoter deletion is associated with 3.9-fold differential transcription of ovine CCR5 and reduced proviral level of ovine progressive pneumonia virus. Anim. Genet. 2009, 40, 583–589. [Google Scholar] [CrossRef]
  25. Bowles, D.; Carson, A.; Isaac, P. Genetic Distinctiveness of the Herdwick Sheep Breed and Two Other Locally Adapted Hill Breeds of the UK. PLoS ONE 2014, 9, e87823. [Google Scholar] [CrossRef] [PubMed]
  26. Heaton, M.P.; Clawson, M.L.; Chitko-Mckown, C.G.; Leymaster, K.A.; Smith, T.P.; Harhay, G.P.; White, S.N.; Herrmann-Hoesing, L.M.; Mousel, M.R.; Lewis, G.S.; et al. Reduced lentivirus susceptibility in sheep with TMEM154 mutations. PLoS Genet. 2012, 8, e1002467. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  27. Heaton, M.P.; Kalbfleisch, T.S.; Petrik, D.T.; Simpson, B.; Kijas, J.W.; Clawson, M.L.; Chitko-McKown, C.G.; Harhay, G.P.; Leymaster, K.A.; International Sheep Genomics Consortium. Genetic testing for TMEM154 mutations associated with lentivirus susceptibility in sheep. PLoS ONE 2013, 8, e55490. [Google Scholar] [CrossRef] [Green Version]
  28. Clawson, M.L.; Redden, R.; Schuller, G.; Heaton, M.P.; Workman, A.; Chitko-McKown, C.G.; Smith, T.P.; Leymaster, K.A. Genetic subgroup of small ruminant lentiviruses that infects sheep homozygous for TMEM154 frameshift deletion mutation A4Δ53. Vet. Res. 2015, 46, 22. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Sider, L.H.; Heaton, M.P.; Chitko-McKown, C.G.; Harhay, G.P.; Smith, T.P.; Leymaster, K.A.; Laegreid, W.W.; Clawson, M.L. Small ruminant lentivirus genetic subgroups associate with sheep TMEM154 genotypes. Vet. Res. 2013, 44, 64. [Google Scholar] [CrossRef] [Green Version]
  30. Larruskain, A.; Minguijón, E.; García-Etxebarria, K.; Moreno, B.; Arostegui, I.; Juste, R.A.; Jugo, B.M. MHC class II DRB1 gene polymorphism in the pathogenesis of Maedi-Visna and pulmonary adenocarcinoma viral diseases in sheep. Immunogenetics 2010, 62, 75–83. [Google Scholar] [CrossRef]
  31. Yaman, Y.; Keleş, M.; Aymaz, R.; Sevim, S.; Sezenler, T.; Önaldı, A.T.; Kaptan, C.; Baskurt, A.; Koncagül, S.; Öner, Y.; et al. Association of TMEM154 variants with visna/maedi virus infection in Turkish sheep. Small Rumin. Res. 2019, 177, 61–67. [Google Scholar] [CrossRef]
  32. Leymaster, K.A.; Chitko-McKown, C.G.; Clawson, M.L.; Harhay, G.P.; Heaton, M.P. Effects of TMEM154 haplotypes 1 and 3 on susceptibility to ovine progressive pneumonia virus following natural exposure in sheep. J. Anim. Sci. 2013, 91, 5114–5121. [Google Scholar] [CrossRef]
  33. Abendaño, N.; Esparza-Baquer, A.; Bernales, I.; Reina, R.; de Andrés, D.; Jugo, B.M. Gene Expression Profiling Reveals New Pathways and Genes Associated with Visna/Maedi Viral Disease. Animals 2021, 11, 1785. [Google Scholar] [CrossRef] [PubMed]
  34. Panei, C.J.; Takeshima, S.N.; Omori, T.; Nunoya, T.; Davis, W.C.; Ishizaki, H.; Matoba, K.; Aida, Y. Estimation of bovine leukemia virus (BLV) proviral load harbored by lymphocyte subpopulations in BLV-infected cattle at the subclinical stage of enzootic bovine leucosis using BLV-CoCoMo-qPCR. BMC Vet. Res. 2013, 9, 95. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Ravazzolo, A.P.; Nenci, C.; Vogt, H.R.; Waldvogel, A.; Obexer-Ruff, G.; Peterhans, E.; Bertoni, G. Viral load, organ distribution, histopathological lesions, and cytokine mRNA expression in goats infected with a molecular clone of the caprine arthritis encephalitis virus. Virology 2006, 350, 116–127. [Google Scholar] [CrossRef] [PubMed]
  36. Herrmann-Hoesing, L.M.; Noh, S.M.; White, S.N.; Snekvik, K.R.; Truscott, T.; Knowles, D.P. Peripheral ovine progressive pneumonia provirus levels correlate with and predict histological tissue lesion severity in naturally infected sheep. Clin. Vaccine Immunol. 2009, 16, 551–557. [Google Scholar] [CrossRef] [Green Version]
  37. Crespo, H.; Bertolotti, L.; Proffiti, M.; Cascio, P.; Cerruti, F.; Acutis, P.L.; de Andrés, D.; Reina, R.; Rosati, S. Low proviral small ruminant lentivirus load as biomarker of natural restriction in goats. Vet. Microbiol. 2016, 192, 152–162. [Google Scholar] [CrossRef]
  38. Stonos, N.; Wootton, S.K.; Karrow, N. Immunogenetics of small ruminant lentiviral infections. Viruses 2014, 6, 3311–3333. [Google Scholar] [CrossRef]
  39. Olech, M.; Kuźmak, J. Molecular Characterization of Small Ruminant Lentiviruses of Subtype A5 Detected in Naturally Infected but Clinically Healthy Goats of Carpathian Breed. Pathogens 2020, 9, 992. [Google Scholar] [CrossRef]
  40. Andrews, S. FastQC: A Quality Control Tool for High Throughput Sequence Data. 2010. Available online: http://www.bioinformatics.babraham.ac.uk/projects/fastqc (accessed on 22 June 2020).
  41. Dodt, M.; Roehr, J.T.; Ahmed, R.; Dieterich, C. FLEXBAR-Flexible Barcode and Adapter Processing for Next-Generation Sequencing Platforms. Biology 2012, 1, 895–905. [Google Scholar] [CrossRef] [Green Version]
  42. Trapnell, C.; Pachter, L.; Salzberg, S.L. TopHat: Discovering splice junctions with RNA-Seq. Bioinformatics 2009, 25, 1105–1111. [Google Scholar] [CrossRef]
  43. Anders, S.; Pyl, P.T.; Huber, W. HTSeq—A Python framework to work with high-throughput sequencing data. Bioinformatics 2015, 31, 166–169. [Google Scholar] [CrossRef] [PubMed]
  44. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Huang da, W.; Sherman, B.T.; Lempicki, R.A. Bioinformatics enrichment tools: Paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009, 37, 1–13. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  46. Pfaffl, M.W. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001, 29, e45. [Google Scholar] [CrossRef]
  47. Sahu, A.R.; Wani, S.A.; Saxena, S.; Rajak, K.K.; Chaudhary, D.; Sahoo, A.P.; Khanduri, A.; Pandey, A.; Mondal, P.; Malla, W.A.; et al. Selection and validation of suitable reference genes for qPCR gene expression analysis in goats and sheep under Peste des petits ruminants virus (PPRV), lineage IV infection. Sci. Rep. 2018, 8, 15969. [Google Scholar] [CrossRef] [Green Version]
  48. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2013; Available online: http://www.R-project.org/ (accessed on 22 June 2020).
  49. Swaminathan, S.; Qiu, J.; Rupert, A.W.; Hu, Z.; Higgins, J.; Dewar, R.L.; Stevens, R.; Rehm, C.A.; Metcalf, J.A.; Sherman, B.T.; et al. Interleukin-15 (IL-15) Strongly Correlates with Increasing HIV-1 Viremia and Markers of Inflammation. PLoS ONE 2016, 11, e0167091. [Google Scholar] [CrossRef]
  50. Jarczak, J.; Słoniewska, D.; Kaba, J.; Bagnicka, E. The expression of cytokines in the milk somatic cells, blood leukocytes and serum of goats infected with small ruminant lentivirus. BMC Vet. Res. 2019, 15, 424. [Google Scholar] [CrossRef]
  51. Di Paolo, N.C.; Shayakhmetov, D.M. Interleukin 1α and the inflammatory process. Nat. Immunol. 2016, 17, 906–913. [Google Scholar] [CrossRef] [Green Version]
  52. Gabay, C.; Lamacchia, C.; Palmer, G. IL-1 pathways in inflammation and human diseases. Nat. Rev. Rheumatol 2010, 6, 232–241. [Google Scholar] [CrossRef]
  53. Zhang, Z.; Watt, N.J.; Hopkins, J.; Harkiss, G.; Woodall, C.J. Quantitative analysis of maedi-visna virus DNA load in peripheral blood monocytes and alveolar macrophages. J. Virol. Methods. 2000, 86, 13–20. [Google Scholar] [CrossRef]
  54. Mikula, I.; Bhide, M.; Pastorekova, S.; Mikula, I. Characterization of ovine TLR7 and TLR8 protein coding regions, detection of mutations and Maedi Visna virus infection. Vet. Immunol. Immunopathol. 2010, 138, 51–59. [Google Scholar] [CrossRef]
  55. Olech, M.; Ropka-Molik, K.; Szmatoła, T.; Piórkowska, K.; Kuźmak, J. Single Nucleotide Polymorphisms in Genes Encoding Toll-Like Receptors 7 and 8 and Their Association with Proviral Load of SRLVs in Goats of Polish Carpathian Breed. Animals 2021, 11, 1908. [Google Scholar] [CrossRef]
  56. Ge, Y.; Mansell, A.; Ussher, J.E.; Brooks, A.E.; Manning, K.; Wang, C.J.; Taylor, J.A. Rotavirus NSP4 Triggers Secretion of Proinflammatory Cytokines from Macrophages via Toll-Like Receptor 2. J. Virol. 2013, 87, 11160–11167. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  57. Lester, S.N.; Li, K. Toll-like receptors in antiviral innate immunity. J. Mol. Biol. 2014, 426, 1246–1264. [Google Scholar] [CrossRef] [PubMed]
  58. Oliveira-Nascimento, L.; Massari, P.; Wetzler, L.M. The Role of TLR2 in Infection and Immunity. Front. Immunol. 2012, 3, 79. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  59. Hernández, J.C.; Stevenson, M.; Latz, E.; Urcuqui-Inchima, S. HIV type 1 infection upregulates TLR2 and TLR4 expression and function in vivo and in vitro. AIDS Res. Hum. Retroviruses 2012, 28, 1313–1328. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  60. Lechner, F.; Vogt, H.R.; Seow, H.F.; Bertoni, G.; Cheevers, W.P.; von Bodungen, U.; Zurbriggen, A.; Peterhans, E. Expression of cytokine mRNA in lentivirus-induced arthritis. Am. J. Pathol. 1997, 151, 1053–1065. [Google Scholar]
  61. Lechner, F.; Machado, J.; Bertoni, G.; Seow, H.F.; Dobbelaere, D.A.; Peterhans, E. Caprine arthritis encephalitis virus dysregulates the expression of cytokines in macrophages. J. Virol. 1997, 71, 7488–7497. [Google Scholar] [CrossRef] [Green Version]
  62. Gómez-Lucía, E.; Collado, V.M.; Miró, G.; Doménech, A. Effect of type-I interferon on retroviruses. Viruses 2009, 1, 545–573. [Google Scholar] [CrossRef] [Green Version]
  63. White, S.N.; Mousel, M.R.; Reynolds, J.O.; Herrmann-Hoesing, L.M.; Knowles, D.P. Deletion variant near ZNF389 is associated with control of ovine lentivirus in multiple sheep flocks. Anim. Genet. 2014, 45, 297–300. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  64. Zhu, Y.; Chen, G.; Lv, F.; Wang, X.; Ji, X.; Xu, Y.; Sun, J.; Wu, L.; Zheng, Y.T.; Gao, G. Zinc-finger antiviral protein inhibits HIV-1 infection by selectively targeting multiply spliced viral mRNAs for degradation. Proc. Natl. Acad. Sci. USA 2011, 108, 15834–15839. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Figure 1. Proviral load distribution among goats with low (LPL) and high (HPL) proviral load. HPL—high proviral load, LPL—low proviral load.
Figure 1. Proviral load distribution among goats with low (LPL) and high (HPL) proviral load. HPL—high proviral load, LPL—low proviral load.
Viruses 13 02054 g001
Figure 2. Venn diagram showing the number of overlapping DEGs between C vs. LPL, C vs. HPL, and LPL vs. HPL.
Figure 2. Venn diagram showing the number of overlapping DEGs between C vs. LPL, C vs. HPL, and LPL vs. HPL.
Viruses 13 02054 g002
Figure 3. The heatmap of the most significantly deregulated genes between the uninfected goats and goats with high (HPL) and low (LPL) proviral loads. Animals with high (HPL) proviral load: b15, b25, b26 and b30; animals with low (LPL) proviral load: a8, a31, a32, and a33; Uninfected animals: c44, c45, and c46.
Figure 3. The heatmap of the most significantly deregulated genes between the uninfected goats and goats with high (HPL) and low (LPL) proviral loads. Animals with high (HPL) proviral load: b15, b25, b26 and b30; animals with low (LPL) proviral load: a8, a31, a32, and a33; Uninfected animals: c44, c45, and c46.
Viruses 13 02054 g003
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).
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).
Viruses 13 02054 g004
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.
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.
Viruses 13 02054 g005
Figure 6. Tree graph represented the top 20 (10 up and 10 down) most regulated GO terms. Orange and blue color represented down and upregulated GO term, respectively. (WebGestalt, WEB-based GEne SeT AnaLysis Toolkit).
Figure 6. Tree graph represented the top 20 (10 up and 10 down) most regulated GO terms. Orange and blue color represented down and upregulated GO term, respectively. (WebGestalt, WEB-based GEne SeT AnaLysis Toolkit).
Viruses 13 02054 g006
Table 1. The significant GO terms enrichment between uninfected and infected goats.
Table 1. The significant GO terms enrichment between uninfected and infected goats.
GOAccession NumberNumber of GenesFDRIdentified Genes
MyD88-dependent toll-like receptor signaling pathwayGO: 000275540.019TLR2, TLR4, TLR8, TLR7
Response to growth factorGO: 007084850.050HHIP, BMPR1B, ADAMTS3, PDGFRB, BMP6
Regulation of cytokine secretionGO: 000181730.050TLR2, TLR4, TLR8
Inflammatory responseGO: 000695440.050TLR2, TLR4, TLR8, TLR7
Transmembrane signaling receptor activityGO: 000488830.010TLR2, TLR4, TLR7
Innate immune responseGO: 004508730.010TLR2, TLR4, TLR8
Extracellular regionGO: 004508760.010KITLG, BMP6, INHBB, INSL3, IGFBP1, IGFBP3
Integral component of membraneGO: 0016021130.010ATP6, KITLG, ND3, ND4, BMPR1B, CALCRL, COX1, DGAT2, SLC11A1, SCD, TLR2, TLR8, TLR7
FDR—false discovery rate (p-value adjusted), GO—gene ontology.
Table 2. Significantly overrepresented pathways involved genes associated with SRLV infection.
Table 2. Significantly overrepresented pathways involved genes associated with SRLV infection.
Biological PathwaysNumber of Genes UpregulatedUpregulated GenesNumber of Genes DownregulatedDownregulated GenesFDR
Toll-like receptor signaling pathway6LY96, PIK3R5, TLR2, TLR4, TLR7, TLR81MAPK120.050
Rheumatoid arthritis4ATP6V1A, TLR2, TLR4, TGFB20-0.010
Ras signaling pathway5GAB2, EFNA4, KDR, PIK3R5, PDGFRB5KITLG, ANGPT2, PAK6, PLA2G1B, RGL10.010
PI3K-Akt signaling pathway8CSF3R, EFNA4, GYS1, KDR, PIK3R5, PDGFRB, TLR2, TLR44KITLG, ANGPT2, COMP, COL1A10.010
TNF signaling pathway3TNFRSF1A, PIK3R5, SOCS31MAPK120.010
Phagosome9ATP6V1A, RAB7B, LAMP2, MSR1, NCF1, TLR2, TLR4, TUBB1, TUBB4A1COMP0.020
Cytokine-cytokine receptor interaction9TNFRSF1A, CSF2RB, CSF3R, CRLF2, IFNGR1, IFNGR2, IL1R2, PPBP, TGFB25CCR10, TNFRSF13C, TNFRSF6B, AMH, BMPR1B, 0.004
False discovery rate (p-value adjusted).
Table 3. The significant GO terms of downregulated genes in group HPL goats.
Table 3. The significant GO terms of downregulated genes in group HPL goats.
GOAccession NumberNumber of GenesFDRIdentified Genes
execution phase of apoptosisGO: 0097194100.0240PTGIS, ENDOG, BOK, ERN2, SHARPIN, CIDEC, CXCR3, SIRT2, RPS3, BAX
tetrapyrrole metabolic processGO: 003301390.00679ALAS2, ALAD, UROD, CYP1A2, UROS, HNF1A, MMAB, ATP5IF1, BLVRB
response to toxic substanceGO: 0009636430.00691HBM, CHRNB2, LYPD1, GNRH1, MPO, CHRNA6, DRD3, CHRND, LTC4S, IL6
regulation of signaling receptor activityGO: 001046953<0.0000FGF16, CCL2, CGA, INHBE, CLEC11A, CXCL5, TNFSF11, AVP, LYPD1, FOXH1, GNRH1, NPY, CCL17, OXT, GHRH, NOG, EDA, CXCL9, MIF, RETN
peptidyl-arginine modificationGO: 0018195100.01500ART1, PADI6, PADI3, PADI1, KRTCAP2, COPRS, PARK7, PRMT7, PRMT2, PRMT1
amine transportGO: 0015837150.0060CHRNB2, CHRNA6, DRD3, ACE2, SNCG, SYT2, SLC22A16, SYT1, SLC18A1, DTNBP1
organic cation transportGO: 001569530.0290SLC22A1, SLC22A16, SLC18A3
ammonium transportGO: 001569615<0.0000RHAG, CHRNB2, CHRNA6, DRD3, SLC6A2, SNCG, SYT2, ADCYAP1, SLC22A16, SYT1
NADH dehydrogenase complex assemblyGO: 0010257210.0010NDUFA13, NDUFS5, NDUFS7, NDUFB9, NDUFA9, NDUFA8, BCS1L, NDUFB2, NDUFA5, NDUFA2
FDR—false discovery rate (p-value adjusted), GO—gene ontology, HPL—high proviral load.
Table 4. The significant GO terms of upregulated genes in group HPL goats.
Table 4. The significant GO terms of upregulated genes in group HPL goats.
GOAccession NumberNumber of GenesFDRGenes
peptidyl-cysteine modificationGO: 001819870.0020GOLGA7B, MAP6D1, TMX3, RAB3D, ZDHHC23, RAB6A, ZDHHC22
cytokine productionGO: 000181695<0.000ADAMTS3, IL15, CMKLR1, NOD2, ABCA1, GDF2, CCR2, HFE, RSAD2, CHUK, ITK, LBP, ADCY7, IL6R, IL1A, RAB7B, IFIH1, C3AR1, LPL, CD274
vacuole organizationGO: 0007033340.0351ABCA1, PIP4K2B, UBXN2A, RALB, RAB7B, PPT1, GAA, ACP2, TBC1D14, AKTIP
vesicle organizationGO: 0016050380.0300SEC16B, RAB8B, ABCA1, SAMD9, VPS39, RAB7B, DYSF, STX19, EXOC8, BCL2
glutamate receptor signaling pathwayGO: 000721550.0600CACNG3, GRM2, CRHBP, HOMER2, PRNP
tissue remodelingGO: 0048771140.0600IHH, DCSTAMP, RASSF2, CCR2, DLL4, RAB3D, IL1A, EFNA2, NOTCH2, CLDN18
FDR—false discovery rate (p-value adjusted), GO—gene ontology, HPL—high proviral load.
Table 5. The selected genes in goats that were significantly upregulated with a high number of provirus copies (HPL) involved in multiple biological processes of cytokine production (GO: 0001816; FDR < 0.000).
Table 5. The selected genes in goats that were significantly upregulated with a high number of provirus copies (HPL) involved in multiple biological processes of cytokine production (GO: 0001816; FDR < 0.000).
GeneProtein NameFCAdjpvalProtein Function
TLR4Toll-like receptor 41.490.02Acts via MYD88, TIRAP, and TRAF6, leading to NF-kappa-B activation, cytokine secretion, and the inflammatory response.
TLR2Toll-like receptor 21.600.03Related to mediating the innate immune response to bacterial lipoproteins or lipopeptides, related to cytokine secretion and the inflammatory response.
TLR6Toll-like receptor 61.660.04Acts via MYD88 and TRAF6, leading to NF-kappa-B activation, cytokine secretion, and the inflammatory response.
CHUKInhibitor of nuclear factor kappa-B kinase subunit alpha2.440.03Plays an essential role in the NF- kappa-B signaling pathway activated by multiple stimuli also by viral products.
CSF1RMacrophage colony-stimulating factor 1 receptor1.800.03Controlling the proliferation and differentiation of hematopoietic precursor cells, especially mononuclear phagocytes, such as macrophages and monocytes.
IRF1Interferon regulatory factor 11.410.02Regulation of IFN and IFN-inducible genes, host response to viral and bacterial infections.
NRLP3NACHT, LRR, and PYD domain-containing protein 3 Plays a crucial role in innate immunity and inflammation.
IFIH1Interferon-induced helicase C domain-containing protein 12.140.05Plays 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.
TBK1Serine/threonine-protein kinase TBK11.800.02Regulation of transcriptional activation of proinflammatory and antiviral genes including IFNA and IFNB.
CD14Monocyte differentiation antigen CD141.670.05Acts via MyD88, TIRAP, and TRAF6, leading to NF-kappa-B activation, cytokine secretion, and the inflammatory response.
MYD88Myeloid differentiation primary response protein MyD881.820.01Acts via toll-like receptor and IL-1 receptor signaling pathway in the innate immune response.
FC—fold change.
Table 6. Significantly overrepresented pathways involved genes for which expressions were associated with SRLV copy numbers.
Table 6. Significantly overrepresented pathways involved genes for which expressions were associated with SRLV copy numbers.
Biological PathwaysNumber of Genes UpregulatedUpregulated GenesNumber of Genes DownregulatedDownregulated GenesadjP *
B-cell receptor signaling pathway8LYN, CHUK, DAPP1, GRB2, PIK3CA, PIK3CB, PIK3AP1, PIK3R14NFKBIB, HRAS, CD79B, CD79A0.018
Fc gamma R-mediated phagocytosis7CRKL, LYN, ARPC5, PIK3CA, PIK3CB, PIK3R1, PRKCD5DOCK2, LAT, PLPP2, PRKCG, RPS6KB20.034
Apoptosis6FAS, CHUK, PIK3CA, PIK3CB, PIK3R1, TNFSF104AIFM1, ENDOG, IL3RA, NTRK10.034
Natural killer cell-mediated cytotoxicity8FAS, GRB2, IFNAR1, PIK3CA, PIK3CB, PIK3R1, PTK2B, TNFSF104HRAS, HCST, LAT, PRKCG0.057
T-cell receptor signaling pathway7CHUK, GRB2, MAPK14, PIK3CA, PIK3CB, PIK3R1, TEC3HRAS, NFKBIB, LAT0.036
mTOR signaling pathway5EIF4E, PIK3CA, PIK3CB, PIK3R1, ULK22PRKCG, RPS6KB20.037
FoxO signaling pathway13EP300, CHUK, CDKN1A, GADD45A, GRB2, MAPK14, PIK3CA, PIK3CB, PIK3R1, PRKAB1, STAT3, TGFBR2, TNFSF103HRAS, FOXO1, G6PC30.034
TNF signaling pathway9FAS, TNFRSF1B, CHUK, IL15, MAPK14, NOD2, PIK3CA, PIK3CB, PIK3R10-0.044
Toll-like receptor signaling pathway8TBK1, CHUK, IFNAR1, MAPK14, PIK3CA, PIK3CB, PIK3R1, TLR60-0.044
Regulation of actin cytoskeleton10CRKL, GIT1, IQGAP1, ARPC5, CYFIP1, PXN, PIK3CA, PIK3CB, PIK3R1, PPP1CB4HRAS, FGFR2, PFN1, PPP1CA0.017
HIF-1 signaling pathway8EP300, CUL2, CDKN1A, EIF4E, PIK3CA, PIK3CB, PIK3R1, STAT34TIMP1, FLT1, PRKCG, RPS6KB20.054
Inflammatory mediator regulation of TRP channels8ADCY7, MAPK14, PIK3CA, PIK3CB, PIK3R1, PLA2G4A, PRKCD, PPP1CB4CALM3, NTRK1, PRKCG, PPP1CA0.050
Signaling pathways regulating pluripotency of stem cells8JAK1, JAK2, GRB2, MAPK14, PIK3CA, PIK3CB, PIK3R1, STAT32HRAS, FGFR20.049
Ras signaling pathway11RAP1B, RALB, TBK1, CSF1R, CHUK, GRB2, PIK3CA, PIK3CB, PIK3R1, PLA2G4A, RALBP16HRAS, CALM3, FGFR2, FLT1, LAT, PRKCG0.033
P * value adjusted using Benjamini–Hochberg correction.
Table 7. The correlation coefficients obtained for RNA-seq data validation using qPCR and between qPCR results and provirus copy number.
Table 7. The correlation coefficients obtained for RNA-seq data validation using qPCR and between qPCR results and provirus copy number.
Correlation Coefficient
Gene SymbolqPCR vs. RNA-seq 1qPCR vs. Provirus Copy Number 2
CCL20.8500.778 ***
IL150.919 **−0.242 ns
CXCR30.794 *0.359 ns
MIF−0.248 ns−0.270 ns
NOD20.433 *−0.470 *
CCR0.441 ns0.515 **
BCL2−0.107 ns0.673 **
CXCL5−0.232 ns0.759 ***
ITK−0.589 ns0.478 *
1 correlation between qPCR and RNA-seq data; 2 correlation between qPCR data estimated for all goats tested by qPCR and provirus copy number; * p-value < 0.05; ** p-value < 0.01; *** p-value < 0.0001, −p-value < 0.1, ns—not significant.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Olech, M.; Ropka-Molik, K.; Szmatoła, T.; Piórkowska, K.; Kuźmak, J. Transcriptome Analysis for Genes Associated with Small Ruminant Lentiviruses Infection in Goats of Carpathian Breed. Viruses 2021, 13, 2054. https://doi.org/10.3390/v13102054

AMA Style

Olech M, Ropka-Molik K, Szmatoła T, Piórkowska K, Kuźmak J. Transcriptome Analysis for Genes Associated with Small Ruminant Lentiviruses Infection in Goats of Carpathian Breed. Viruses. 2021; 13(10):2054. https://doi.org/10.3390/v13102054

Chicago/Turabian Style

Olech, Monika, Katarzyna Ropka-Molik, Tomasz Szmatoła, Katarzyna Piórkowska, and Jacek Kuźmak. 2021. "Transcriptome Analysis for Genes Associated with Small Ruminant Lentiviruses Infection in Goats of Carpathian Breed" Viruses 13, no. 10: 2054. https://doi.org/10.3390/v13102054

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop