ABSTRACT: Purpose: RNA-sequencing (RNA-seq) was used to identify the changes in gene expression profile in horse skin fibroblast cell line expressing one of the BPV genes encoding transforming proteins (BPV-E4 and BPV-E4^E1). The whole transcriptome sequencing allows us to establish the influence of BPV gene transfection on host cells and indicates which virus genes can be responsible for the neoplastic process. It is the first study concerning a cell model of horse skin transfected with BPV genes. Such a model could contribute to a more accurate understanding of changes inside cells during viral infection. Methods: The equine fibroblast cell lines were transfected using the nucleofection method with BPV-E4 and BPV-E4^E1 constructs. RNA was isolated from cell lines directly from culture dishes with PureLink™ RNA mini kit. The quality and quantity of obtained libraries (the TruSeq RNA Kit v2 kit (Illumina, San Diego, CA, USA) ) were assessed using Qubit 2.0 (Qubit™ dsDNA BR AssayKit, Invitrogen, Waltham, MA, USA ) and TapeStation 2200 (D100 screencaps, Agilent). Next, the cDNA libraries were sequenced on NextSeq 500 Illumina platform (Illumina, San Diego, CA, USA) and NextSeq 500/550 High Output KIT v 2.5 (75 cycles) according to protocol. The qPCR validation was used to confirm RNA-seq data. Then, the quality of raw reads was checked with FastQC software followed by removal of adapters and reads of low quality (under phred quality of 20) and reads under length of 36 (Flexbar software). Next, the filtered reads were mapped to EquCab3 genome with STAR software and reads were counted to specific gene thresholds provided in Ensembl gtf file version 100 by htseq-count software. Differential expresion extimation was followed with the use of Deseq2 software. Results: The transcriptome profiling allowed us to perform a comparison of the whole expression profile between the control and both BPV-E4 and BPV-E4^E1 groups. According to the comparison of control and BPV-E4 groups, the 1640 DEGs were identified, of which 624 were up-regulated and 1016 down-regulated in the BPV-E4 group. The highest number of DEGs – 3328 were detected from the control and BPV-E4^E1 comparison. Among them, 1626 genes were up-regulated and 1602 down-regulated in the BPV-E4^E1 group compared to control samples. The most overrepresented Gene Ontology Term between control and BPV-E4 groups were negative regulation of cell proliferation – 34 DEGs (FDR<0.002); positive regulation of cell migration – 24 DEGs (FDR<0.0001); cell adhesion and cell migration; 21 DEGs in both GO terms (FDR<0.003 and 0.0005; respectively). The enrichment analysis performed for differentially expressed genes between control and BPV-E4^E1 samples showed the focal adhesion GO terms were overrepresented by 99 DEGs of which 65 were up-regulated and 34 down-regulated (FDR<0.0001). For cells transfected BPV-E4 and BPV-E4^E1 constructs, the most significant were pathways: of the cell cycle, cytoskeleton, and ECM-matrix remodelling; regulation of actin cytoskeleton; focal adhesion and ECM-receptor interaction. The Pathways in cancer were identified uniquely for cells transfected by BPV-E4 construct, while only for BPV-E4^E1 cells the FoxO, Rap1, and TNF-signaling path-ways and Proteoglycans in cancer. Conclusions: The results obtained showed that both transfection types significantly affected cells transcriptome towards cancer transformation. Nevertheless, cells transfected with BPV-E4^E1 genes seem to be closer to the molecular modification that occurred in vivo in equine sarcoids. The present study showed how equine fibroblast cells could be modified, at the molecular level, in the presence of BPV E4 gene. These findings broaden the knowledge about the possible interaction of BPV viruses on host cells.
