The checkpoint of time for sampling in RNA-seq
First of all, the optimal time points for extracting total RNA to perform RNA-seq were identified. Innate immunity can be activated in the early stage of infection (≤ 2 h). [20]. The process of viral DNA replication represents major target for host antiviral activity, which is a critical bottleneck for successful infection [21]. Poxvirus DNA replication starts at approximately 2 hpi [4]. Therefore, to investigate the LSDV-induced innate immunity response, one of the checkpoints in RNA-seq was set at 4 hpi, corresponding to the early stage of infection. Typical pox-like lesions in LSDV-infected MDBK cells (MOI = 0.5) were observed at 24 hpi, characterized by cell clustering and formation of strip-like structures (Fig. 1D). These lesions worsened by 48 hpi, with cells beginning to detach and crumble by 72 hpi (Fig. 1E and F). No lesion was observed in PBS-treated cells during the monitoring period (Fig. 1A, B, and C), and there were also no observed lesions in LSDV-infected cells at 4 hours post-infection (data not shown). The virus titers showed an increasing trend over time (up to 72 h) (Fig. 1G). The p35 protein, which is expressed in the late stage of the viral replication cycle, was first detectable at 24 hpi (Fig. 1H) [22, 23]. Therefore, to investigate cellular transcriptional changes during the nascent period of pox-like lesions, another checkpoint in RNA-seq was set at 24 hpi, corresponding to the late stage of infection.
Sequencing data quality
The RNA-seq results revealed an average of 91,876,711 raw reads in Group A, 91,478,195 in Group B, and 88,921,758 in Group C (Table 2). Valid reads [24] for each group were identified, averaging 83,959,887 in Group A, 83,990,021 in Group B, and 79,158,401 in Group C, respectively. The Q30% value, which indicates a sequencing error rate of less than 0.001, was ≥ 94.97% across all nine libraries (Table 2). The average rates of mapped read in Groups A, B, and C were 93.51%, 95.20%, and 78.01%, respectively (data not shown). The observed decrease in mapping rates is likely attributable to an increase in virus transcripts. Pearson’s correlation coefficients, with a minimum value of ≥ 0.98, demonstrated high consistency in sequencing data across the nine libraries (data not shown). These results proved the reliability of sequencing data and the effect of viral infection on cellular transcription.
Table 2
Data statistics and quality control
Sample a | Raw reads | Raw data size | Valid reads | Valid data size | Valid Ratio | Q20% b | Q30% c | GC content% |
A1 | 90084168 | 13.51G | 82840092 | 12.43G | 91.96 | 99.96 | 95.77 | 49.5 |
A2 | 89830190 | 13.47G | 81903184 | 12.29G | 91.18 | 99.96 | 95.92 | 50.5 |
A3 | 95715776 | 14.36G | 87136386 | 13.07G | 91.04 | 99.96 | 96.03 | 50.5 |
B1 | 93880978 | 14.08G | 86085200 | 12.91G | 91.70 | 99.96 | 95.94 | 51.5 |
B2 | 91831440 | 13.77G | 84712310 | 12.71G | 92.25 | 99.97 | 95.90 | 51.5 |
B3 | 88722168 | 13.31G | 81172554 | 12.18G | 91.49 | 99.96 | 95.72 | 51.5 |
C1 C2 C3 | 91556164 92074642 83134468 | 13.73G 13.81G 12.47G | 82680404 82504296 72290504 | 12.40G 12.38G 10.84G | 90.31 89.61 86.96 | 99.96 99.96 99.96 | 94.97 95.08 95.22 | 51.0 51.5 49.5 |
a: the 4-hour non-infected control Group B (B1, B2, B3), the 4-hour post-infection (hpi) LSDV-infected experimental Group A (A1, A2, A3), and the 24 hpi LSDV-infected experimental Group C (C1, C2, C3). |
b: Q20 represents reads with a 1% error rate. |
c: Q30 represents reads with a 1‰ error rate. |
Sequences characterization
From 4 hpi to 24 hpi, there was an increase in the rate of exon transcript counts, while the rate of intron transcript sequences experienced a decline (Fig. 2A) [14]. Moreover, lincRNA (unknown and intergenic transcript lncRNA) showed a dominant rate in nine libraries (data not shown). Sense overlap lncRNA (generic exonic overlap with a reference transcript lncRNA) and sense intron lncRNA (a transfrag falling entirely within a reference intron lncRNA) showed an upward trend. In contrast, antisense lncRNA (exonic overlap with reference on the opposite strand) presented a downward trend. These findings suggest that viral infections impact the transcription processes of host cells.
The majority of the transcripts were predominantly longer than 1000 nucleotides (nt) (Fig. 2B). Additionally, the quantity of exons found within lncRNAs was fewer than those in mRNA (Fig. 2C). Similarly, both the expression levels and the number of lncRNA items were lower compared to those of mRNAs (Fig. 2D).
Differentially expressed analysis
The differentially expressed lncRNA (delncRNA) and mRNA (demRNA) play a crucial role in mediating host-virus interactions, influencing the regulation of various signaling pathways. For lncRNA, 108 lncRNAs were differentially expressed at 4 hpi, with 60 lncRNAs upregulated and 48 lncRNAs downregulated (Fig. 3A). 1101 lncRNAs were differentially expressed at 24 hpi, with 528 lncRNAs upregulated and 573 lncRNAs downregulated (Fig. 3A). For mRNA, the number of demRNAs was 798 at 4 hpi, with 432 mRNAs upregulated and 366 mRNAs downregulated (Fig. 3B). 5295 mRNAs were differentially expressed at 24 hpi, with 2260 mRNAs upregulated and 3035 mRNAs downregulated (Fig. 3B). The volcano plots in Fig. 3C to 3F highlighted the notable increase in both the quantity and significance of the delncRNAs and demRNAs from 4 hpi to 24 hpi.
DelncRNA/mRNA cis-regulation networks and PPI networks
The two principal mechanisms of lncRNA-mediated transcriptional regulation are cis-regulation, which affects neighboring genes, and trans-regulation, which influences genes across chromosomes. Because the substantial data volume associated with trans-regulation predictions was usually too large to validate efficiently, we only predicted cis-regulation in this study. A total of 2253 lncRNAs and 5275 mRNAs were predicted to have cis-regulation relationships.
During the early stage of infection (4 hpi), five pairs exhibited a co-expression relationship, including CD46/ENSBTAT00000079804, CHCHD6/MSTRG.9077.4, MLLT11/MSTRG.11757.3, DYNC2I1/MSTRG.12938.2, and ENSBTAG00000048990/MSTRG.6257.3 (Fig. 4A). Cis-regulation network was constructed based on delncRNA/mRNA pair due to the limited number of co-expression relationship pairs (Fig. 4A). MSTRG.13499.1 and MSTRG.13499.3, along with MSTRG.7172.2 and MSTRG.7172.3, are interconnected due to their close proximity within the genome. The protein-protein interaction (PPI) network, constructed based on cis-regulation targets of delncRNAs, highlighted that LENG1 and PRPF31 served as central hubs in orchestrating this network (Fig. 4B). Subsequently, a comprehensive PPI network was established using all the demRNAs identified at 4 hpi, with the six most interactive clusters being further expanded (Fig. 5A). The impacted biological processes associated with each cluster were identified as follows: cytokine-induced inflammation (Fig. 5B), ECM-receptor interaction (Fig. 5C), citrate cycle (Fig. 5D), apoptosis induced by quercetin and NF-κB/AP-1 (Fig. 5E), p130Cas linkage to MAPK signaling for integrins (Fig. 5F), and postmitotic nuclear pore complex reformation via SUMOylation (Fig. 5G).
During the late stage of infection (24 hpi), co-expression relationships (n = 615) were identified between 380 delncRNAs and 475 demRNAs. Instead of presenting the delncRNA/demRNA cis-regulation network in a single figure, it was divided into Fig. 6A and Fig. 6B for clearer visualization. 241 delncRNAs/demRNA pairs with interconnection relationships between different pairs were shown in Fig. 6A. The top five genes with the highest number of interconnections, listed in descending order, were WDR35 (n = 6), JTB (n = 5), ENSBTAG00000024487 (n = 5), PMEPA1 (n = 5), and RAB42 (n = 4). For instance, WDR35, which has the most interconnections (n = 6), is associated with delncRNAs including MSTRG.2530.1, MSTRG.2530.2, MSTRG.2530.4, MSTRG.2530.7, MSTRG.2530.8, and MSTRG.2530.9, all of which are closely located on the genome. In contrast, 374 delncRNA/demRNA pairs that lack interconnections between different pairs were shown in Fig. 6B. Furthermore, a PPI network based on co-expression relationship pairs was constructed (Fig. 7A). The top six interactive clusters, displayed in Figs. 7B to 7G, include functions related to nucleotide excision repair in xeroderma pigmentosum (Fig. 7B), mRNA processing (Fig. 7C), ubiquitin protein ligase binding (Fig. 7D), mitochondrial translation termination (Fig. 7E), assembly of mitochondrial complex III (Fig. 7F), and the Ubl conjugation pathway (Fig. 7G). Unfortunately, the PPI network based on 5295 demRNAs at 24 hpi was too complex for detailed presentation in an image, and therefore, its network was not constructed in this study.
Functional enrichment analyses
At 4 hpi, delncRNAs and demRNAs enrichment analyses were performed, respectively. For delncRNA, GO functional enrichment analysis highlighted significant enrichment in processes related to glycometabolism, growth, apoptosis, morphology, and host cell immunity (Fig. 8A). Similarly, KEGG pathway analysis revealed enrichment in glycometabolism processes, including starch and sucrose metabolism, cAMP signaling pathway, glucagon signaling pathway, and insulin resistance signaling (Fig. 8B). These findings revealed that glycometabolism control could be an important interaction point during the early infection stage under lncRNA regulation. For demRNA, GO functional enrichment analysis identified processes associated with the regulation of transcription, cellular responses to lipopeptides of pathogens, cell differentiation, and apoptosis as primarily enriched (Fig. 8C). The KEGG pathway analysis suggested that immunity-related signaling pathways were mainly enriched, including inflammatory cytokines, cell immunity, and Toll-like receptor-mediated innate immunity (Fig. 8D). These insights revealed that immunity response might be central to mRNA regulation during the early infection stage (4 hpi). At this time, cells actively engage in glycometabolism and lipopeptide-induced immune response against the invading poxvirus.
At 24 hpi, for delncRNAs, the GO functional enrichment analysis suggested that lncRNA might play important roles in various biological processes, including mechanical stimulus, RNA polymerase II preinitiation complex assembly, Golgi cisterna membrane, positive regulation of mitotic cell cycle, neuron maturation, and iron ion homeostasis (Fig. 9A). The KEGG pathway analysis further indicated an impact on cell metabolism, particularly in glycometabolism (Fig. 9B). Subsequently, to explore pivotal lncRNA/mRNA pairs, co-expression KEGG networks were constructed, focusing on life activities during virus infection, such as cell proliferation, death, and immunity (Fig. S1). Notably, cis-regulations involving MSTRG.10125.1/MAPK3 and MSTRG.10127.1/MAPK3 were central to these networks, linking to most pathways. For demRNAs, GO functional enrichment analysis indicated that oxidation-reduction process, antigen processing and presentation, and actin-mediated molecular transportation were affected (Fig. 9C). Particularly noteworthy is the enrichment of actin-mediated molecular transportation, which aligns with the observed phenomenon in reality [5]. The Arp2/3-mediated actin-based motility was abused by several pathogens, including poxvirus [25]. Meanwhile, KEGG pathways analysis showed that pathways related to disease were predominantly enriched during the early infection stage (Fig. 9D). At this time, cells were distressed and had been transformed into LSDV replication factories.
Validation for RNA-seq and identification of antiviral lncRNAs
To confirm the accuracy of results tested by RNA-seq, the total RNA of all samples was further used to detect the expression level of transcripts by RT-qPCR. The relative expression of five lncRNAs (Fig. 10A) and five mRNAs (Fig. 10C) tested with RT-qPCR corresponded to those in RNA-seq. The linear regression analyses indicated favorable correlation between the outcomes obtained from RNA-seq and RT-qPCR (Fig. 10B and 10D).
Two delncRNAs with relatively high expression levels at 4 hpi were selected to identify whether they were host antiviral effectors or not. Overexpression of these two recombinant plasmids was validated via RT-qPCR (data not shown). In vitro experiments demonstrated that overexpression of the lncRNAs AV5191 (AV: antiviral) or AV1557 can inhibit virus replication at both the mRNA level (Fig. 11A) and titer level (Fig. 11B). Based on the prediction of cis-regulation, CCDC62, ENSBTAG00000050187, HCAR1, and DENR are posited to serve as intermediaries in the lncRNA AV5191-mediated suppression of virus replication. Similarly, POLE3, ALAD, HDHD3, C8H9orf43, BSPRY, RGS3, WDR31, and ENSBTAG00000055315 are proposed as potential downstream targets of lncRNA AV15574. These findings underscore the significance of host lncRNA differential expressions during the initial phase of infection in influencing antiviral immunity or virus replication.