- 1Heyuan Branch, Guangdong Provincial Laboratory of Lingnan Modern Agricultural Science and Technology & Guangdong Provincial Key Lab of Agro-Animal Genomics and Molecular Breeding & Key Laboratory of Chicken Genetics, Breeding and Reproduction, Ministry of Agriculture, College of Animal Science, South China Agricultural University, Guangzhou, China
- 2South China Collaborative Innovation Center for Poultry Disease Control and Product Safety, Guangzhou, China
- 3Guangdong Engineering Research Center for Vector Vaccine of Animal Virus, Guangzhou, China
- 4United States Department of Agriculture (USDA), Agriculture Research Service, Avian Disease and Oncology Laboratory, East Lansing, MI, United States
Avian Leukosis Virus Subgroup J (ALV-J) is a tumorigenic virus with high morbidity and rapid transmission. N6-methyladenosine (m6A) is a common epigenetic modification that may be closely related to the pathogenicity of ALV-J. Currently, there are no reports on whether m6A modification is related to ALV-J induced tumor formation. In this study, we used methylated RNA immunoprecipitation sequencing (MeRIP-seq) and RNA sequencing (RNA-seq) to examine the differences in m6A methylation and gene expression in normal livers and ALV-J-induced tumor livers systematically, with functional enrichment and co-expression analysis. The results identified 6,541 m6A methylated peaks, mainly enriched in CDS, and more than 83% of the transcripts contained 1-2 m6A peaks. For RNA-seq, 1,896 and 1,757 differentially expressed mRNAs and lncRNAs were identified, respectively. Gene enrichment analysis indicated that they may be involved in biological processes and pathways such as immunology-related and apoptosis. Moreover, we identified 17 lncRNAs, commonly existing in differently expressed methylome and transcriptome. Through co-expression analysis, 126 differentially expressed lncRNAs, and 18 potentially m6A-related methyltransferases were finally identified and connected, suggesting that m6A modifications might affect gene expression of lncRNAs and play a role in ALV-J induced tumor formation. This study provides the first comprehensive description of the m6A expression profile in tumor livers induced by ALV-J infection in chickens, which provides a basis for studying the role of m6A modification in ALV-J induced tumorigenesis. This study provides clues for studying the epigenetic etiology and pathogenesis of ALV-J.
Introduction
Avian Leukosis Virus (ALV) belongs to the Alpharetrovirus genus of the Retrvoviridae family that causes a variety of neoplastic diseases in chickens (1). There are seven subgroups of ALV in chickens, including subgroups A, B, C, D, E, J and K (2). ALV-J is an oncogenic retrovirus that infection can induce erythroblastosis and myelocytomatosis, having the greatest pathogenicity and transmission ability within this class of viruses (3). ALV-J has caused significant economic loss in the poultry industry due to the increased tumors and mortality (4). ALV infection with a concomitant enhanced secondary infection is likely the outcome of immunosuppression (5, 6). ALV-J can be transmitted vertically and horizontally, and horizontal transmission of ALV-J is more efficient than other ALV subgroups (7). ALV-J infection has caused substantial economic losses in the poultry industry worldwide (8–10).
N6-methyladenosine (m6A), the methylation modification at the sixth N atom of adenine, is the most common post-transcriptional modification on mRNA, mediating over 60% RNA methylation (11, 12). The modification process of m6A RNA is mainly manipulated by three categories of proteins, including “writers”, “erasers” and “reader”, and exert their function by either directly being recognized by m6A-binding proteins or readers or indirectly by tuning the structure of the modified RNA to regulate RNA reader-protein interactions (13). RNA modification by m6A results in changes in properties such as charge, base pairing, secondary structure and protein-RNA interactions, which in turn affect the transport, localization, translation, and degradation of RNA, ultimately performing the function of regulating gene expression (14, 15). In recent years, many studies demonstrated that m6A methylation in mRNA is significantly associated with tumor proliferation, migration, invasion, and metastasis during cancer development and progression (16–18).
LncRNAs are a class of low protein-coding potential RNA with transcripts longer than 200 bp (19). LncRNAs play important roles in transcriptional, post-transcriptional regulation and chromatin modification by regulating gene expression, involving in a variety of biological processes in eukaryotes. Their abnormal expression is closely related to the malignancy of tumors, including tumor proliferation, differentiation, and apoptosis (20). Recent studies have reported that m6A-related lncRNAs are associated with the development of tumors and that aberrant lncRNA expression can serve as diagnostic markers in tumors (21). Therefore, exploring the function of lncRNA in ALV-J-induced tumors can help us better understand the pathogenesis of ALV-J.
Up to now, it has not been reported that m6A modification is involved in tumor formation induced by ALV-J in chickens. Given the indispensable role of RNA m6A modification in various tumors, it is reasonable to surmise that the dysregulation of m6A modification might be associated with the tumor formation induced by ALV-J infection in chickens. In this study, we reported the transcriptome-wide m6A analysis and the expression of lncRNAs in normal livers and ALV-J induced tumor livers in chickens, and a comprehensive analysis of m6A modifications of lncRNAs in ALV-J induced tumor liver tissues. We also analyzed the signaling pathways and potential regulated methyltransferases associated with the m6A modification of lncRNAs. These results indicated that m6A modification may play an important regulatory role in the tumor formation induced by ALV-J in chickens. We hope that this study will help to further investigate the potential role of lncRNAs with m6A modifications in the pathogenesis of ALV-J.
Materials and Methods
Animals and Tissue Collection
Specific pathogen-free (SPF) chickens were purchased from Guangdong Wen’s Food Group Co., Ltd. (Yunfu, China), fed in a negative-pressured biosafety isolation chamber with free water and commercial feed. 110 one-day-old SPF chickens were randomly divided into two groups, where each group contains 55 chickens. The first group (positive group, ALV-J induced tumor livers) was inoculated intraperitoneally with 103.7 TCID50/0.2 mL ALV-J strain NX0101 at one day of age; the other group (negative group, normal livers) was inoculated with the same volume of nutrient solution, and the negative group served as a control group. 110 days after infection, three chickens in each group were randomly euthanized and necropsied. The liver tissue samples were collected and then quickly frozen in liquid nitrogen and stored at -80°C environment.
The NX0101 strain of ALV-J used in the study was obtained from Professor Cui, Shandong Agricultural University, People’s Republic of China. This experiment strictly adhered to institutional and national guidelines for using and caring of laboratory animals. The use of animals in this study was approved by the South China Agricultural University Committee of Animal Experiments (approval ID: SYXK 2019-0136).
RNA Extraction and Quality Control
Total RNA of the livers was extracted using TRIzol Reagent (Invitrogen Corporation, CA, USA) according to the manufacturer’s protocols. Three biological replicates were used for positive and negative groups, respectively. The extracted RNA was first detected by PCR for ALV-J infection with the following primers (forward primer: 5′-AGCAACAAGCAAGAAAGACC-3′; reverse primer: 5′-CCGAACCAAAGGTAACACAC-3′). Then, the RNA quantification and quality assurance of each sample were qualified using Qubit 3.0 instrument. RNA integrity and gDNA contamination was identified by denaturing agarose gel electrophoresis.
MeRIP Sequencing and Analysis
Methylated RNA immunoprecipitation sequencing (MeRIP-seq) was performed on the basis of previously reported methods with some modifications (22). m6A RNA immunoprecipitation was performed with the GenSeq™ m6A RNA IP Kit (GenSeq Inc., China). Briefly, fragment RNA was incubated with an anti-m6A antibody (202003, Synaptic Systems, Germany) in immunoprecipitation buffer at 4°C for 2 hours. The mixture was further immunoprecipitated by incubation using protein A magnetic beads (Thermo Fisher, USA) at 4°C for 2 hours. Then, bound RNA on the magnetic beads was then eluted with N6-methyladenosine (BERRY & ASSOCIATES, PR3732), and then extracted with TRIzol reagent (Thermo Fisher, USA). Purified RNA was constructed by RNA-seq library construction using NEBNext® Ultra II Directional RNA Library Prep Kit (New England Biolabs, Inc., USA). After qualifying with Agilent 2100 Bioanalyzer system (Agilent Technologies Inc., CA, USA), the library sequencing was performed on Illumina NovaSeq 6000 sequencer with 150 bp paired-end.
After sequencing, paired-end reads were generated and quality controlled by Q30, followed by removing adaptors and low-quality reads using the Cutadapt program (v1.9.3) (23). Then, clean reads of input libraries were aligned to the reference genome (GCF_000002315.6_GRCg6a) using STAR, and circRNAs were identified by DCC software (24, 25). After that, clean reads of all libraries were mapped to the reference genome by Hisat2 (v2.0.4) (26). Differential methylation genes and sites were identified using Model-based Analysis of ChIP-Seq (MACS) and diffReps program, respectively (27, 28). These peaks identified by two software overlapping with exon of mRNA, lncRNA and circRNA were screened and annotated. Differentially expressed RNAs were determined by paired t-test algorithm. The thresholds for hyper- or hypomethylation were set as absolute fold change (FC) > 2.0 and P < 0.05. The gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of methylated associated genes using clusterProfiler package in R (29). The GO was categorized as three types: cellular component (CC), biological process (BP), and molecular function (MF). The P value < 0.05 denotes the statistically significant.
Transcriptome Library Construction and Sequencing
To investigate whether methylated modifications affect gene expression and their molecular mechanisms, we performed transcriptome sequencing for six samples and detected mRNAs, lncRNAs and circRNAs. Briefly, RNA was extracted using the TRIzol method, rRNAs were removed from total RNA using the NEBNext® rRNA Depletion Kit (New England Biolabs, Inc., Massachusetts, USA) following the manufacturer’s instructions, and then, RNA sequencing libraries were constructed using rRNA-depleted RNAs with TruSeq Stranded Total RNA Library Prep Kit (Illumina, USA) after quality control. Library quality control and quantification were performed using the BioAnalyzer 2100 system (Agilent Technologies, USA). Subsequent sequencing of 150 bp paired-end reads on the Illumina Novaseq 6000 sequencer. Raw data was quality controlled by Q30, and the Cutadapt program was used to remove adaptors and low-quality reads. The high-quality reads were aligned to the chicken reference genome (GCF_000002315.6_GRCg6a), and then mRNA, lncRNA and circRNA were detected and identified. Subsequently, differential expression analysis and functional enrichment analysis were performed.
Co-Express Analysis of m6A Methyltransferase and m6A-Related lncRNAs
The expression matrixes of 18 m6A-related methyltransferases were retrieved from the MeRIP-seq and RNA-seq, including ALKBH5, FTO, HNRNPA2B1, IGF2BP2, IGF2BP3, METTL3, METTL14, METTL16, RBM15, RBM15B, RBMX, YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, WTAP, and ZC3H13. 126 m6A-related lncRNAs were identified by Pearson’s correlation analysis. The process used the criteria of |Pearson R| > 0.3 and P < 0.05. The co-expression network was constructed by R, showing it as a Sankey plot.
Sequencing Data Validation by RT-qPCR and MeRIP-qPCR
Total RNA was extracted from normal livers and ALV-J induced tumor livers using the TRIzol method, and synthesized into cDNA with reverse transcription using PrimeScriptTM RT reagent Kit with gDNA Eraser (Takara, Japan). Then Reverse transcription quantitative real-time PCR (RT-qPCR) was performed using CFX96 Touch (Bio-Rad, USA) following the manufacturer’s instructions. The relative gene expression was normalized to GAPDH, then calculated using 2−ΔΔCT method. The primers using for RT-qPCR were listed in Table 1. For the detection of the fold enrichment of m6A level, the operation method is similar to RT-qPCR, and the sequences of primers using in MeRIP-qPCR were listed in Table 2.
Statistical Analysis
Data from three or more independent experiments were presented as the mean ± SEM. Statistical analysis were done using GraphPad Prism 8.0 software. Paired Student’s t-tests were performed between normal livers and ALV-J induced tumor livers. Differences with P < 0.05 were defined as the threshold for significance.
Results
Characterization of Normal Livers and ALV-J Induced Tumor Livers
In order to comprehensively analyze the differences of N6-methyladenosine RNA methylomes between normal livers and ALV-J induced tumor livers in chickens, liver tissues were collected from three randomly selected chickens in positive group (ALV-J induced tumor livers) and negative group (normal livers) at 110 days after infection. PCR results showed that the livers of the infected group were positive for ALV-J, while that of the uninfected group was negative for ALV-J (Figure 1A). Compared to the negative group, most chickens in the ALV-J infected group showed gradual emaciation. Tumor livers induced with ALV-J showed hepatomegaly and myeloma symptoms compared with the normal livers in the uninfected group (Figure 1B).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/868892/fimmu-13-868892-HTML/image_m/fimmu-13-868892-g001.jpg)
Figure 1 Characterization of normal livers and ALV-J induced tumor livers. (A) Electrophoretic analysis showed the infection between the two groups based on the PCR method. (B) Tumor liver induced by ALV-J infection in SPF chickens at 110 days. Negative represents chickens without ALV-J infection (normal livers); Positive represents chickens with ALV-J infection (tumor livers).
Transcriptome-Wide Detection of m6A Methylation in Normal Livers and ALV-J Induced Tumor Livers
MeRIP-seq generated 215,876,754 to 227,827,726 raw reads from IP or input samples from normal livers and ALV-J induced tumor livers. After filtering out low-quality data, more than 96% of high-quality reads from each sample were aligned to the Gallus reference genome. Over 86% of the clean reads from all the samples were uniquely mapped to the chicken reference genome (Supplementary Table S1). MeRIP-seq analysis of RNA isolated from chicken livers indicated that, 14,584 and 14,057 expressed mRNA transcripts were detected in positive (ALV-J induced tumor livers) and negative groups (normal livers), respectively. After mapping the methylated RNA fragments to the transcriptome, a total of 5,848 m6A sites were identified among the 3,597 coding transcripts in the positive group and 6,223 m6A sites were identified among 3,637 coding transcripts in the negative group (Supplementary Table S2). The proportion of methylated transcripts was 40.10% and 44.27% in positive and negative groups, respectively.
To determine how the m6A modification distributed throughout the chicken transcriptome, we classified the methylated transcripts based on the number of m6A peaks contained in each transcript. We found that the number of m6A modified sites varies among different transcripts, while more than 83% of modified transcripts contained one or two m6A peaks, and about 7.69% of the methylated transcripts contained four or more peaks (Figure 2A). In total, the chicken liver transcriptome on average contains 1.71 and 1.61 m6A peaks per methylated transcript in positive and normal negative groups, respectively.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/868892/fimmu-13-868892-HTML/image_m/fimmu-13-868892-g002.jpg)
Figure 2 Overview of m6A methylation in the normal livers and ALV-J induced tumor livers. (A) The proportion of genes harboring different numbers of m6A peaks in two groups. (B, C) The percentage of m6A peaks in five nonoverlapping segments of transcripts in two groups. startC, start codon; CDS, coding sequence; stopC, stop codon; 3’ UTR, the three prime untranslated regions; 5’ UTR, the five prime untranslated regions. (D) Accumulation of m6A peaks along with transcript in positive and negative groups, respectively. Each transcript is divided into three parts, including 5’ UTR, CDS, and 3’ UTR. (E, F) The top three motifs enriched across m6A peaks in positive (E) and negative group (F). Negative represents chickens without ALV-J infection (normal livers); Positive represents chickens with ALV-J infection (tumor livers).
To determine the distributed localization of m6A in the mRNA transcripts, m6A peaks were divided into five non-overlapping segments: 5′UTR, start codon segment (startC), coding sequence (CDS), stop codon segment (stopC) and 3′UTR (Figures 2B, C). Our results showed that m6A was most commonly present in CDS, with some located near the start and stop codons, and some differences between the two experimental groups, which is consistent with the patterns identified in the mouse, pig and chicken (30–32). The m6A peaks near the startC were 20.09% and 19.86% in positive and negative samples, respectively. m6A peaks near the CDS reduced 1.62% from negative to the positive group. Furthermore, the abundance of m6A peaks near the stopC increased 1.01% in the positive group (ALV-J induced tumor livers). The distribution of m6A across the transcriptome was verified by m6A reads along the transcripts. Consistent with the distribution of m6A peaks, m6A reads were distributed throughout the mRNA transcript, where reads were elevated in the CDS and peaked at the 3′UTR. Metagene profiling of the m6A peaks showed that they were primarily enriched in CDS, near the start and stop codons, and close to the beginning of 3′UTR and the ending of 5′UTR, which differs from the pattern identified in mammals, and same to chicken (Figure 2D) (22, 30). Together, the results reveal that m6A is dynamic in ALV-J induced tumor livers in chickens.
Motif enrichment analysis of 1000 peaks within mRNAs with the highest scores (-10×log10(p-value)) obtained from three biological replicates (1,000 peaks per replicate) revealed the conservative RRACH motif (R represents purine, A represents m6A, and H represents a nonguanine base), as well as other motifs, in the positive and negative groups, respectively. The results showed that GGACU was significantly enriched and consistently considered to be the beat motif in both groups, suggesting that the RRACH motif employed in the genesis of ALV-J is conserved in the presence or absence of ALV-J infection, which means that GGACU is a commonly present m6A modification sequence (Figures 2E, F).
Distribution of Differentially Methylated m6A Sites in Normal Livers and ALV-J Induced Tumor Livers
To investigate the differences in m6A modifications between the positive group (ALV-J induced tumor livers) and negative group (normal livers), we initially identified m6A peaks and genes specific to both groups. We found 6,541 peaks that were methylated in both positive and negative groups, and 2,222 and 1,651 peaks that were specifically hypermethylated and hypomethylated in the mRNA of positive against negative group, respectively (Figure 3A). Of which, 68.48% (4,858/7,094) were significantly hypermethylation among the differentially expressed methylation sites (positive vs. negative, Figure 3B). Similarly, we screened 99 methylation peaks in lncRNAs that were present simultaneously in two groups, and 67 specific hypermethylated peaks and 50 specific hypomethylated peaks were identified in positive group against negative group (Supplementary Figure S1). We also identified 561 differentially methylated m6A sites within 172 lncRNA transcripts, of which 71.30% (400/561) were hypermethylation sites (positive vs. negative). A total of 138 differentially methylated m6A sites within 95 circRNA genes, of which 73.19% (101/138) were hypermethylated (positive vs. negative). Table 3 showed the top 20 hyper- and hypomethylated m6A sites within mRNAs with the highest values of fold change, of which were >65-fold.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/868892/fimmu-13-868892-HTML/image_m/fimmu-13-868892-g003.jpg)
Figure 3 Global m6A modification changes in the normal livers and ALV-J induced tumor livers. (A) Venn diagram of m6A methylation sites identified in mRNAs from two groups. (B) The Volcano plot showed the identification of significant hyper- and hypomethylated m6A peaks. (C, D) GO enrichment analysis of biological processes involved in hypermethylated (C) and hypomethylated (D) m6A genes in ALV-J infected positive group compared with negative group. (E, F) Pathway analysis of hypermethylated (E) and hypomethylated (F) m6A genes in mRNAs of ALV-J infected positive group compared to negative group. Negative represents chickens without ALV-J infection (normal livers); Positive represents chickens with ALV-J infection (tumor livers).
Differentially Methylated mRNAs Involved in Important Biological Pathways in ALV-J Induced Tumor Livers
To reveal the function of m6A in ALV-J induced tumor livers in chickens, protein-coding genes containing differentially methylated m6A sites were collected for GO enrichment and KEGG pathway analysis. As for the BP category of mRNA, genes with hypermethylated m6A sites were significantly enriched in immune related terms (13/193) such as regulation of immune response/process, immune system development, adaptive immune response, while genes with hypomethylated m6A sites were highly involved in regulation of cell death/apoptosis in ALV-J induced tumor livers (7/86, P < 0.05, Figures 3C, D). For the CC category, hypermethylated m6A genes were mainly enriched in membrane and chromosome. For the MF category, hypermethylation of m6A in the positive group was notably enriched in transferase activity, kinase regulator activity (Supplementary Table S3). Compared with the negative group of chickens which were not infected by ALV-J, hypermethylation in tumor livers of chickens were more involved in regulating the immune system, while hypomethylation in tumor livers of chickens were more participated in cell apoptosis. Furthermore, genes with hypermethylated m6A sites were found to be significant in associated pathways, such as Toll-like receptor, MAPK, apoptosis, VEGF, and JAK-STAT signaling pathway, while hypomethylated m6A sites were highly related to PPAR and Wnt signaling pathway (P < 0.05, Figures 3E, F).
Overview of RNA Expression Profiles in Normal Livers and ALV-J Induced Tumor Livers
To understand whether m6A modifications could affect gene expression in ALV-J induced tumor livers, RNA-seq data were firstly used to investigate the differential expression of genes between the two groups. The expression patterns of RNA were investigated by RNA-seq, from which the expression level of 15,363 mRNAs, 31,279 lncRNAs, and 4,705 circRNAs were found in the six samples, respectively. Differently expressed genes between the positive (ALV-J induce tumor livers) and negative (normal livers) groups were further screened. Compared with the negative group, 1,896 mRNAs, 339 lncRNAs, and 43 circRNAs with differential expression in ALV-J induced tumor livers were obtained, including 1,687 up-regulated mRNAs, 295 up-regulated lncRNAs, and 29 up-regulated circRNAs, respectively. The hierarchical clustering and scatter plot of the RNA-seq data were shown in Figures 4A, B, indicating that these RNAs have different expression patterns in two groups.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/868892/fimmu-13-868892-HTML/image_m/fimmu-13-868892-g004.jpg)
Figure 4 Screening of differentially expressed genes in normal livers and ALV-J induced tumor livers by RNA-seq. (A) Heatmap of RNA-seq data in two groups. Rows: mRNAs; columns: ALV-J infected positive and negative groups. Red, black, and green indicate the up-regulation, unchanged, and down-regulation of mRNAs, respectively. (B) Differentially expressed mRNAs in ALV-J infected positive and negative groups. Red dots represent up-regulated genes in the positive group compared with the negative group, green dots represent down-regulated genes, and blue dots represent genes with no significant differences. (C) GO enrichment of biological processes involved in up-regulated genes in the positive group compared to the negative group. (D) Pathway analysis of up-regulated genes in the positive group. Negative represents chickens without ALV-J infection (normal livers); Positive represents chickens with ALV-J infection (tumor livers).
To further understand the effects of lncRNAs, we annotated the differentially expressed lncRNAs and analyzed their GO enrichment and KEGG pathways. GO analysis revealed that positive up-regulated genes in ALV-J induced tumor livers were significantly enriched in biological processes involving regulation of the cellular process, regulation of cellular metabolic process, and negative regulation of RNA biosynthetic process (Figure 4C). Moreover, KEGG analysis showed that MAPK, cytokine-cytokine receptor interaction, Toll-like receptor, the intestinal immune network for IgA production, and Wnt pathways were significantly altered in up-regulated lncRNAs (Figure 4D). Together, these results indicated that m6A modification might influence RNA expression, thereby regulating the tumor formation induced by ALV-J infection.
Screening the Common Differential Genes Between lncRNAs Methylation and Gene Expression of lncRNAs in ALV-J Induced Tumor Livers
To obtain differentially expressed lncRNAs shared in MeRIP and RNA-seq, a few lncRNAs with both differential methylation and expression levels in ALV-J-induced tumor livers were identified in two sequencing results. Thus, 17 differentially m6A-methylated lncRNAs were screened (Figure 5A). Besides, their associated genes were found, including HDAC9, CMPK2, MINDY4B, ATOH8, SOX7, LOC100859478, etc. Among them, five lncRNAs (LOC107052235, ATOH8, LOC101751658, SOX7, and LOC100859478) were hypomethylated and low expressed in ALV-J induced tumor livers, while the other lncRNAs like HDAC9, CMPK2, MINDY4B, etc., was hypermethylated and highly expressed in ALV-J induced tumor livers. The details were shown in Table 4.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/868892/fimmu-13-868892-HTML/image_m/fimmu-13-868892-g005.jpg)
Figure 5 Conjoint analysis of m6A methyltransferases and lncRNA genes in ALV-J induced tumor livers. (A) Venn diagram showing the relationship between m6A methylation and lncRNA expression level. (B) Positive correlation between m6A methylation and mRNA expression level (Pearson R=0.75; p<1e-10). (C) Sankey relational diagram for 18 m6A genes and 126 m6A-related lncRNAs.
Co-Expression Analysis of m6A Methyltransferases and lncRNAs in ALV-J Induced Tumor Livers
By cross-analysis of the MeRIP-seq and RNA-seq data, we discovered a positive correlation of differentially methylated m6A sites and gene expression levels in normal livers and ALV-J induced tumor livers (Figure 5B). To explore the potential target lncRNAs associated with differentially expressed m6A methyltransferases in ALV-J induced tumor livers, we analyzed the expression correlation of all differentially expressed mRNAs with them. According to the previous studies (33), the expression matrixes of 18 m6A methyltransferases were obtained from the MeRIP-seq and RNA-seq data, including the expression of writers METTL3, METTL14, METTL16, RBM15, RBM15B, ZC3H13, and WTAP), readers (IGF2BP2, IGF2BP3, YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPA2B1, and RBMX), and erasers (ALKBH5, and FTO). After integrating these genes, a co-expression network was constructed for 18 m6A methyltransferases and 126 differentially expressed genes of lncRNAs to visually manifest the relationship (Figure 5C, Supplementary Table S4). From the Sankey diagram, it was clear that RBM15, YTHDC1, and YTHDC2 were associated with more lncRNAs ― over 37% (147/397) ― indicating that these methylation-related enzymes have important catalytic and recognition roles in the occurrence of the tumor formation induced by ALV-J infection. Likewise, we detected expressed genes LOC112530919, RAD50, LOC107054512, FBXO40, ARGLU1, ELOC, RBM39, KHDRBS1, and so on, that were connected with more methylation signals, involved in the methylation process of m6A modifications as major factors during disease onset and progression. Together, m6A modification can regulate the function of lncRNAs through the binding of reader proteins, which can affect the progression of tumor formation.
Differentially Expressed lncRNAs and m6A Modification Levels of lncRNAs Were Confirmed by RT-qPCR and MeRIP-qPCR
To confirm the accuracy of RNA-seq, five lncRNAs were selected and confirmed by RT-qPCR. RT-qPCR analysis revealed that the lncRNA levels of HDAC9, CMPK2, MINDY4B were up-regulated in the ALV-J induced tumor livers (positive group) compared to the normal livers in the negative group, whereas lncRNA levels of ATOH8 and LOC100859478 were down-regulated in the positive group (Figure 6A). Furthermore, to confirm the accuracy of MeRIP-seq, we selected three transcripts (HDAC9, MINDY4B, SOX7) and the m6A modification in these three transcripts were further verified by MeRIP-qPCR. MeRIP-qPCR results showed that the relative level of m6A modification in HDAC9, MINDY4B was significantly increased in ALV-J induced tumor livers, while the relative level of m6A modification of SOX7 in ALV-J induced tumor livers was significantly decreased (Figure 6B). These results of RT-qPCR and MeRIP-qPCR were consistent with the data of RNA-seq (Figure 6C) and MeRIP-seq (Figure 6D), indicating that the sequencing data and the analysis of this study are credible and valid. Thus, these data suggested that m6A mediates dynamic gene expression of ALV-J induced tumor formation.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/868892/fimmu-13-868892-HTML/image_m/fimmu-13-868892-g006.jpg)
Figure 6 Validation of differentially expressed lncRNAs and m6A modification levels of lncRNAs in normal livers and ALV-J induced tumor livers using RT-qPCR and MeRIP-qPCR. (A) The relative lncRNA expression level of CMPK2, HDAC9, MINDY4B, ATOH8, and LOC100859478 in ALV-J induced tumor livers (positive) and normal livers (negative) was detected by RT-qPCR. (B) The m6A level of HDAC9, MINDY4B, and SOX7 was validated by MeRIP-qPCR. Expression levels were normalized by GAPDH. (C) The lncRNA fold change of CMPK2, HDAC9, MINDY4B, ATOH8, and LOC100859478 in ALV-J induced tumor livers (positive) and normal livers (negative) in RNA-Seq results. (D) The m6A level of HDAC9, MINDY4B, and SOX7 in MeRIP-seq results. Data are presented as mean ± SEM. **P < 0.05; *P < 0.01. Negative represents chickens without ALV-J infection (normal livers); Positive represents chickens with ALV-J infection (tumor livers).
Discussion
ALV is a highly infectious and oncogenic disease that threatens poultry farming countries worldwide (34). Like HIV, which belongs to the same retrovirus, different subtypes of ALV have different clinical features, and the pathogenesis of ALV-J is particularly complex (35, 36). Epigenetic changes associated with RNA reversible chemical modifications play an essential role in the life cycle of the virus. A growing number of evidence discusses that m6A modifications are closely associated with the development and progression of different cancers, anticipating the use of m6A methylation as a new cancer diagnostic biomarker and therapeutic target (33, 37). Han et al. discovered that METTL3 affects the proliferation of bladder cancer by regulating m6A modifications in non-coding RNAs (38). It was found that dynamic changes in m6A of lncRNA play a vital regulatory role in chicken Marek’s disease virus replication by constructing a transcriptome-wide m6A modification map of Marek’s disease (19). With these examples, we suspected that m6A modifications also contribute to tumor formation in ALV-J infected chickens by affecting gene expression. Here, we report that the dynamic presence of m6A modifications in chicken liver transcripts with and without ALV-J infection affects gene expression.
Infection with ALV-J leads to dynamic changes in m6A and accumulation in tumor liver tissues, providing evidence for disease-induced changes in m6A status. In this study, a chicken liver model with or without ALV-J infection was established to assess m6A modification and revealed its differences in the tumor livers of chickens after ALV-J infection, supporting the dynamic characterization of m6A modifications. To our knowledge, this is the first report about the high-throughput study of RNA methylation in liver tissues of chickens infected with ALV-J. The data showed that extensive mRNA hypermethylation and hypomethylation occurred in liver tissues during the occurrence of tumor formation induced by ALV-J infection, further suggesting that m6A modification may influence ALV-J induced tumor formation by regulating gene expression.
LncRNA is a type of transcript with a length of more than 200 nucleotides and has no protein-coding function, which significantly impacts chromatin organization, transcription, and post-transcriptional regulation (39, 40). To investigate the effect of m6A modification in lncRNA on tumor formation induced by ALV-J infection, our results showed an increase in m6A modification of lncRNA in ALV-J induced tumor livers. The m6A modifications were positively correlated with mRNA expression in the tumor liver samples, and we further confirmed that there was a close relationship between lncRNA and m6A modification in the formation of tumor livers induced by ALV-J infection. m6A is prevalent in mRNA transcripts from chicken liver tissue, with m6A reaching its highest value near stop codon and then decreasing at the 3′UTR. This distribution pattern is consistent with previous findings in mice (31). GO and KEGG analysis identified the differential expression of m6A and lncRNAs associated with specific signaling pathways in ALV-J induced tumor livers. This study predicted the Toll-like, MAPK, cytokine-cytokine receptor interaction, ErbB pathways (ALV-J induced tumor livers vs. normal livers) and found that lncRNAs with m6A modifications were common in tumor livers induced by ALV-J infection; hence, we suspected that lncRNAs with m6A modifications might mediate tumor formation and mediate viral immunosuppression. It was reported that high expression of Toll-like receptors contributes to epithelial ovarian carcinogenesis, inhibits the development of bladder cancer, and serves as a potential therapeutic target for colorectal cancer (41, 42). Similarly, the ErbB signaling pathway regulates human malignancies like breast, lung, and bladder cancer, allowing for cancer prevention and treatment (43–45). Nevertheless, m6A differences were investigated to reveal potential effects. However, further analysis is still needed to validate these results. These suggest that m6A-modified lncRNAs affect the occurrence and development of tumor formation inducted by ALV-J infection through biological processes, cellular composition, molecular functions, and signaling pathways.
According to the correlation analysis of MeRIP-seq and RNA-seq, we identified 17 m6A-modified lncRNAs, most of which are present in the intergenic region, which is very similar to the definition of lncRNAs that cannot encode proteins but may cause functional changes in mRNAs through structural changes (40, 46). Furthermore, we performed co-expression network construction using 18 methyltransferases and 126 lncRNAs. We identified the crucial functional factor and key lncRNAs, which will provide a theoretical basis for future studies on m6A modification of lncRNA to induce tumor formation induced by ALV-J infection in chickens. Regulation of m6A modifications may become a strategy for treating the disease induced by ALV-J in chickens in the future. Currently, there are no reports of m6A modification of lncRNAs on HIV. We trust that the results of this study may provide new insights into the pathogenesis of HIV, especially that m6A modification of lncRNAs may be a factor affecting HIV pathogenesis and explore possible prevention and control methods.
Overall, this study used transcriptomic and m6A methylationome data to integrally analyze differential lncRNA with m6A modifications in ALV-J-induced tumor livers. Comparing ALV-J-induced tumor livers with normal livers, we identified a few functional lncRNAs and methyltransferases, revealing a tight relationship between m6A methylation and chicken tumor formation induced by ALV-J infection. Therefore, the results of this study provide a new vision and theoretical basis for the prevention and control of subgroup J avian leukosis.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, GSE192892.
Ethics Statement
This experiment strictly adhered to institutional and national guidelines for using and caring of laboratory animals. The use of animals in this study was approved by the South China Agricultural University Committee of Animal Experiments (approval ID: SYXK 2019-0136).
Author Contributions
QX and XZ designed the study. QZ and XZ wrote the manuscript. ZY, LC, YH, and ZX performed the experiments. QZ and ZY analyzed the results. HZ, WL, FC, and QX supervised the project and critically revised the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by the Key Research and Development Program of Guangdong Province (2020B020222001), the National Natural Science Foundation of China (Grant no. 31902252, 31972659), Guangdong Basic and Applied Basic Research Foundation (2019A1515012006), Chief expert Project of Agricultural Industry Technology system in Guangdong Province (2019KJ128), and Special Project of National Modern Agricultural Industrial Technology System (Grant no. CARS-41).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
We sincerely appreciate members of Cloudseq Biotech Inc. (Shanghai, China) who provided in sequencing and data analysis.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.868892/full#supplementary-material
References
1. Cui Z, Sun S, Zhang Z, Meng S. Simultaneous Endemic Infections With Subgroup J Avian Leukosis Virus and Reticuloendotheliosis Virus in Commercial and Local Breeds of Chickens. Avian Pathol (2009) 38:443–8. doi: 10.1080/03079450903349188
2. Zhou D, Xue J, Zhang Y, Wang G, Feng Y, Hu L, et al. Outbreak of Myelocytomatosis Caused by Mutational Avian Leukosis Virus Subgroup J in China, 2018. Transbound Emerg Dis (2019) 66:622–6. doi: 10.1111/tbed.13096
3. Li Y, Cui S, Li W, Wang Y, Cui Z, Zhao P, et al. Vertical Transmission of Avian Leukosis Virus Subgroup J (ALV-J) From Hens Infected Through Artificial Insemination With ALV-J Infected Semen. BMC Vet Res (2017) 13:204. doi: 10.1186/s12917-017-1122-4
4. Dong X, Zhao P, Xu B, Fan J, Meng F, Sun P, et al. Avian Leukosis Virus in Indigenous Chicken Breeds, China. Emerg Microbes Infect (2015) 4:e76. doi: 10.1038/emi.2015.76
5. Liu D, Dai M, Zhang X, Cao W, Liao M. Subgroup J Avian Leukosis Virus Infection of Chicken Dendritic Cells Induces Apoptosis via the Aberrant Expression of microRNAs. Sci Rep (2016) 6:20188. doi: 10.1038/srep20188
6. Wang F, Wang X, Chen H, Liu J, Cheng Z. The Critical Time of Avian Leukosis Virus Subgroup J-Mediated Immunosuppression During Early Stage Infection in Specific Pathogen-Free Chickens. J Vet Sci (2011) 12:235–41. doi: 10.4142/jvs.2011.12.3.235
7. Lin Y, Xia J, Zhao Y, Wang F, Yu S, Zou N, et al. Reproduction of Hemangioma by Infection With Subgroup J Avian Leukosis Virus: The Vertical Transmission is More Hazardous Than the Horizontal Way. Virol J (2013) 10:97. doi: 10.1186/1743-422X-10-97
8. Silva RF, Fadly AM, Hunt HD. Hypervariability in the Envelope Genes of Subgroup J Avian Leukosis Viruses Obtained From Different Farms in the United States. Virology (2000) 272:106–11. doi: 10.1006/viro.2000.0352
9. Thapa BR, Omar AR, Arshad SS, Hair-Bejo M. Detection of Avian Leukosis Virus Subgroup J in Chicken Flocks From Malaysia and Their Molecular Characterization. Avian Pathol (2004) 33:359–63. doi: 10.1080/0307945042000220435
10. Xu B, Dong W, Yu C, He Z, Lv Y, Sun Y, et al. Occurrence of Avian Leukosis Virus Subgroup J in Commercial Layer Flocks in China. Avian Pathol (2004) 33:13–7. doi: 10.1080/03079450310001636237a
11. Du K, Zhang L, Lee T, Sun T. M(6)A RNA Methylation Controls Neural Development and Is Involved in Human Diseases. Mol Neurobiol (2019) 56:1596–606. doi: 10.1007/s12035-018-1138-1
12. Wang S, Sun C, Li J, Zhang E, Ma Z, Xu W, et al. Roles of RNA Methylation by Means of N(6)-Methyladenosine (M(6)A) in Human Cancers. Cancer Lett (2017) 408:112–20. doi: 10.1016/j.canlet.2017.08.030
13. Brocard M, Ruggieri A, Locker N. M6a RNA Methylation, A New Hallmark in Virus-Host Interactions. J Gen Virol (2017) 98:2207–14. doi: 10.1099/jgv.0.000910
14. Du Y, Hou G, Zhang H, Dou J, He J, Guo Y, et al. SUMOylation of the M6a-RNA Methyltransferase METTL3 Modulates its Function. Nucleic Acids Res (2018) 46:5195–208. doi: 10.1093/nar/gky156
15. Zhao BS, Roundtree IA, He C. Post-Transcriptional Gene Regulation by mRNA Modifications. Nat Rev Mol Cell Biol (2017) 18:31–42. doi: 10.1038/nrm.2016.132
16. Liu ZX, Li LM, Sun HL, Liu SM. Link Between M6a Modification and Cancers. Front Bioeng Biotechnol (2018) 6:89. doi: 10.3389/fbioe.2018.00089
17. Zhao W, Qi X, Liu L, Ma S, Liu J, Wu J. Epigenetic Regulation of M(6)A Modifications in Human Cancer. Mol Ther Nucleic Acids (2020) 19:405–12. doi: 10.1016/j.omtn.2019.11.022
18. Liu S, Li Q, Chen K, Zhang Q, Li G, Zhuo L, et al. The Emerging Molecular Mechanism of M(6)A Modulators in Tumorigenesis and Cancer Progression. BioMed Pharmacother (2020) 127:110098. doi: 10.1016/j.biopha.2020.110098
19. Sun A, Zhu X, Liu Y, Wang R, Yang S, Teng M, et al. Transcriptome-Wide N6-Methyladenosine Modification Profiling of Long Non-Coding RNAs During Replication of Marek's Disease Virus In Vitro. BMC Genomics (2021) 22:296. doi: 10.1186/s12864-021-07619-w
20. Yu J, Mao W, Sun S, Hu Q, Wang C, Xu Z, et al. Identification of an M6a-Related lncRNA Signature for Predicting the Prognosis in Patients With Kidney Renal Clear Cell Carcinoma. Front Oncol (2021) 11:663263. doi: 10.3389/fonc.2021.663263
21. Wang H, Meng Q, Ma B. Characterization of the Prognostic M6a-Related lncRNA Signature in Gastric Cancer. Front Oncol (2021) 11:630260. doi: 10.3389/fonc.2021.630260
22. Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR. Comprehensive Analysis of mRNA Methylation Reveals Enrichment in 3' UTRs and Near Stop Codons. Cell (2012) 149:1635–46. doi: 10.1016/j.cell.2012.05.003
23. Martin M. Cutadapt Removes Adapter Sequences From High-Throughput Sequencing Reads. EMBnet J (2011) 17.1:10–2. doi: 10.14806/ej.17.1.200
24. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: Ultrafast Universal RNA-Seq Aligner. Bioinformatics (2013) 29:15–21. doi: 10.1093/bioinformatics/bts635
25. Cheng J, Metge F, Dieterich C. Specific Identification and Quantification of Circular RNAs From Sequencing Data. Bioinformatics (2016) 32:1094–6. doi: 10.1093/bioinformatics/btv656
26. Kim D, Langmead B, Salzberg SL. HISAT: A Fast Spliced Aligner With Low Memory Requirements. Nat Methods (2015) 12:357–60. doi: 10.1038/nmeth.3317
27. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-Based Analysis of ChIP-Seq (MACS). Genome Biol (2008) 9:R137. doi: 10.1186/gb-2008-9-9-r137
28. Shen L, Shao NY, Liu X, Maze I, Feng J, Nestler EJ. Diffreps: Detecting Differential Chromatin Modification Sites From ChIP-Seq Data With Biological Replicates. PloS One (2013) 8:e65598. doi: 10.1371/journal.pone.0065598
29. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. Clusterprofiler 4.0: A Universal Enrichment Tool for Interpreting Omics Data. Innovation (N Y) (2021) 2:100141. doi: 10.1016/j.xinn.2021.100141
30. Tao X, Chen J, Jiang Y, Wei Y, Chen Y, Xu H, et al. (6) -Methyladenosine Methylome Profiling of Porcine Muscle and Adipose Tissues Reveals a Potential Mechanism for Transcriptional Regulation and Differential Methylation Pattern. BMC Genomics (2017) 18:336. doi: 10.1186/s12864-017-3719-1
31. Luo Z, Zhang Z, Tai L, Zhang L, Sun Z, Zhou L. Comprehensive Analysis of Differences of N(6)-Methyladenosine RNA Methylomes Between High-Fat-Fed and Normal Mouse Livers. Epigenomics-Uk (2019) 11:1267–82. doi: 10.2217/epi-2019-0009
32. Cheng B, Leng L, Li Z, Wang W, Jing Y, Li Y, et al. Profiling of RNA N (6) -Methyladenosine Methylation Reveals the Critical Role of M(6)A in Chicken Adipose Deposition. Front Cell Dev Biol (2021) 9:590468. doi: 10.3389/fcell.2021.590468
33. He L, Li H, Wu A, Peng Y, Shu G, Yin G. Functions of N6-Methyladenosine and Its Role in Cancer. Mol Cancer (2019) 18:176. doi: 10.1186/s12943-019-1109-9
34. Lee HJ, Park KJ, Lee KY, Yao Y, Nair V, Han JY. Sequential Disruption of ALV Host Receptor Genes Reveals No Sharing of Receptors Between ALV Subgroups A, B and J. J Anim Sci Biotechnol (2019) 10:23. doi: 10.1186/s40104-019-0333-x
35. Yu H, Xu H, Yan C, Zhu S, Lan X, Lu Y, et al. gga-miR-148a-5p-Targeting PDPK1 Inhibits Proliferation and Cell Cycle Progression of Avain Leukosis Virus Subgroup J (ALV-J)-Infected Cells. Front Cell Dev Biol (2020) 8:587889. doi: 10.3389/fcell.2020.587889
36. Swanstrom R, Graham WD, Zhou S. Sequencing the Biology of Entry: The Retroviral Env Gene. Curr Top Microbiol Immunol (2017) 407:65–82. doi: 10.1007/82_2017_35
37. Sun T, Wu R, Ming L. The Role of M6a RNA Methylation in Cancer. BioMed Pharmacother (2019) 112:108613. doi: 10.1016/j.biopha.2019.108613
38. Han J, Wang JZ, Yang X, Yu H, Zhou R, Lu HC, et al. METTL3 Promote Tumor Proliferation of Bladder Cancer by Accelerating Pri-Mir221/222 Maturation in M6a-Dependent Manner. Mol Cancer (2019) 18:110. doi: 10.1186/s12943-019-1036-9
39. Ma S, Chen C, Ji X, Liu J, Zhou Q, Wang G, et al. The Interplay Between M6a RNA Methylation and Noncoding RNA in Cancer. J Hematol Oncol (2019) 12:121. doi: 10.1186/s13045-019-0805-7
40. Dai D, Wang H, Zhu L, Jin H, Wang X. N6-Methyladenosine Links RNA Metabolism to Cancer Progression. Cell Death Dis (2018) 9:124. doi: 10.1038/s41419-017-0129-x
41. Ohadian MS, Nowroozi MR. Toll-Like Receptors: The Role in Bladder Cancer Development, Progression and Immunotherapy. Scand J Immunol (2019) 90:e12818. doi: 10.1111/sji.12818
42. Jiang Y, Wan Y, Gong M, Zhou S, Qiu J, Cheng W. RNA Demethylase ALKBH5 Promotes Ovarian Carcinogenesis in a Simulated Tumour Microenvironment Through Stimulating NF-kappaB Pathway. J Cell Mol Med (2020) 24:6137–48. doi: 10.1111/jcmm.15228
43. Britten CD. Targeting ErbB Receptor Signaling: A Pan-ErbB Approach to Cancer. Mol Cancer Ther (2004) 3:1335–42. doi: 10.1158/1535-7163.1335.3.10
44. Normanno N, De Luca A, Bianco C, Strizzi L, Mancino M, Maiello MR, et al. Epidermal Growth Factor Receptor (EGFR) Signaling in Cancer. Gene (2006) 366:2–16. doi: 10.1016/j.gene.2005.10.018
45. Wang Z. ErbB Receptors and Cancer. Methods Mol Biol (2017) 1652:3–35. doi: 10.1007/978-1-4939-7219-7_1
Keywords: ALV-J, m6A, MeRIP-seq, RNA-seq, co-expressed
Citation: Zhao Q, Yao Z, Chen L, He Y, Xie Z, Zhang H, Lin W, Chen F, Xie Q and Zhang X (2022) Transcriptome-Wide Dynamics of m6A Methylation in Tumor Livers Induced by ALV-J Infection in Chickens. Front. Immunol. 13:868892. doi: 10.3389/fimmu.2022.868892
Received: 03 February 2022; Accepted: 28 March 2022;
Published: 22 April 2022.
Edited by:
Aijun Sun, Texas A&M University, United StatesReviewed by:
Cao Yong Chang, Sun Yat-sen University, ChinaShuang Chang, Shandong Agricultural University, China
Copyright © 2022 Zhao, Yao, Chen, He, Xie, Zhang, Lin, Chen, Xie and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Qingmei Xie, qmx@scau.edu.cn; Xinheng Zhang, xhzhang@scau.edu.cn
†These authors have contributed equally to this work and share first authorship