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 10^3.7 TCID₅₀/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).

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.

Methylated RNA immunoprecipitation sequencing (MeRIP-seq) was performed on the basis of previously reported methods with some modifications (22). m⁶A RNA immunoprecipitation was performed with the GenSeq™ m⁶A RNA IP Kit (GenSeq Inc., China). Briefly, fragment RNA was incubated with an anti-m⁶A 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.

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.

The expression matrixes of 18 m⁶A-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 m⁶A-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.

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 m⁶A level, the operation method is similar to RT-qPCR, and the sequences of primers using in MeRIP-qPCR were listed in Table 2 .

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.

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 ).

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 m⁶A sites were identified among the 3,597 coding transcripts in the positive group and 6,223 m⁶A 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 m⁶A modification distributed throughout the chicken transcriptome, we classified the methylated transcripts based on the number of m⁶A peaks contained in each transcript. We found that the number of m⁶A modified sites varies among different transcripts, while more than 83% of modified transcripts contained one or two m⁶A 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 m⁶A peaks per methylated transcript in positive and normal negative groups, respectively.

To determine the distributed localization of m⁶A in the mRNA transcripts, m⁶A 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 m⁶A 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 m⁶A peaks near the startC were 20.09% and 19.86% in positive and negative samples, respectively. m⁶A peaks near the CDS reduced 1.62% from negative to the positive group. Furthermore, the abundance of m⁶A peaks near the stopC increased 1.01% in the positive group (ALV-J induced tumor livers). The distribution of m⁶A across the transcriptome was verified by m⁶A reads along the transcripts. Consistent with the distribution of m⁶A peaks, m⁶A reads were distributed throughout the mRNA transcript, where reads were elevated in the CDS and peaked at the 3′UTR. Metagene profiling of the m⁶A 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 m⁶A 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 m⁶A, 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 m⁶A modification sequence ( Figures 2E, F ).

To investigate the differences in m⁶A modifications between the positive group (ALV-J induced tumor livers) and negative group (normal livers), we initially identified m⁶A 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 m⁶A sites within 172 lncRNA transcripts, of which 71.30% (400/561) were hypermethylation sites (positive vs. negative). A total of 138 differentially methylated m⁶A 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 m⁶A sites within mRNAs with the highest values of fold change, of which were >65-fold.

To reveal the function of m⁶A in ALV-J induced tumor livers in chickens, protein-coding genes containing differentially methylated m⁶A sites were collected for GO enrichment and KEGG pathway analysis. As for the BP category of mRNA, genes with hypermethylated m⁶A 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 m⁶A 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 m⁶A genes were mainly enriched in membrane and chromosome. For the MF category, hypermethylation of m⁶A 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 m⁶A sites were found to be significant in associated pathways, such as Toll-like receptor, MAPK, apoptosis, VEGF, and JAK-STAT signaling pathway, while hypomethylated m⁶A sites were highly related to PPAR and Wnt signaling pathway (P < 0.05, Figures 3E, F ).

To understand whether m⁶A 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.

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 m⁶A modification might influence RNA expression, thereby regulating the tumor formation induced by ALV-J infection.

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 m⁶A-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 .

By cross-analysis of the MeRIP-seq and RNA-seq data, we discovered a positive correlation of differentially methylated m⁶A 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 m⁶A 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 m⁶A 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 m⁶A 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 m⁶A modifications as major factors during disease onset and progression. Together, m⁶A modification can regulate the function of lncRNAs through the binding of reader proteins, which can affect the progression of tumor formation.

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 m⁶A modification in these three transcripts were further verified by MeRIP-qPCR. MeRIP-qPCR results showed that the relative level of m⁶A modification in HDAC9, MINDY4B was significantly increased in ALV-J induced tumor livers, while the relative level of m⁶A 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 m⁶A mediates dynamic gene expression of ALV-J induced tumor formation.