Identification of aberrantly expressed circRNAs in subgroup J avian leucosis virus induced tumor livers by RNA sequencing

ALV-J (subgroup J avian leucosis virus) is a kind of oncogenic exogenous retrovirus and diseases associated with ALV-J have caused severe reproduction problems in the poultry industry worldwide. However, the pathogenesis of ALV-Jinduced tumor formation is still unclear. In recent years, circRNAs are validated to be related to cancer, which prompts us to explore whether aberrant expression of circRNAs exists in ALV-J-infected tumor samples. In this study, a total of 2822 differentially expressed circRNAs were detected by RNA sequencing in the tumor livers of ALV-J-susceptible chickens (n = 3) and the normal liver samples of ALV-J-resistant chickens (n = 3), among which 2808 circRNAs in the tumor liver samples were detected as unregulated compared to those in the normal liver samples. As circRNAs are considered as miRNA sponge, we predicted top 5 miRNA targets for 24 selected upregulated circRNAs in the ALV-J-induced tumor livers. CircRNA-miRNA-mRNA regulatory axes were also constructed to provide clear interaction relationship among them. Furthermore, Gene Ontology and KEGG pathway analyses were performed for miRNA target genes, the result of which implied that metabolism plays a key role in circRNAs alterations in ALV-J-induced tumor livers. In summary, we identified that circRNAs are involved in ALV-J-induced tumor formation in chickens. However, the pathogenicity mechanism still needs to be further explored.


INTRODUCTION
Circular RNA (circRNA) is a novel kind of endogenous noncoding RNA that form circular structures without 5′-3′ polarities and polyadenylated tails [1].In 2017, two reports identified that circRNAs had the translation function [2,3].CircRNAs are produced from pre-mRNA back-splicing and according to the circRNA biogenesis, intronic circRNA, exonic circRNA, sense overlapping circRNA, antisense circRNA and intergenic circRNA are the five types that have been reported [4][5][6][7].CircRNAs are enriched, conserved and stable in mammalian cells [7,8].Recently, it was reported that circRNAs function as microRNA (miRNA) sponges that naturally sequester and competitively suppress miRNA activity [9].It was demonstrated that circRNAs are involved in the development of several types of diseases, such as atherosclerosis, myocardial Infarction and neuron injury [10][11][12].Indeed, there has been a special outpouring of concern for immense amounts of research shown that circRNAs are associated with many human cancers, such as hepatocellular carcinoma [13], esophageal squamous cell carcinoma [14], human gastric cancer [15], colorectal cancer [16].Our previous study showed that circRNA alterations are involved in resistance to avian leucosis virus subgroup-J-induced tumor formation in chickens [17], however, no reports showed that circRNAs play a key role in chicken cancer which was induced by avian leucosis virus subgroup-J (ALV-J).
The avian leucosis virus belongs to the Alpharetrovirus genus of the Retrvoviridae family and induces sarcoma and leucosis in chickens [18].Up to now, there are seven subgroups of ALV in chickens: ALV-A, -B, -C, -D, -E, -J and -K [19,20].ALV-J was first detected in meat-type chickens in 1987 as a new subgroup [21] and has spread to many countries such as America [22], UK [23], China [24] and so on.Clinical infection is associated with immune tolerance, high mortality, delayed growth, inducing different kinds of tumors [25,26].ALV-J can be transmitted vertically and horizontally, and horizontal transmission of ALV-J is more efficient than for other ALV subgroups [27].ALV-J has caused great economic losses in the poultry industry worldwide [28].Developed countries and China have made efforts to eradicate ALV-J from chicken breeding flocks in the poultry industry, unfortunately, an effective vaccine against ALV-J infection is not currently available [29,30].Breeding ALV-J-resistance chickens may be a new technique to reduce the severe economic cost which is caused by ALV-J and we have successfully cultivated ALV-J-resistance chickens [31].
Recently, we have identified that circRNAs alterations are involved in resistance to ALV-J-induced tumor formation in chickens [17], however, no reports showed that circRNAs aberrant expression also plays a key role in ALV-J-induced tumor formation.High-throughput RNA sequencing (RNA-seq) has been widely used to infer host gene expression in response to virus infection [32].Herein, we confirm that enriched and differentially expressed circRNAs by RNA-seq in tumor liver samples of ALV-J-susceptible chickens.Further biological function prediction analysis was also conducted to develop a better understanding of the circRNAs function in ALV-J-induced tumor chickens.

Identification of differentially upregulated circRNAs in tumor liver samples of ALV-Jsusceptible chickens
CircRNAs were sequenced in three tumor liver samples of ALV-J-susceptible chickens and three normal liver samples of ALV-J-resistance chickens.A total of 13419 circRNAs were detected after RNA-seq.Among those 13419 circRNAs, 2822 differentially expressed circRNAs (fold change ≥ 2, P ≤ 0.05, FDR < 0.05) were detected in both tumor liver samples and normal liver samples by Volcano Plot, including 2808 differentially upregulated circRNAs that were identified in tumor liver samples (Figure 1A).Scatter plot analysis made the variation of circRNAs expressions of two groups be visualized in these plots (Figure 1B).Hierarchical clustering analysis of randomly selected 24 differentially expressed circRNAs clearly screened them into two groups (Figure 2).After comparison with the genome sequence of chicken, total five types of circRNAs including exonic, intronic, sense overlapping, antisense and intergenic circRNAs were found, comprising 40.62%, 6.86%, 21.36%, 13.62% and 17.53%, respectively (Figure 3).

MiRNA response elements analysis of differentially upregulated circRNAs in tumor liver samples
For better understanding the detailed original information of circRNAs, 24 upregulated circRNAs in tumor liver samples of ALV-J-susceptible chickens were selected and ranked by fold change (FC), the 24 circRNAs were annotated in the gallus genome and the names of the best-predicted transcript for each of this 24 circRNAs were shown in Table 1.

The homology of chicken circRNAs
Comparing chicken circRNAs with those of mouse and humans, homology analysis results showed that 8564 (63.82%) chicken circRNAs were homologous with humans, 8278 (61.69%) were homologous with the mouse.Importantly, a total of 8024 chicken circRNAs (59.80%) showed homology with mouse and humans (Figure 5).

MiRNA binding sites prediction and circRNAtargeted miRNA-gene network identification
Because that circRNAs as an active miRNA sponge regulator, we then listed top 5 possible miRNA targets for each circRNA which was up-regulated in tumor liver samples of ALV-J-susceptible chickens in Table 2. CircRNA-targeted miRNA-gene network was constructed using Targetscan and miRanda.CircRNAs contain many miRNA response elements (MREs) [33], we then selected the top 5 predicted miRNA targets of 3 randomly differentially upregulated circRNAs (chr5: 51521674-51524828, chr3: 3137166-3138170 +, chr4: 59505289-59510641-) in tumor liver samples, meanwhile, the top 10 target genes for each miRNA were collected.Finally, the integrated circRNA-targeted miRNA-gene network was determined (Figure 6).

Functional annotations for miRNA target genes
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of the predicted miRNA target genes which were predicted to be targeted by the 24 differentially upregulated circRNAs in tumor liver samples of ALV-J-susceptible chickens showed that those miRNA target genes were involved in cysteine and methionine metabolism pathway (Figure 7).Furthermore, Gene Ontology (GO) annotation of biological process (BP) of those miRNA target genes participated in mRNA metabolic process, cell development, mRNA catabolic process, phospholipid metabolic process, regulation of cell projection organization, regulation of cellular component size, regulation of cellular component organization, regulation of transmembrane transporter activity, regulation of cell projection assembly and mRNA processing (Figure 8A).GO annotation of cellular component (CC) of those miRNA target genes participated in cell junction and collagen trimer (Figure 8B), while for GO annotation of molecular function (MF) of those miRNA target genes participated in nucleoside phosphate binding, nucleotide binding, small molecule binding, DNA helicase activity, adenyl nucleotide binding, adenyl ribonucleotide binding, binding, cytoskeletal protein binding, ATP binding and RNA binding (Figure 8C).

DISCUSSION
CircRNA expression has the cell-type-specific, tissue-specific or developmental-stage specific manner [34,35].To date, several studies have proved that circRNAs may regulate gene expression by manipulating miRNAs through its miRNA sponges function [36], and circRNAs also compete with mRNAs for miRNA binding in the cytoplasm and interfere in gene regulation [8].MiRNAs control a large set of biological processes through direct interaction with their target mRNA, so circRNA sponge activity will also affect these pathways [37,38].Importantly, circRNAs have become the new markers for diseases including cancers [10,39].
RNA-seq was used for discovering new transcript and genes expression levels including circRNAs expression in many cancer research [40,41,42].ALV-J mainly induces myelocytomatosis in chickens and can cause other reproduction problems in the poultry industry worldwide [43].In the present study, we want to explore whether circRNAs are involved in ALV-J-induced tumor formation or not.Interestingly, totally 13419 circRNAs were detected using RNA-seq in tumor liver samples of ALV-J-susceptible chickens and normal liver samples of ALV-J-resistance chickens, totally 2822 differentially expressed circRNAs (fold change ≥ 2, P ≤ 0.05, FDR < 0.05) were detected in both of tumor liver samples and normal liver samples, while 2808 differentially upregulated circRNAs were identified in tumor liver samples, the different expression patterns implied circRNAs may play a key role in tumor formation which is induced by ALV-J.QRT-PCR was a necessary tool to valid the RNA-seq data in many studies [44,45].Here, in our study, we confirmed the differentially expressed circRNAs from RNA-seq data by qRT-PCR.As expected, qRT-PCR results of randomly selected 4 differentially expressed circRNAs had the same expression level trend compared with that from RNA-seq data, then we were sure that the RNA-seq results were reliable.
For better understanding the function of miRNA target genes, the additional functional analysis was conducted using the GO annotation and KEGG pathway.Statistical analysis of miRNA target genes was used to determine which GO terms or KEGG pathways were significantly enriched in the tumor livers samples as compared to the normal livers samples.These enriched categories along with the number of miRNA target genes identified for biological process, cellular component, molecular function three aspects.Interestingly, GO term of biological process and KEGG pathway analysis both showed that miRNA target genes are involved in the metabolic process, that implied metabolism make an important function in tumor liver formation induced by ALV-J.We suspected that maybe those up-regulated circRNAs changed the status and then caused the physical disorder and finally result in tumor formation.We also conducted homology analysis of circRNAs among chickens, mouse and humans to explore whether differentially expressed circRNAs of tumorigenesis in chickens may also involve in tumorigenesis in humans.It is worth mentioning that homology analysis results showed that 63.82% and 61.69% of chicken circRNAs were homologous with humans and mouse respectively, while a total of 59.80% chicken circRNAs showed homology with mouse and humans, that implies those differentially expressed circRNAs of tumor livers in chickens may also involve in tumor formation in humans.
In recent years, circRNA-miRNA-gene regulatory axes were confirmed that involved in pathogenic processes or cancer [46,47].In order to show a more intuitionistic result, three up-regulated circRNAs (chr5: 51521674-51524828-, chr3: 3137166-3138170+, chr4: 59505289-59510641-) in tumor liver samples of ALV-J-susceptible chickens were selected and five target miRNAs and 10 miRNA target genes were showed in the circRNA-miRNA-gene network which provides a foundation for further pathogenic mechanism research.In conclusion, our study for the first time reveals that circRNAs are involved in ALV-J-induced tumor formation in chickens.We totally identified 2808 differentially expressed upregulated circRNAs in tumor livers of ALV-J-susceptible chickens.Based on miRNA sponge function, top 5 miRNA targets for 24 selected upregulated circRNAs in ALV-J-induced tumor livers were listed according to the match score.We also constructed the circRNA-miRNAgene network, furthermore, GO terms and KEGG pathway analyses to point out that metabolism plays a key role in circRNAs alterations of ALV-J-induced tumor livers.Homology analysis results showed that totally 59.80% chicken circRNAs showed homology with mouse and humans, that opens a scenario for the investigation of the role of those differentially expressed chicken circRNA in human tumors.However, more underlying mechanism research still needs to be done to reveal which circRNAs alteration are directly involved in ALV-J-induced tumor formation in chickens.

Ethics statement
All experiments were carried out in strict accordance with the recommendations of the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health.The use of animals in this study was approved by South China Agricultural University Committee of Animal Experiments (approval ID: 201004152).

RNA-seq specimens
Tumor livers of ALV-J-susceptible chickens (the model of ALV-J-induced tumor chickens) and normal livers of ALV-J-resistance chickens (the model of ALV-J-infected negative chickens) in the F3 generation at 20 weeks after challenge with ALV-J NX0101 strain which induced myeloid tumors were collected for RNA-seq [31].In each group, three livers were collected as three duplicates.

RNA library synthesis and RNA-seq
Total RNAs were extracted from tumor livers and normal livers using TRizol reagent (Invitrogen, Karlsruhe, Germany) according to the manufacturer's instructions.The purity and concentration of the RNA samples were determined by assessing the OD260/280 using a NanoDrop ND-2000 instrument (Thermo, Waltham, MA, USA).To enrich circRNAs, RNase R (Epicentre, Madison, WI) was used to digest the extracted total RNAs to remove linear RNAs [48].RNA libraries were constructed from the treated RNAs using the TruSeq Stranded Total RNA Library Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer's instructions.Libraries were controlled for quality and quantity according to the BioAnalyzer 2100 system (Agilent Technologies, Inc., Santa Clara, CA, USA) and then was subjected to paired-end sequencing for 150 cycles by HiSeq 4000 Sequencer according to the manufacturer's instructions (Illumina, SanDiego, CA).The paired-end reads were harvested and subject to quality control using Q30 [49].After 3′ adaptor-trimming and the removal of low quality reads using Cutadapt software [50], the high quality trimmed reads were used for analysis of circRNAs.The differential expressions of circRNAs between the two groups were detected through foldchange filtering and Student's t-test.CircRNAs exhibiting fold changes ≥ 2.0 and P-values ≤ 0.05 were selected as significantly differentially expressed circRNAs.

Quantitative real-time PCR
QRT-PCR was achieved using PrimeScript RT Reagent Kit (Perfect Real Time; TaKaRa, Osaka, Japan) on a Mx3005P real-time PCR System (Stratagene, La Jolla, CA) with SYBR Green qPCR SMix (ROX; Roche) following the manufacturer's instructions.Divergent primers, rather  Biotech (Shanghai, China) and were listed in Table 3.The 2 − ΔΔCt method was used to analyze the qRT-PCR results [46].

Bioinformatics analysis
CircRNA and miRNA interactions were predicted using customized Arraystar miRNA target prediction software based on TargetScan and miRanda [51,52]; the top five putative target miRNAs were identified for upregulated circRNAs in tumor liver samples of ALV-Jsusceptible chickens.The putative target genes of these miRNAs were identified using Targetscan [51].Cytoscape software was used to construct the circRNA-miRNA-gene networks [53].Homology analysis of chicken circRNA was conducted using blastn (v2.2.27) from the alignment between human and mouse [54].Gene ontology and KEGG pathway analysis were performed for the differentially expressed circRNA-associated genes [55,56].

Statistical analysis
Data were processed using GraphPad Prism (version 5.0) and are expressed as the mean ± SE.The Student's t-test was used to assess differences among groups, where P ≤ 0.05 was considered to show significant differences between the groups.

Figure 1 :
Figure 1: Bioinformatics analysis of differentially expressed circRNAs in tumor livers of ALV-J-susceptible chickens and normal livers of ALV-J-resistance chickens.(A) Volcano plot of differentially expressed circRNAs.(B) Scatter plot of differentially expressed circRNAs.

Figure 2 :
Figure 2: Hierarchical cluster analysis of circRNAs.circRNA expression profile of six RNA-seq samples was showed in red (upregulated expression) and green (downregulated expression) scale.The dendrogram shows the relationships among the expression levels of samples.

Figure 3 :
Figure 3:The percentage of five different types of circRNAs in chickens.

Figure 4 :
Figure 4: Validation of the expression level of four selected differentially expressed circRNAs by qRT-PCR.(A) Relative circRNA levels as indicated by RNA-seq; (B) Fold change for circRNA levels as indicated by qRT-PCR.Bars are mean ± SD (n = 3); ** P < 0.05.

Figure 5 :
Figure 5: Chicken circRNAs shared homology with mouse and humans.

Figure 6 :
Figure 6: circRNA-miRNA-gene network for three selected upregulated circRNAs in tumor livers of ALV-J-susceptible chickens.Yellow represents circRNAs, green represents miRNAs and pink represents gene.

Figure 7 :Figure 8 :
Figure 7: KEGG pathway analysis for the miRNA target genes of 24 upregulated circRNAs in tumor livers of ALV-Jsusceptible chickens.The yellow boxes indicate the genes targeted by differentially expressed circRNAs.

Table 3 : The primers were used for qRT-PCR in this study
R: AGGCATACAGGGACAGCACA www.impactjournals.com/oncotarget