A comprehensive analysis of avian lymphoid leukosis-like lymphoma transcriptomes including identification of LncRNAs and the expression profiles

Abstract

Avian lymphoid leukosis-like (LL-like) lymphoma has been observed in some experimental and commercial lines of chickens that are free of exogenous avian leukosis virus. Reported cases of avian lymphoid leukosis-like lymphoma incidences in the susceptible chickens are relatively low, but the apathogenic subgroup E avian leukosis virus (ALV-E) and the Marek's disease vaccine, SB-1, significantly escalate the disease incidence in the susceptible chickens. However, the underlying mechanism of tumorigenesis is poorly understood. In this study, we bioinformatically analyzed the deep RNA sequences of 6 lymphoid leukosis-like lymphoma samples, collected from susceptible chickens post both ALV-E and SB-1 inoculation, and identified a total of 1,692 novel long non-coding RNAs (lncRNAs). Thirty-nine of those novel lncRNAs were detected with altered expression in the LL-like tumors. In addition, 13 lncRNAs whose neighboring genes also showed differentially expression and 2 conserved novel lncRNAs, XLOC_001407 and XLOC_022595, may have previously un-appreciated roles in tumor development in human. Furthermore, 14 lncRNAs, especially XLOC_004542, exhibited strong potential as competing endogenous RNAs via sponging miRNAs. The analysis also showed that ALV subgroup E viral gene Gag/Gag-pol and the MD vaccine SB-1 viral gene R-LORF1 and ORF413 were particularly detectable in the LL-like tumor samples. In addition, we discovered 982 novel lncRNAs that were absent in the current annotation of chicken genome and 39 of them were aberrantly expressed in the tumors. This is the first time that lncRNA signature is identified in avian lymphoid leukosis-like lymphoma and suggests the epigenetic factor, lncRNA, is involved with the avian lymphoid leukosis-like lymphoma formation and development in susceptible chickens. Further studies to elucidate the genetic and epigenetic mechanisms underlying the avian lymphoid leukosis-like lymphoma is indeed warranted.

Fig 1. RNA-seq reads mapped to genomes of AF227 and SB-1 viruses.

(A) Clean RNA-seq reads that failed mapping to the chicken reference genome were further mapped to the genomes of AF227 (NCBI # MF817820) and/or SB-1 (NCBI # HQ840738) viruses. The total number of mapped reads were counted for each virus in each sample. (B) Distribution of mapped reads of the LL-like tumor samples across the genome of AF227 and (C) SB-1 virus. The Y axis indicates the average depth of each location of the six tumor samples.

Fig 2. Cluster analysis of the samples.

(A) A total of 12,457 genes including 10,601 protein-coding genes, 874 known lncRNAs and 982 novel lncRNAs with raw reads counts >10 in at least 4 samples from at least one group were used for clustering and differential expression analysis. (B) Principal component analysis (PCA) of all samples based on the 12,457 expressed genes. Samples from the same group are indicated by a same color. (C) Hierarchical clustering results of all samples based on the 12,457 expressed genes. Spearman’s correlation was used as distance metric with the average linkage algorithm.

Fig 3. Differentially expressed protein-coding genes and function enrichment analysis.

(A) Volcano plot showing the numbers of significantly up- and down-regulated protein-coding genes in tumor samples as compared to normal bursal and (B) splenic B cell controls. (C) Venn diagram showing the number of common significantly up- and down-regulated protein-coding genes identified by both comparison pairs of tumor samples vs normal bursal controls and tumor samples vs splenic B cell controls. (D) Significantly enriched KEGG pathways and (E) GO terms of biological process for significant up- and down-regulated PCGs. The top 10 most significantly enriched GO terms of biological processes for down-regulated protein-coding genes were shown.

Fig 4. Analysis of differentially expressed lncRNAs.

(A) Volcano plots showing the numbers of significantly up- and down-regulated lncRNAs in the LL-like tumor samples as compared to normal bursal and (B) splenic B cell controls. (C) Number of significantly up- and down-regulated commonly known/novel lncRNAs identified by both comparisons between the LL-like tumor samples and normal bursal controls, and LL-like tumor samples and splenic B cell controls. (D) Correlation analysis between the expression of significantly differentially expressed (DE) or non-DE lncRNAs and their protein-coding gene partners. (E) Significantly dysregulated lncRNA with their protein-coding gene partners exhibited significantly differential expression in the LL-like tumors as compared to the normal controls as well. (F) Synteny and sequence conservation of the novel lncRNA XLOC_0022595 and (G) XLOC_001407. The PhastCon plot is relative to the loci of human genome.

Fig 5. Construction of lncRNA-miRNA-mRNA regulatory network.

(A) Predicted lncRNA-miRNA-mRNA regulatory network. LncRNAs with at least 10 predicted binding sites for at least 1 miRNA associated with ALV invasion were included. (B) Top 10 most significantly enriched KEGG pathways for the 178 differentially expressed protein coding genes identified in the lymphoma samples and potentially regulated by lncRNA-miRNA axis. (C) Sankey diagram showing the regulatory network mediated by the novel lncRNA XLOC_004542 that harbors the largest number of predicted miRNA binding sites. The width of lines linking XLOC_004542 and miRNAs is proportional to the putative number of miRNA binding sites, which is labeled on the lines.

Fig 6. Protein coding genes and lncRNAs that were significantly correlated with virus (A) SB-1 or (B) AF227 expression.

The number of RNA-seq reads mapped to each virus genome in each of the tumor samples were normalized to CPM and Pearson correlation between SB-1 and AF227 abundance and all the significantly differentially expressed genes was calculated in 6 tumor samples. p <0.05 was used as the threshold to determine statistical significance. The lncRNAs are highlighted in red.

Fig 7. Validation of RNA-seq results by ddPCR.

(A) Correlation analysis indicating that the expression levels of the 7 selected lncRNA/protein-coding gene pairs identified by RNA-seq and ddPCR were significantly positively correlated with a R² = 0.7475 and p value < 0.0001. (B) The abundance of the selected lncRNAs and their protein-coding gene partners in individual samples from each of the bursa (n = 3), B cell (n = 3) and the LL-like tumor (n = 3) groups were confirmed by ddPCR (*p < 0.05 in an un-paired t test).