RNA-Seq and miRNA-Seq data from Epstein-Barr virus-infected tree shrews reveal a ceRNA network contributing to immune microenvironment regulation

ABSTRACT

Epstein-Barr virus (EBV) infection in humans is ubiquitous and associated with various diseases. Remodeling of the immune microenvironment is the primary cause of EBV infection and pathogenesis; however, the underlying mechanism has not been fully elucidated. In this study, we used whole-transcriptome RNA-Seq to detect mRNAs, long non-coding RNAs (lncRNA), and microRNA (miRNA) profiles in the control group, 3 days, and 28 days after EBV infection, based on the tree shrew model that we reported previously. First, we estimated the proportion of 22 cell types in each sample using CIBERSORT software and identified 18 high-confidence DElncRNAs related to immune microenvironment regulation after EBV infection. Functional enrichment analysis of these differentially expressed lncRNAs primarily focused on the autophagy, endocytosis, and ferroptosis signalling pathways. Moreover, EBV infection affects miRNA expression patterns, and many miRNAs are silenced. Finally, three competing endogenous RNA regulatory networks were built using lncRNAs that significantly correlated with immune cell types, miRNAs that responded to EBV infection, and potentially targeted the mRNA of the miRNAs. Among them, MRPL42-AS-5 might act as an hsa-miR-296-5p “sponge” and compete with target mRNAs, thus increasing mRNA expression level, which could induce immune cell infiltration through the cellular senescence signalling pathway against EBV infection. Overall, we conducted a complete transcriptomic analysis of EBV infection in vivo for the first time and provided a novel perspective for further investigation of EBV-host interactions.

KEYWORDS:

• Epstein-Barr virus
• tree shrew
• long non-coding RNAs
• microRNA
• immune microenvironment

Introduction

The Epstein-Barr virus (EBV) is one of the most ubiquitous γ-herpes viruses that establish lifelong persistent latent infection in over 90% of adults worldwide [1]. EBV is associated with several diseases, including acute or chronic inflammatory diseases and various human malignancies of lymphoid or epithelial origin [2]. Due to its co-evolution with its host, EBV has developed sophisticated strategies to evade human immune surveillance and ensure its replication and survival [3]. These processes are often accompanied by the remodelling of the microenvironment, making it more conducive to cellular transformation and tumour progression [4]. Non-coding RNA (ncRNAs) have received increasing attention due to their diverse functions and participation in many biological pathways. Herpesviruses have evolved strategies to evade host immune responses, including using microRNAs (miRNAs) to inhibit and hijack long non-coding RNA (lncRNAs) in vitro experiments. For example, in EBV-positive gastric carcinoma and nasopharyngeal carcinoma cell lines, EBV-encoded miR-BART6-3p inhibits cancer cell metastasis and invasion by targeting lncRNA LOC553103 [5]. LncRNAs can be implicated in transcription regulatory processes through the competing endogenous RNA (ceRNA) mechanism, which indirectly inhibits miRNAs by regulating their target mRNA [6]. However, no systematic study has been done on the interaction between lncRNAs and miRNAs in remodelling the microenvironment after EBV infection, owing to the lack of suitable animal models. We previously reported that the tree shrew (Tupaia belangeri chinensis) is a suitable animal model for evaluating the mechanisms of EBV infection [7].

In this study, we acquired RNA-Seq and miRNA-Seq data from the peripheral blood of the tree shrew model at pre-infection, acute infection (3 days post EBV infection), and chronic infection (28 days post EBV infection). First, we estimated the proportion of different cell types in each sample from RNA-Seq data to show the altered landscape of immune cells during different stages of infection. Furthermore, we extracted lncRNA expression profiles from all RNA-Seq data and filtered out lncRNAs with high expression abundance and significant correlations with changes in the immunological microenvironment. Consensus clustering was used to identify the expression patterns of infection-related miRNAs using the K-means method. Finally, based on these regulatory relationships, we built a ceRNA regulatory network using immune microenvironment regulation-related lncRNAs, miRNAs with significant changes during infection, and potential target mRNA of the miRNAs.

Results

Dynamic changes of immune microenvironment characteristics in whole blood samples among healthy controls and EBV-infected tree shrew

The experimental workflow is shown in Figure 1a. First, we constructed a tree shrew model of EBV infection by intravenous inoculation of EBV suspensions through the femoral vein. Tree shrews developed variable levels of viraemia within a week post-infection, and EBV DNA copies decreased over time, allowing the virus to enter a latent state gradually (Figure 1b). Moreover, immunohistochemistry against EBV latent membrane protein 1 (LMP-1) showed a detectable viral protein in the spleen at 28 days post-infection (dpi), and EBV-encoded small RNAs (EBERs) were detected by in situ hybridization (ISH) (Figure 1c).

Figure 1. Dynamic immune microenvironment characteristic changes in whole blood samples among healthy controls (ctrl) and Epstein-Barr virus (EBV)-infected tree shrew. (a) Flow chart of the study design. (b) The EBV DNA loads in blood at different time points are shown. The red dotted lines indicate the limit of detection of the assay. (c) Detection of EBV latent membrane protein 1 (LMP-1) and EBV-encoded small RNAs (EBERs) in tree shrew spleens using immunohistochemistry and in situ hybridization at the endpoint (28 days post-infection [dpi]) of animal experimentation. Positive immunoreactivity of LMP-1 protein is indicated by brown staining in the first row. The red color represents the positive staining of EBERs in the second row. (d) Boxplot showing the fraction of each cell type in each group. The significant difference in the immune cell fractions among the ctrl, 3 dpi, and 28 dpi samples was calculated using the Student’s t-test. *, p ≤ 0.05; **, p ≤ 0.01. (e) Fractions of immune cells rank ordered based on decreasing values of the relative frequency ratio at the 3 dpi group relative to the ctrl group (the leftmost panel), 28 dpi group relative to the ctrl group (middle picture), and 28 dpi group relative to the 3 dpi group (the rightmost panel).

Human microbial infections typically induce local and systemic immune responses in the host, and the cell components in the immune microenvironment changes during disease progression [8–10]. In this study, CIBERSORT software [11] was used to estimate relative proportions of 22 types of infiltrated immune cells present in each sample of the healthy control (Ctrl group), 3 dpi, and 28 dpi groups (fractions of different cells estimated by CIBERSORT in each sample are shown in Supplementary Figure S1). After EBV infection, dynamic changes in the microenvironment composition occurred over time (Figure 1d). Notably, the relative proportion of naive B cells gradually declined, and the differences became significant at 28 dpi (compared with the control group, p < 0.05). Similar changes were observed in the relative proportion of resting NK cells (compared to the Ctrl group, p < 0.01). The naïve CD4 T cell level markedly increased at 3 dpi (compared to the Ctrl group, p < 0.05) and then rapidly decreased at 28 dpi (compared to the 3 dpi group, p < 0.01). The relative proportion of immune cells was log-transformed to present the data clearly (Figure 1e). At 3 dpi, monocytes, naïve CD4 T cells, and M2 macrophages increased, and M0 macrophages increased substantially compared to the Ctrl group. In the 28 dpi group, increases were observed in monocytes and neutrophils, and naïve B, naïve CD4 T, naïve T, and resting NK cells decreased compared to the Ctrl group.

Identification of lncRNAs in response to EBV infection in tree shrews

To obtain a comprehensive list of lncRNAs associated with infiltrating immune cells, we first performed transcriptome assembly and identified novel lncRNAs (Supplementary Figure S2A shows the bioinformatics analysis pipeline). LncRNAs must have a read per kilobase of exon per million reads mapped (FPKM) value of ≥ 0.2 in at least one sample to be considered expressed and to reduce analytical errors. After trimming, 273, 308, and 321 lncRNAs remained in the Ctrl, 3 dpi, and 28 dpi groups, respectively, 198 of which were shared among the three groups (Figure 2a). Differential analyses of the mRNAs (DEmRNAs) and lncRNAs (DElncRNAs) were conducted using DEseq2 (details are provided in the Methods). The number of DEmRNAs or DElncRNAs was at a maximum at 3 dpi and gradually recovered on day 28 dpi (Figure 2b).

Figure 2. Identification of lncRnas in response to EBV infection in tree shrews. (A) Venn diagram showing the number of detected lncRnas. The known and novel predicted lncRnas were combined, and more than one sample with reads per kilobase of exon per million reads mapped (FPKM) ≤ 0.2 was supposed to be detected. (B) Bar plot showing the number of upregulated and downregulated annotated genes (DEG, left panel) and lncRNA (DElncRNA, right panel). DESeq2 was used for differential expression analysis with criteria of false discovery rate (FDR) ≤ 0.05 and fold change (FC) ≥ 2 or ≤ 0.5. (C) Principal component analysis (PCA) based on FPKM values of all DElncRNAs. The ellipse for each group is the confidence ellipse. (D) Compared within any two groups, DElncRNAs with an absolute FPKM value ≥ 1 of the difference between the two groups were screened out. Heatmap showing the expression profile of these lncRnas. (E) Co-expression analysis of DElncRNA from (D) with at least one immune cell with p ≤ 0.01 and Pearson’s correlation value ≥ 0.8 or ≤ −0.8. *, p ≤ 0.05; **, p ≤ 0.01. (F) Co-expression analysis was performed on the 18 DElncRNA in (E) and DEmRNAs. Cutoffs of P-value ≤ 0.01 and Pearson’s coefficient ≥ 0.6 or ≤ −0.6 were used to identify the co-expression pairs. The bar plot exhibits the most enriched KEGG pathways of co-expressed DEmRNA genes. dpi, days post-infection.

Moreover, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation and pathway enrichment analyses of the upregulated and downregulated DEmRNAs. Upregulated mRNAs were primarily enriched in pathways related to programmed cell death functions, such as autophagy, endocytosis, and ferroptosis signalling pathways. In contrast, downregulated mRNAs were predominantly enriched in the NF-kappa B and viral protein interaction with cytokine and cytokine receptor signalling pathways (Supplementary Table S2 and Supplementary Figures S2F and G). Our finding suggested that transcriptional repression of the NF-kB signalling pathway might be one of the main mechanisms of action against EBV infection in tree shrews.

Notably, within the three groups, the principal component analysis (PCA) of the DElncRNAs separated them from each other (Figure 2c). These DElncRNAs may be one of the major factors responsible for EBV infection in the tree shrew model. Among these DElncRNAs, we further screened 18 DElncRNAs (for example: ADCY5-AS-5, CCDC51-AS-1, DKC1-AS-2 and so on) with high relative expression (Figure 2d, absolute FPKM of the difference ≥ 1) and significant co-expression patterns. A significant correlation between these DElncRNAs and various immune cell types (Pearson’s correlation value ≥ 0.8 or ≤ −0.8, p < 0.01) was observed (Figure 2E). Finally, to link the lncRNA to its associated protein-coding gene, lncRNA co-expressed mRNAs were extracted (Pearson’s correlation value ≥ 0.6 or ≤ −0.6, p ≤ 0.01), and KEGG pathway analysis was performed (http://www.genome.jp/kegg). Notably, similar to our previous enrichment results, autophagy, endocytosis, and ferroptosis signalling pathways were significantly enriched (Figure 2F).

Furthermore, we used transmission electron microscopy (TEM) and Western blotting (WB) to examine the autophagy and ferroptosis signalling pathways. The cristae were well preserved, and the electron density of the mitochondrial matrix appeared regular in the Ctrl group (Figure 3A). In contrast, cells undergoing ferroptosis in the 28 dpi group had several distinguishing features, including mitochondrial atrophy, higher mitochondrial membrane density, reduced mitochondrial cristae, and outer mitochondrial membrane rupture. Double-membrane autophagosomes were observed in the 28 dpi group. WB analysis showed that the expression of the autophagy marker LC3II/I and Beclin-1 was significantly increased, and that of p62 was significantly decreased after EBV infection in the 28 dpi group (compared to the Ctrl group, p < 0.05). Simultaneously, EBV infection significantly reduced the protein expression of the ferroptosis biomarker proteins SLC7A11, NRF2, and GPX4, whereas that of FTH1 (a major intracellular iron storage protein) was significantly increased (28 dpi vs. Ctrl group, p < 0.05) (Figure 3B).

Figure 3. Autophagy and ferroptosis signalling pathways are activated in animals infected with Epstein-Barr virus (EBV). (A) Transmission electron microscopy examinations. Black arrows: normal mitochondria; green arrows: mitochondria with varying cristae densities; red arrows: mitochondria with swelling and cristolysis; yellow arrows: autophagosome. (B) Western blotting (WB) analysis of autophagy and ferroptosis signalling pathways related proteins. ACTIN served as the internal control for the WB assay. *, p ≤ 0.05; **, p ≤ 0.01.

Dynamic miRNA changes before and after EBV infection in tree shrews

The pairwise differential expression analysis was first implemented to identify miRNAs with significant differential expression (DEmiRNAs, fold change (FC) value ≥ 1.5 or ≤ −1.5, p < 0.05). These results are summarized in the Upset plot, and only the DEmiRNAs with significant variance were retained in the chosen comparison (Figure 4A). miRNAs were clustered based on K-means, resulting in the five groups in Figure 4B and Supplementary Figure S3. Specifically, the K3 group (252 miRNAs) was predominately characterized by downregulation and contained the largest number of miRNAs, whereas the K2 group (28 miRNAs) showed a constant downregulation pattern. The K4 group (30 miRNAs) showed little change at 3 dpi compared to that of the Ctrl animals but significantly increased at 28 dpi. miRNA expression in the K5 group (16 miRNAs) was upregulated at first (3 dpi) and then returned to Ctrl levels by 28 dpi. Because a sample (F12) with strong heterogeneity was present in the K1 group (17 miRNAs), this group was excluded from subsequent analyses (Figure 4C).

Figure 4. Dynamic miRNA changes before and after Epstein-Barr virus (EBV) infection in tree shrews. (A) Upset visualization of the DEmiRNA intersection size from different data sets (black points connected by a line). (B – C) DEmiRNAs among control (ctrl), 3 days post-infection (dpi), and 28 dpi were clustered using K-means. The heatmap in (B) presents five clusters of DE miRnas identified by K-means clustering, and the heatmap in (C) shows expression patterns according to clustered miRnas. K2, K3, K4, and K5, EBV-related clusters, were labelled red.

Construction of ceRNA network mediated by immune microenvironment-related lncRnas during EBV infection

According to the ceRNA hypothesis, lncRNAs can serve as ceRNAs to compete for miRNA response elements, thereby regulating the expression of target genes [6]. A ceRNA regulatory network was constructed using 18 lncRNAs related to immune cell changes from Figure 2E, the miRNAs regulated by EBV infection identified in Figure 4C (except for the K1 cluster), and miRNA-targeted mRNAs, to gain insight into the function of immune-related lncRNAs during EBV infection. The interaction between lncRNAs and miRNAs was identified using miRanda [12] (score >150) and Rnahybrid [13] (p < 0.05), and the results were considered as the intersection. We identified three lncRNA-miRNA pairs: MRPL42-AS-5/hsa-miR-296-5p, GATA1-AS-1/hsa-miR-423-5p, and PDZK1IP1-AS-1/hsa-miR-671-5p (Figure 5B,C).

Figure 5. Construction of competing endogenous RNA (ceRNA) network mediated by immune microenvironment-related lncRnas during Epstein-Barr virus (EBV) infection. (A) The ceRNA regulatory network was constructed using 18 lncRnas related to immune cell changes from figure 2E, the miRNA regulated by EBV infection identified in figure 3C (except the K1 cluster), and miRNA-targeted mRnas. Cytoscape was used to plot the network. (B) Boxplot showing the expression profile of three lncRnas in the ceRNA network shown in (A). (C) Boxplot showing the expression profiles of three miRnas in the ceRNA network shown in (A). (D) The boxplot shows the expression profile of target mRnas of hsa-miR-296-5p/MRPL42-AS-5. (E) Reverse transcription quantitative polymerase chain reaction (RT-qPCR) was used to validate the seven targeted genes expression of hsa-miR-296-5p/MRPL42-AS-5. (F) RNA fluorescence in situ hybridization (FISH) to detect miR-296-5p expression level in the first row in the figure; the red fluorescence indicates positive signals, and the blue fluorescence shows the nuclei. In the second row, immunofluorescence experiments were performed. Green fluorescence shows the CCND3 protein, and blue fluorescence shows the nuclei using immunofluorescence. The plots on the right side of the figure show the miR-296-5p and CCND3 scores separately, determined by counting the ratio of positive cells from five random regions in each section. *, p ≤ 0.05; **, p ≤ 0.01. Ctrl, control.

First, through co-expression analysis, the miRNA-mRNA co-expression pair with a correlation coefficient of ≤ −0.8 (negative correlation) and a p value of ≤ 0.01 was acquired to ascertain the relationship between miRNA and mRNA. On the one hand, we overlapped the targets of three miRNAs (hsa-miR-296-5p/hsa-miR-423-5p/hsa-miR-671-5p) in Targetscan [14] or mirDB database (http://mirdb.org/) with the miRNA-mRNA co-expression pair. On the other hand, miRanda [12] and Rnahybrid [13] were also used to identify potential new target relationships and overlapped with the miRNA-mRNA co-expression pair. Finally, the above results were combined to determine the potential mRNA-target relationships of the miRNAs (Supplementary Figure S4B).

LncRNA-miRNA and miRNA-mRNA pairs sharing the same miRNA were merged into a lncRNA-miRNA-mRNA interaction as a candidate ceRNA regulatory network (Figure 5A). A significant increase in MRPL42-AS-5 expression after infection was observed (Figure 5B). However, hsa-miR-296-5p, which belongs to the same ceRNA network, exhibited the opposite response (Figure 5C). Moreover, mRNA co-expressed with hsa-miR-296-5p showed similar trends to MRPL42-AS-5 changes (Figure 5D). The target genes in the ceRNA network were primarily associated with the cellular senescence signalling pathway (Supplementary Table S4).

To validate the RNA-Seq data, we performed reverse transcription quantitative polymerase chain reaction (RT-qPCR) to detect the seven target genes (Figure 5E). In addition, we performed fluorescence in situ hybridization (FISH) to detect the expression of hsa-miR-296-5p, and the expression of the target protein CCND3 was assessed by immunofluorescence (IF)(Figure 5F).

Discussion

Viral infections trigger a broad host immune response that controls viral replication [15]. However, our understanding of the immunological mechanisms involved in EBV infection is limited, particularly under in vivo conditions, because there are no suitable animal models. Previously, we found that tree shrews are susceptible to EBV infection and exhibit characteristics similar to those of humans [7,16]. Using the tree shrew model, we found dramatic changes within the blood microenvironment with the progression of EBV infection. Notably, the number of naïve CD4+ T cells significantly increased at 3 dpi and then decreased to normal levels at 28 dpi. This finding might result from the proliferation and differentiation of this subgroup of T cells due to stimulation with viral antigens at an early stage of infection. However, not all these changes reached statistical significance, possibly due to the small sample size and individual heterogeneity of the groups.

Increasing evidence suggests that ncRNA is widely involved in immune regulation during viral infection [17–19]. LncRNAs regulate diverse biological functions (e.g. cell survival and cytokine production) in innate immunity and are involved in adaptive immunity (e.g. polarization and function of CD4+ T cells) [20–22]. To analyse DElncRNAs between healthy and infectious conditions, we compared the subgroups separately to the Ctrl group and found that up- or downregulated DElncRNAs showed a much more dramatic shift at 3 dpi. In addition, the PCA showed clear group segregation at various time points post-infection. This finding implies that DElncRNAs play multiple roles during the different stages of EBV infection. Notably, more discrete 3 dpi group samples were presented in the PCA graph, indicating that higher lncRNA-regulated heterogeneity and drastic changes in transcript levels may have amplified this heterogeneity during the initial stages of infection.

Next, we screened 18 highly expressed DElncRNAs strongly associated with changes in the proportion of immune cells through a combined analysis of correlation and relative expression abundance. CD4+ T cell subgroups and M2 macrophages were positively or negatively correlated with multiple DElncRNAs, indicating the involvement of lncRNAs in immune cell mobilization. Co-expressed mRNAs were extracted, and KEGG pathway analysis was performed to probe the function of these 18 DElncRNAs further. Primarily enriched pathways were related to programmed cell death, such as ferroptosis, autophagy, mitophagy, and endocytosis. TLR/RLR and cGAS-STING activation by viral infection may help induce the autophagy pathway, which could enhance interferon production and thus restrict viral replication [23]. While many in vitro studies have been published on the induction of ferroptosis and autophagy signalling pathways by EBV, little is known about whether this occurs in vivo [24–26]. Our results showed that programmed cell death-related signalling pathways play an essential role in the innate immunity of animals by initiating the cell death process to defend against primary EBV infection. Moreover, the enrichment analysis of DEmRNAs showed similar results to the 18 DElncRNA co-expressed mRNAs, which implies a potentially crucial biological function for DElncRNAs involved in regulating programmed cell death.

The ability of EBV to manipulate host miRNAs to escape the immune system allows the infection to persist lifelong without being eliminated [27,28]. This is a dominant feature of EBV infection in the tree shrew model, in which miRNA expression is downregulated. Similar observations were made for Caprine parainfluenza virus type 3 [29] and HIV-1 infections [30]. miRNA expression is linked to type I interferon stimulation [31,32]. We speculate that the immune escape mechanism in EBV infection likely involves inhibiting the NF-kB signalling pathway by these miRNAs, which can affect the interferon level. However, further studies are required to investigate this hypothesis.

The antiviral immune response induces a complex regulatory network involving diverse non-coding and coding RNAs. LncRNAs can be implicated in transcription regulatory processes through the ceRNA mechanism, which indirectly inhibits miRNAs by regulating their target mRNA [6]. We constructed three regulatory networks to explore the ceRNA regulatory mechanisms of EBV infection. Notably, one of these networks showed typical ceRNA features with a significant increase in the MRPL42-AS-5 expression after infection, whereas hsa-miR-296-5p exhibited the opposite response, and the mRNA co-expressed with hsa-miR-296-5p showed similar trends to MRPL42-AS-5 changes. The target genes in the ceRNA network were primarily associated with the cellular senescence signalling pathway, a cellular stress response that recruits immune cells [33]. Therefore, it is interesting to speculate that MRPL42-AS-5 might act as an hsa-miR-296-5p “sponge” and compete with target mRNAs, thus increasing mRNA expression level, which could induce immune cell infiltration through the cellular senescence signalling pathway against EBV infection.

In conclusion, this study provides the first comprehensive analysis of the lncRNA and miRNA profiles of EBV infection at the acute and chronic stages, which could have functional roles in remodelling the immune microenvironment. Importantly, our report revealed ceRNA networks in EBV infection that regulate EBV-host interactions. For example, MRPL42-AS-5 might act as an hsa-miR-296-5p “sponge” and compete with target mRNAs, thus increasing mRNA expression level, which could induce immune cell infiltration through the cellular senescence signalling pathway against EBV infection. Our research provides clear indications of the roles lncRNAs and miRNAs play in the modulation of host immune defence mechanisms against EBV infection. These findings open avenues for developing novel therapeutic strategies and diagnostic tools aimed at combating EBV and its associated conditions. Future studies are recommended to focus on the characterization and functional analysis of these identified RNA molecules to confirm their roles and explore their therapeutic potential in greater depth.

Material and methods

Animal experiments and ethics statement

The experimental workflow is shown in Figure 1a. Six tree shrews (F1 generation, production approval number: SCXK (Dian) 2020–0004, 8 ± 2 months, male and female, 126.3 ± 5.2 g) were acquired from the Kunming Institute of Zoology, Chinese Academy of Sciences. The protocols used in this study were approved by the Animal Ethics Review Committee of the Guangx i Medical University (approval number: 202005015). The experimental protocols strictly followed the “Guiding Principles for the Use and Care of Experimental Animals” issued by the Ministry of Science and Technology of China. The total number of animals and their suffering were minimized according to the 3 R principle.

The animals were randomly assigned to either the experimental (n = 3) or control group (n = 3). The control group was inoculated with 200 µl 1640 Medium, whereas the experimental group was inoculated subcutaneously with 200 µl EBV suspensions (1 × 10⁸ copies/ml) via the femoral vein. Blood samples were collected at 3 dpi and 28 dpi for EBV infection, and the animals were euthanized at 28 dpi with an overdose of sodium pentobarbital. Spleen samples were collected immediately after euthanization. Subsequently, the blood samples (1.5 ml) were used for transcriptome sequencing. Residual blood samples and tissues were used for further laboratory experiments, including in situ hybridization, immunohistochemistry, and transmission electron microscopy. The specific details of each experiment are provided in Supplemental Experimental Procedures.

RNA-Seq data process

Clean data were obtained by removing reads containing adapters, poly-N, and low-quality reads from raw data. The clean data was higher than 12.92 Gb. The percentage of bases whose base recognition accuracy reaches the Q30 standard (the accuracy of base recognition ≥ 99.9%) was > 94.51%. All downstream analyses were based on high-quality clean data.

Reads alignment and differential gene expression analysis

Clean data were aligned to the tree shrew TS3.0 genome (http://www.treeshrewdb.org/download.html) using HISAT2 [34]. Uniquely mapped reads were used to calculate the read number and FPKM for each gene. The expression levels of genes were evaluated using FPKM. DEseq2 was used for differential gene expression analysis [35], which models the original reads and uses a scale factor to explain differences in library depth. Subsequently, DEseq2 estimates the gene dispersion and reduces these estimates to generate more accurate dispersion estimates to model the read count. Finally, the negative binomial distribution model was fitted using DEseq2, and the hypothesis was tested using the Wald or likelihood ratio tests. DEseq2 can be used to analyse the differential expression between two or more samples, and the results can be used to determine whether a gene is differentially expressed by FC (FC ≤ 2, or ≤ 0.5) and false discovery rate (FDR, FDR ≤ 0.05).

Cell-type quantification

The CIBERSORT algorithm (v1.03) [11] with default parameters was used to estimate immune cell fractions using the FPKM values of each expressed gene. A total of 22 immune cell phenotypes were analysed in this study, including seven T cell types (CD8 T, naïve CD4 T, memory CD4 resting T, memory CD4 activated T, T follicular helper, regulatory T cells [Tregs] and T cells gamma delta), naïve and memory B cells, plasma cells, resting and activated NK cells, monocytes, M0, M1, and M2 macrophages, resting and activated dendritic cells, resting and activated mast cells, eosinophils, and neutrophils.

LncRNA prediction and direction identification

We used a pipeline for lncRNA identification similar to that previously reported [36] to systematically analyse the lncRNA expression pattern, which was constructed using StringTie 1.3.3 software [37]. We used four software packages to predict the coding potentials of the lncRNAs: CPC [38], LGC [39], CNCI [40], and CPAT [41].

Differentially expressed miRnas

Differentially expressed miRNAs between paired groups were analysed using edgeR [42] in the R package. A significant p value and FDR were obtained for each miRNA based on the negative binomial distribution model. FCs in gene expression were estimated using the EdgeR statistical package. The criterion for DEmiR was set as |FC| > 1.5 and p ≤ 0.05. The potential miRNA targets were obtained from TargetScan (v8.0) and miRDB (http://mirdb.org/).

Co-expression analysis

We calculated Pearson’s correlation coefficients (PCCs) to explore the regulatory mode between lncRNAs/miRNAs and mRNAs and classified their relationship into three classes: positively correlated, negatively correlated, and non-correlated based on the PCCs value.

Other statistical analyses

The pheatmap package (https://cran.r-project.org/web/packages/pheatmap/index.html) in R was used to perform clustering based on Euclidean distance. The Student’s t-test was used for comparisons between the two groups.

Author contributions

W. X.: Conceptualization, Methodology, Investigation, Visualization, and Writing-Original Draft. N. S.: Resources and Investigation. A. Z. T. and C.Q.L.: Supervision, Project administration, and funding acquisition.

Acknowledgements

We thank the Experimental Animal Center of Guangxi Medical University for providing technical support for animal feeding. We thank the members of our research groups for providing technical assistance and participating in discussions.

Disclosure statement

No potential conflicts of interest were stated by the author(s).

Data availability statement

Raw data have been deposited at the SRA databank (www.ncbi.nlm.nih.gov/sra) and are available under the accession number PRJNA981676.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/21505594.2024.2306795

Additional information

Funding

This study was supported by the National Natural Science Foundation of China (Grant Nos. 32060132, 82260519, and U21A20371), the Guangxi Clinic Medicine Research Center of Nasopharyngeal Carcinoma (Grant No. GuikeAD20297078).