Chang Fan
Yanzhen Ma1,2
Hui Jiang- 1Experimental Center of Clinical Research, The First Affiliated Hospital of Anhui University of Chinese Medicine, Hefei, China
- 2School of Pharmacy, Anhui University of Chinese Medicine, Hefei, China
- 3Key Laboratory of Xin’an Medicine of the Ministry of Education, Anhui University of Chinese Medicine, Hefei, China
- 4Department of Pharmacy, The First Affiliated Hospital of Anhui University of Chinese Medicine, Hefei, China
- 5Department of Pharmacy, The First Affiliated Hospital of USTC, Division of Life Sciences and Medicine, University of Science and Technology of China, Hefei, China
N6-Methyladenosine (m6A), a unique and common mRNA modification method in eukaryotes, is involved in the occurrence and development of many diseases. Liver fibrosis (LF) is a common response to chronic liver injury and may lead to cirrhosis and even liver cancer. However, the involvement of m6A methylation in the development of LF is still unknown. In this study, we performed a systematic evaluation of hepatic genome-wide m6A modification and mRNA expression by m6A-seq and RNA-seq using LF mice. There were 3,315 genes with significant differential m6A levels, of which 2,498 were hypermethylated and 817 hypomethylated. GO and KEGG analyses illustrated that differentially expressed m6A genes were closely correlated with processes such as the endoplasmic reticulum stress response, PPAR signaling pathway and TGF-β signaling pathway. Moreover, a total of 90 genes had both a significant change in the m6A level and mRNA expression shown by joint analysis of m6A-seq and RNA-seq. Hence, the critical elements of m6A modification, including methyltransferase WTAP, demethylases ALKBH5 and binding proteins YTHDF1 were confirmed by RT-qPCR and Western blot. In an additional cell experiment, we also observed that the decreased expression of WTAP induced the development of LF as a result of promoting hepatic stellate cell (HSC) activation. Therefore, this study revealed unique differential m6A methylation patterns in LF mice and suggested that m6A methylation was associated with the occurrence and course of LF to some extent.
Introduction
N6-Methyladenosine (M6A) is a posttranscriptional modification found in eukaryotic messenger RNA (mRNA), which is similar to DNA methylation and histone modification and is regulated by a variety of methyltransferases (Bushkin et al., 2019; Gu et al., 2019; Berulava et al., 2020). Methyltransferase complexes are composed of METTL3 (methyltransferase-like 3), METTL14 and their additional linker molecules such as WTAP (Wilms tumor associated protein) and KIAA1429, which can catalyze mRNA m6A methylation. The m6A methylation site on RNA is recognized by m6A-binding proteins, including YTHDC1/2 (1ap2 containing YTH domain), YTHDF1/2/3 (YTH family proteins 1–2–3) and IGF2BP1/2/3 (insulin-like growth factor 2 mRNA binding protein 1/2/3), which can bind to methylated m6A sites and perform specific functions. In addition, demethyltransferase FTO (fat mass and obesity related protein) and ALKBH5 (alkyl B homolog 5) reduce m6A modified RNA to original RNA (Du et al., 2018; Zhang Z. et al., 2020; Mapperley et al., 2021). The combined action of these methyltransferases makes m6A modification a dynamic and reversible process (Lu et al., 2020). It has been confirmed that m6A modification affects the control of key cellular processes, including RNA stability (Wang et al., 2014), translation efficiency (Wang et al., 2015), secondary structure (Liu et al., 2015), subcellular localization (Meyer and Jaffrey, 2014), splicing and transport (Yang et al., 2018), and plays important roles in a variety of diseases (Zhang B. et al., 2020; Liu et al., 2020).
Liver fibrosis (LF) is defined as excessive deposition of extracellular matrix (ECM) in response to various cases of liver injury, which is a reversible abnormal tissue response, and excessive activation of hepatic stellate cells (HSCs) is central to its pathogenesis (Bataller and Brenner, 2005; Zhang et al., 2017; Smith-Cortinez et al., 2020). LF is the most common pathological consequence of liver diseases and may lead to liver cirrhosis and liver cancer, and even develop into liver failure in severe cases (Wang Q. et al., 2020). Existing studies have found that m6A methylation plays an extremely important role in a variety of physiological and pathological processes of the liver (Lin et al., 2020; Ondo et al., 2021). Zhong et al. (2019) found that the m6A binding protein YTHDF2 can inhibit tumor proliferation and growth by reducing the stability of EGFR mRNA in hepatocellular carcinoma. Ma et al. (2017) found that the methyltransferase METTL14 can inhibit the metastasis of hepatocellular carcinoma by regulating the methylation of microRNAs. However, as a preliminary process in these severe liver diseases, m6A methylation in LF is rarely described.
The purpose of this study was to establish the expression profile of m6A modification in mice with LF and to explore the potential regulatory mechanism of m6A methylation on LF. Therefore, we used m6A-seq and RNA-seq, to analyze the difference in gene methylation modification and mRNA expression levels after LF at the full transcriptional level, and verified the change in methylase expression and its regulatory role in LF (Figure 1). In conclusion, this study revealed that RNA m6A methylation can play a key role in the pathogenesis of LF by regulating the mRNA expression level of related transcripts. Moreover, methylase affects the occurrence and development of LF by regulating the process of m6A methylation, which could represent an important factor in the process of LF.
FIGURE 1. A schematic diagram of m6A-seq and RNA-seq analyses of mice with LF. LF was induced by subcutaneous injection of CCl4 in mice, and extracted total RNA from liver. Then, RNA was fragmented, and the m6A RNA was separated by immunoprecipitation magnetic beads specifically recognizing m6A sites. Subsequently, the m6A-seq and RNA-seq library were constructed and sequenced.
Materials and Methods
Animal
SPF male C57 BL/6 mice (6–8 weeks old, 20 ± 2 g) were purchased from the Experimental Animal Center of Anhui Province. All mice were raised in the animal facility of the First Affiliated Hospital of Anhui University of Chinese Medicine with an indoor temperature of 18–22°C and humidity of 40–60%, under 12 h alternate dark/light cycles. All mice were allowed food and water freely. Following 1 week of adaptive feeding, a model of LF was established by subcutaneous injection of 0.01 ml/g 20% carbon tetrachloride (CCl4) in an olive oil solution in the back flank of the mice twice a week for 12 weeks, as described in our previous study (Fan et al., 2020). The number of samples was three per group for control mice and LF model mice. The experimental design was approved by the Animal Ethics Committee of Anhui University of Chinese Medicine (AHUCM-mouse-2020032).
Histopathological Analysis
Twelve weeks after establishing the model, the mice were sacrificed by cervical dislocation and the liver samples were taken for histopathological analysis under white light, and hematoxylin and eosin and Masson staining.
Another part of the fresh liver sample was fixed in 2.5% glutaraldehyde and incubated overnight at 4°C. The sample was then fixed in 2% osmium tetroxide for 1 h and dehydrated to 100% through a fractionated series of ethanol (Jiang et al., 2018). Then the sample was embedded in the resin and observed under an electron microscope.
M6A Sequencing and RNA Sequencing
Total RNA was isolated from mouse liver tissue using TRIzol reagent (Invitrogen, United Statesa) according to the manufacturer’s protocol. In this study, we used an m6A-specific antibody (Sigma-Aldrich, ABE572) for immunoprecipitation RNA. The m6A RNA-seq service was provided by Shanghai Bohao Biotechnology Corporation (Shanghai, China). Briefly, poly (A) RNA was captured by VAHTS 2X Frag/Prime Buffer. Then one part of the RNA fragment was used to construct the RNA-seq library, and the other part was used for m6A RNA immunoprecipitation through the GenSeqTM m6A-MERIP kit (GenSeq Inc., Cyberjaya, Malaysia), which was used to construct the m6A-seq library. All operations were carried out in accordance with the manufacturer’s instructions. RNA input samples without immunoprecipitation and m6A input samples were used for the generation of RNA-seq libraries. The library quality was evaluated with a Bioptic Qsep100 Analyzer (Bioptic lnc., Taiwan, China). Library sequencing was performed on an Illumina NovaSeq instrument with 150 bp paired-end reads.
Sequencing Data Processing
Cutadapt (v2.5.0) was used to trim adapters and filter for sequences, FastQC (www.bioinformatics.babraham.ac.uk/projects/fastqc) was used to analyze the quality of sequencing data, and the sequencing mass distribution, base content distribution and repeated sequencing fragment proportion were obtained (Garsmeur et al., 2018). Then, the remaining reads were aligned to the human ensemble genome GRCh38 (mouse ensemble genome GRCm38) using Hisat2 aligner (v2.1.0) under the following parameters: -rna-strandness RF. m6A peaks were identified using the exomePeak R package (v2.13.2) under the following parameters: “PEAK_CUTOFF_PVALUE = 0.05, PEAK_CUTOFF_FDR = NA, FRAGMENT_LENGTH = 200”. Identified m6A peaks with a p value < 0.05 were chosen for the de novo motif analysis using homer (v4.10.4) under the following parameters: “-len 6 -rna”. M6A-RNA-related genomic features were visualized using the Guitar R package (v1.16.0). We used the HOMER (http://homer.ucsd.edu/homer/ngs/peakMotifs.html) software to analyze the motifs of the m6A peaks (Heinz et al., 2010). The screening of differential m6A peaks was also carried out by the exomePeak R package, and the filtering threshold was p value <0.05, |fold change| > 2. Moreover, Bam files of sequencing results were visualized using IGV (http://software.broadinstitute.org/software/igv/) (Robinson et al., 2011).
GO and KEGG Analyses
Differential methylated genes and mRNAs screened according to the above filtering threshold p value <0.05, |fold change| > 2 were used for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses (Ashburner et al., 2000). All analyses were performed using the clusterprofile R package (v3.6.0). Then, the top 20 GO terms and pathways were selected for display according to the p value and the degree of enrichment. The figures were generated using OmicShare tools (http://www.omicshare.com/tools).
Protein-Protein Interaction (PPI) Network Analysis
We conducted a joint analysis of genes with differential expression and differential m6A methylation modification and then used the p value and fold change to screen out the genes for PPI analysis. These differentially expressed genes were imported into the STRING database, which contains comprehensive information about interactions between proteins, and was used to determine the interaction relationship between genes (Szklarczyk et al., 2017). The PPI network was constructed based on importing the data into Cytoscape 3.5.1 software, and then, the network was analyzed by Network Analyzer. The genes with interactions with combined scores greater than 0.4 were selected to construct a protein-protein interaction network diagram (Wang X. et al., 2020).
Validation Experiment
RNA m6A Dot Blot Analysis
A dot blot assay was performed to compare the difference in total m6A levels in liver samples between the control group and the model group. According to the manufacturer’s instructions, the total RNA, was isolated from the liver sample with TRIzol (Thermo, 15596018) and the RNA sample was placed on the nitrocellulose filter membrane. The membrane was dried and cross-linked with 200,000 μJ/cm2 UV twice, washed 3 times with PBST for 5 min each time, and blocked at room temperature for 2 h in 5% skimmed milk. The membrane was transferred to a closed solution containing anti-m6A antibody (ab232905, Abcam) at a dilution of 1: 1,000 and incubated overnight at 4°C. Then, the membrane was rinsed again with PBST for 10 min, sealed in a solution of goat antirabbit IgG combined with HRP (Zs-BIO, ZB-2301) at a dilution of 1: 5,000, incubated at room temperature for 1 h and washed with PBST 3 times. The film was developed with ECL (Western Lightning Plus-ECL, Perkin-Elmer) detection reagent (Thermo, 34094), the signal was detected by chemiluminescence, and the bands were analyzed by ImageJ software.
Isolation and Culture of Primary Mice HSCs
Mice were anesthetized by intraperitoneal injection of pentobarbital sodium and fixed on the operating table. A middle incision of the lower abdomen was used to open the abdominal cavity and exposed the liver and portal vein. Then, the liver was perfused with preheated HBSS at a uniform speed, the open vein was cut when the liver turned gray, and then 0.05% type IV collagenase perfusion solution was perfused (Nishanth et al., 2013; Kim et al., 2016). After perfusion, the liver was cut out and placed in a Petri dish to clean the liver surface with PBS. Tweezers were used to tear up the liver, and 0.05% type IV collagenase was added to the 37°C incubator to digest the tissue for 30 min, followed by filtering with a 200-mesh strainer. The filtrate was centrifuged at 80, 50 and 40 g gradients, and the cell precipitate was collected. The cells were resuspended in serum containing DMEM and seeded in plates precoated with rat-tail collagen I (Zhang et al., 2012; Vig et al., 2015; Yang et al., 2019). After 4 h, the cell culture medium was replaced with serum-free DMEM to continue culturing, and the results of HSC identification are shown in Supplementary Figure 1.
Synthesis and Screening of siRNA and Cell Transfection
To suppress the expression of WTAP, the sequence information of WTAP was obtained from the NCBI database, and the specific WTAP small interference RNA (siRNA) sequence was designed and synthesized according to the full-length sequence information. The specific sequence information is shown in Supplementary Table 1. All siRNA sequences were synthesized by Shanghai Jima Biotechnology Co., Ltd. (Shanghai, China). Three dose groups of 50 pmol, 100 and 200 pmol were set for each siRNA to screen the best transfection conditions. The murine HSCs were seeded in 6-well cell culture plates and cultured until the degree of cell fusion reached 60–80% (Wang Z. et al., 2021). Then, WTAP siRNA was transfected into HSCs with Lipofectamine 2000 transfection reagent (Invitrogen). After 24, 48 and 72 h of siRNA transfection, the HSCs were collected and the expression of WTAP was detected by RT-qPCR assay.
Cell Proliferation Assays and Cell Cycle Analysis
The proliferation of HSCs was detected using a CCK-8 assay. In short, HSCs were trypsinized and resuspended in complete medium, and the cell density was adjusted to 1×105. HSCs were inoculated into 96-well plates at 100 μl per well and cultured for 72 h in a 37°C incubator. Then, 10 μl CCK-8 reagent (BestBio, BB-4202-01) was added to each well, and cells were cultured for another 1 h. The absorbance of each well at 450 nm was measured using a microplate reader. Cell cycle was analyzed by flow cytometry. The HSCs of each group were collected and added to PI staining solution (BestBio, BB-4104) and incubated. The percentage of HSCs in each stage was detected by flow cytometry, and the data were analyzed by FlowJo software (Tree Star Inc., United Statesa).
RT-qPCR
RT-qPCR was used to detect the expression level of candidate genes. Total RNA from HSCs was extracted with TRIzol (Thermo, 15596018). An ultramicro spectrophotometer was used to determine the concentration and purity of RNA. Then, cDNA reverse transcription and RT-qPCR reactions were performed using the PrimeScript™ RT reagent Kit with gDNA Eraser (TaKaRa, RR047A) and 2×SYBR Green qPCR Master Mix (High ROX) (Servicebio, G3322-05). The primer information is shown in Supplementary Table 2. Reactions proceeded using the following conditions: 95°C for 30 s, followed by 40 cycles of 95°C for 15 s and 60°C for 30 s.
Western Blot
Total proteins were obtained from HSCs using the radioimmunoprecipitation assay (RIPA) lysis buffer (Beyotime, P0013B) and PMSF (Biosharp, BL507A). The protein contents of the samples were determined by the bicinchoninic acid method. Twenty micrograms of protein samples were separated by 10% SDS-PAGE and transferred to polyvinylidene fluoride membranes. Following blocking with 5% skim milk for 1 h at room temperature, the membranes were incubated with primary antibodies against WTAP (Affinity, DF3282), YTHDF1 (Affinity, DF3422), ALKBH5 (Affinity, DF2585), α-SMA (Affinity, AF1032), and collagen Ⅰ (Affinity, AF7001) overnight at 4°C. The dilution concentrations of the above antibodies were all 1:1,000. After washing with TBST, diluted goat-anti-mouse IgG (1:10,000) antibody (Zs-BIO, ZB-2305) or goat anti-rabbit IgG (1:3,000) antibody (Zs-BIO, ZB-2301) conjugated with horseradish peroxidase was added, and membranes were incubated for 2 h at room temperature. The membranes were developed with an enhanced chemiluminescence detection kit, and the bands were analyzed by ImageJ software.
Statistical Analysis
The experimental data are presented as the mean ± standard deviation (SD). Statistical analysis was performed by using SPSS 23.0 software. Paired Student’s t-tests were used to detect the differences between the two groups. For multiple comparisons, one-way ANOVA was used with Tukey’s multiple comparisons test. When the p value was <0.05, the results were considered to be statistically significant.
Results
Pathologic HE Staining, Sirius Red Staining and Transmission Electron Microscopy of the Liver
Liver morphology and the pathological changes in LF mice were observed by white light, HE staining, Masson staining and transmission electron microscopy. As shown in Figure 2A, after 12 weeks of CCl4 induction, the livers of the control group were red and smooth, while the livers of the model group were relatively swollen and rough, and the color was gray and white. In Figure 2B, the results of HE staining showed that the structure of the hepatic lobules in the control group was clear, and the hepatocyte cords were in their normal arrangement. In contrast, in the model group there were abundant and large lipid droplets in the cytoplasm of hepatocytes, severe steatosis, disordered liver tissue structure, obvious hyperplasia of fibrotic tissue, and unclear structure of some hepatic lobules.
FIGURE 2. Collected livers were subjected to pathological analysis by white light, HE and Sirius red staining and transmission electron microscopy. (A) Liver under white light. (B) HE staining (200-fold). (C) Sirius Red staining (200-fold). (D) Transmission electron microscopy (TEM) analysis (20000-fold).
In Figure 2C, the results of Masson staining showed that there was a large amount of collagen deposition in the liver tissue of the model group compared with the control group. Similarly, obvious changes in the subcellular structure of the liver were observed under an electron microscope (Figure 2D). Hepatocytes in the control group were intact and without morphological signs of degeneration or necrosis, while in the model group, the hepatocytes showed abnormal morphological changes, including disappearance of the cell boundary, rupture of the cell membrane, cytoplasmic turbidity, organelle expansion and nuclear shrinkage.
General Description of m6A Methylation Modification in LF
We compared m6A methylation peaks at each site in hepatic tissues from mice with fibrosis. The differences and overlaps in m6A methylation between the individuals are shown by the Venn diagram in Figure 3A. We found 6,221 m6A methylation modifier genes in the control group and 6,982 m6A methylation modifier genes in the model group, of which 5,111 m6A methylation modifier genes were common between the two groups. Compared with the control group, 1871 m6A methylation modifier genes appeared, and 1,110 m6A methylation modifier genes disappeared in the model group, indicating that there was a significant difference in the m6A modification pattern after LF. Figure 3B shows the level of m6A methylation in different groups. We found an average of 12166 peaks in the control group and 15100 peaks in the model group.
FIGURE 3. Overview of m6A-modified transcripts in LF mice. (A) Venn diagram of m6A-modified genes in the control group and the model group. (B) The average number of m6A peaks in each group. (C) Density of differential m6A peaks along transcripts. Each transcript was divided into three parts: 5′UTR, CDS, and 3′UTR. (D) Pie charts showing the region of m6A peaks in each group. (E) Violin plot of the relative abundance of m6A peaks in each group. (F) Number of peaks per transcript. (G) The most conserved sequence motif of the differential m6A peak region. (H) The distribution patterns of m6A peaks in different chromosomes. (I) The count of m6A peaks in per chromosome.
As shown in Figures 3C,D, m6A methylation of mRNAs occurred mainly in coding sequences (CDSs) and 3′ untranslated regions (3′UTRs). More specifically, approximately 35.7% of m6A peaks were distributed in the CDS region, and 33% of m6A peaks were distributed in the 3′UTR. The violin diagram (Figure 3E) shows the results of the enrichment degree analysis of m6A methylation in each sample. The average logarithmic fold-enrichment of the control group was 4.8, while the average logarithmic fold-enrichment of the model group was 5.3. By means of the distribution of m6A peaks in each gene, we found that approximately 37% of the genes had separate m6A modification sites, and 80% of the genes had one to three m6A modification sites (Figure 3F).
Subsequently, we predicted the m6A motif in LF by the mRNA sequence corresponding to m6A methylation peaks. As shown in Figure 3G, the most significant mRNA methylation occurred at the RRAC motifs. The analysis of the m6A methylation distribution at different chromosome loci found that the m6A peaks of genes in the model group increased, and the chromosomes with the highest m6A methylation frequency were chromosome 7 with 1,119 m6A methylation peaks, chromosome 11 with 993 m6A methylation peaks and chromosome 2 with 940 m6A methylation peaks (Figures 3H,I). By further comparison, we found that there was no significant difference in the distribution number of m6A peaks on chromosomes between the two groups.
Analysis of Differentially Methylated m6A Genes and Their Signaling Pathways
Using the filtering criteria of a p value <0.05 and |fold change| >2, 3,315 genes with differential m6A methylation were identified, of which 2,498 m6A hypermethylated genes and 817 m6A hypomethylated genes were identified (Figures 4A,B). We also visually assessed the enrichment degree and fold change of the top 10 hypermethylated genes and top 10 hypomethylated genes (Figure 4C), as shown in Table 1. Specific information of all differentially methylated m6A genes is presented in Supplementary file 1.
FIGURE 4. Genes with differential m6A methylation modification in LF. (A) Volcano plot representation of microarray data on the differentially expressed m6A methylation genes. The blue and red dots to the left and to the right of the two vertical lines indicate more than a 2-fold change and represent the differentially expressed m6A methylation genes with statistical significance. (B) Hierarchical cluster analysis of differentially expressed m6A methylation genes. Hierarchical clustering shows that the differentially expressed m6A methylation genes ultimately cluster into two major branches, including hypermethylated genes, which are labeled in red, and hypomethylated genes, which are labeled in green. The darker the color, the more significant the difference. (C) The radar map shows the top 10 most significant hypermethylated genes and top 10 hypomethylated genes. (D) GO biological processes enrichment analysis. (E) GO cellular component enrichment analysis. (F) GO molecular function enrichment analysis. (G) KEGG enrichment analysis.
TABLE 1. the top 10 hypermethylation genes and top 10 hypomethylation genes.
Simultaneously, the results of GO and KEGG analyses showed the enrichment of GO functions and pathways of differentially methylated genes. We found 1122 GO terms were significantly enriched in biological processes (Figure 4D), 210 GO terms were significantly enriched in cellular components (Figure 4E), and 476 GO terms were significantly enriched in molecular functions (Figure 4F), especially in the process of transcription, liver development, response of endoplasmic reticulum to unfolded proteins, and protein binding. Similarly, KEGG analysis found that 104 pathways were significantly enriched (Figure 4G), especially protein processing in the endoplasmic reticulum, PI3K-Akt signaling pathway and TGF-β signaling pathway. Specific information on the GO and KEGG pathway enrichment analyses is presented in Supplementary Table 3.
Description of mRNA Expression and Analysis of Differential Genes in LF
In Figure 5A, not only the mRNA distribution and abundance of control samples and LF samples were shown, but also the peak patterns of these samples were visually displayed. The violin diagram in Figure 5B demonstrates a similar result; the average logarithmic fold-enrichment of the control group was 1.2, while the average logarithmic fold-enrichment of the model group was 1.3. The gene distribution pattern of the control group was also different from the gene distribution pattern of the model group, but they were distributed mainly in the CDS region and exon region (Figure 5C).
FIGURE 5. The overall expression of mRNA and the description of differentially expressed mRNAs. (A) Metagene plots reveal the distribution intensity and abundance of mRNA expression after sequence alignment. (B) Violin plot of the relative abundance of mRNA expression in each sample. (C) Regional distribution of mRNA. (D) Volcano plot representation of microarray data on the differentially expressed mRNA genes. (E) Hierarchical cluster analysis of differentially expressed mRNA genes. (F) The radar map shows the top 10 upregulated genes and top 10 downregulated genes. (G) GO biological processes enrichment analysis. (H) GO cellular component enrichment analysis. (I) GO molecular function enrichment analysis. (J) KEGG enrichment analysis.
Then, similar to the screening of differentially methylated genes, a p value <0.05 and |fold change| > 2 were used as screening criteria, and we found 828 differentially expressed genes, including 398 upregulated genes and 430 downregulated genes (Figures 5D,E). Moreover, we also visually compared the expression and corresponding abundance of the top 10 upregulated genes and top 10 downregulated genes (Figure 5F), as shown in Table 2. Specific information of all differentially expressed RNAs is presented in Supplementary File 2. Meanwhile, the results of GO analysis showed that 376 GO terms were significantly enriched in biological processes (Figure 5G), 64 GO terms were significantly enriched in cellular components (Figure 5H), and 136 GO terms were significantly enriched in molecular functions (Figure 5I), particularly in cellular response to hormone stimulus, proteinaceous extracellular matrix, extracellular matrix structural constituent, and more. Similarly, in Figure 4J, the results of KEGG analysis found that 41 pathways were significantly enriched (Figure 4J), particularly the metabolism of xenobiotics by cytochrome P450, retinol metabolism, chemical carcinogenesis, and more. Specific information on the GO and KEGG pathway enrichment analyses is presented in Supplementary Table 4.
TABLE 2. the top 10 up-regulated genes and top 10 down-regulated genes.
Overview of Transcriptome Profiles and Conjoint Analyses of m6A-Seq and RNA-Seq Data
A conjoint analysis was conducted for m6A-seq and RNA-seq data. We found that a total of 8,299 peaks located on 2,353 genes not only had m6A modification but also had altered mRNA levels (Figure 6A). However, not all of them were significant. As shown in Figure 6B, by setting the filter conditions of a p value < 0.05 and |fold change| >2, we found 90 genes that commonly had significant differential m6A methylation levels and significant differentially expressed mRNA levels. Among these genes, there were 4 genes with m6A hypomethylation and downregulated mRNA expression, 51 genes with m6A hypermethylation and downregulated mRNA expression, 26 genes with m6A hypermethylation and upregulated mRNA expression and 9 genes with m6A hypomethylation and upregulated mRNA expression. The specific information on these genes is shown in Supplementary Table 5.
FIGURE 6. Joint analysis of m6A methylation and mRNA expression. (A) Venn diagram of peaks with m6A methylation and mRNA. (B) Four quadrant graph of genes with differential m6A methylation and differentially expressed mRNA levels. (C) Cumulative frequency plot showing that there was a correlation between differential m6A methylation genes and mRNA levels. (D) PPI of genes with differentially expressed m6A methylation and differentially expressed mRNA. (E) GO biological processes enrichment analysis. (F) GO cellular component enrichment analysis. (G) GO molecular function enrichment analysis. (H) KEGG enrichment analysis.
Subsequently, we confirmed the correlation between m6A modification and mRNA levels. The results in Figure 6C show that differential m6A-methylated transcripts do have different mRNA expression levels; that is, the mRNA expression level of hypomethylated transcripts is often higher than the mRNA expression level of hypermethylated transcripts. Based on interactions with combined scores ≥0.4, the PPI network analysis constructed interaction networks for these differential genes, as shown in Figure 6D.
The results of GO analysis showed that 670 GO terms were significantly enriched in biological processes (Figure 6E), 85 GO terms were significantly enriched in cellular components (Figure 6F), and 148 GO terms were significantly enriched in molecular functions (Figure 6G), particularly in lipid biosynthetic process, endoplasmic reticulum correlation, structural constituent of cytoskeleton, and more. Similarly, in Figure 6H, the results of KEGG analysis found that 29 pathways were significantly enriched, particularly steroid hormone biosynthesis, chemical carcinogenesis, gap junction, and more. The specific information of GO and KEGG pathway enrichment analyses is presented in Supplementary Table 6.
Levels of m6A Methylation and Methylase Expression in LF
To further explore the changes in m6A methylation in LF, we performed an m6A dot blot analysis. The results showed that compared with the control group, the m6A methylation abundance of the model group was significantly decreased (Figures 7A,B). Subsequently, considering that the difference in m6A levels in LF was probably caused by m6A regulatory enzymes, we focused on the methyltransferase WTAP, demethylase ALKBH5 and m6A binding protein YTHDF1. IGV visualization analysis was used to show the sequencing results intuitively. At the m6A methylation level, we found that the m6A levels of WTAP and ALKBH5 increased, while the YTHDF1 level decreased in LF (Figures 7C–E).
FIGURE 7. Verification of m6A methylation level and methylase expression in LF. (A) The m6A methylation level in LF. (B) Semiquantitative analysis of m6A methylation. (C) IGV plots of the WTAP m6A level. (D) IGV plots of the ALKBH5 m6A level. (E) IGV plots of the YTHDF1 m6A level. (F) IGV plots of the WTAP expression level. (G) IGV plots of the ALKBH5 expression level. (H) IGV plots of the YTHDF1 expression level. (I) mRNA expression level of WTAP. (J) mRNA expression level of ALKBH5. (K) mRNA expression level of YTHDF1. (L) Protein expression levels of WTAP. (M) Protein expression levels of ALKBH5. (N) Protein expression levels of YTHDF1. (O) Semiquantitative analysis of WTAP protein. (P) Semiquantitative analysis of YTHDF1 protein. (Q) Semiquantitative analysis of ALKBH5 protein. ##p < 0.01 compared with the control group, #p < 0.05 compared with the control group.
Likewise, at the mRNA level, we found that the expression of WTAP, ALKBH5 and YTHDF1 was reduced (Figures 7F–H) by IGV visualization analysis. Then, an RT-qPCR assay was utilized to examine the expression of the above genes. The results showed that the expression levels of WTAP, ALKBH5 and YTHDF1 in the model group were significantly lower than those in the control group, which was consistent with the IGV results (Figures 7I–K). Moreover, we also verified the protein expression levels of WTAP, ALKBH5 and YTHDF1 by Western blot and found that the protein levels of the three genes also decreased significantly in the model group (Figure 7L–Q).
Effects of Methyltransferase WTAP on Proliferation, Cell Cycle and Activation Markers of HSCs
As shown in Figure 8A, we analyzed the expression of WTAP in human LF samples through the GEO database (GSE33650) and found that the expression level of WTAP in high-fibrosis samples was significantly lower than the expression level of WTAP in low-fibrosis samples, which was consistent with our present experimental results. Furthermore, we designed and synthesized small interfering RNA targeting WTAP. As shown in Figures 8B–D, we screened the small interfering RNA sequences, durations and concentrations of WTAP small interfering RNA using RT-qPCR and found that the optimal interference sequence was si-WTAP-1, the optimum time of siRNA treatment for interference was 48 h, and the optimum concentration of siRNA for interference was 100 pmol. Follow-up experiments were carried out according to the above conditions.
FIGURE 8. Effects of methyltransferase WTAP on proliferation, cell cycle and activation markers of HSCs. (A) Expression levels of WTAP in low-fibrosis and high-fibrosis samples derived from the GEO database. (B) Small interfering RNA of WTAP was screened by RT-qPCR assay. (C) Optimal stimulation time of WTAP small interfering RNA was screened by RT-qPCR assay. (D) Optimal stimulation concentration of WTAP small interference RNA was screened by RT-qPCR assay. (E) Cell proliferation was detected by CCK8 assay. (F) The phase of the cell cycle was detected by flow cytometry. a, control group. b, model group. c, si-WTAP group. d, si-NC group. (G) Quantification of the cell cycle results. (H) mRNA expression level of α-SMA. (I) The mRNA expression level of collagen Ⅰ. (J) Protein expression levels of α-SMA. (K) Protein expression levels of collagen Ⅰ. (L) Semiquantitative analysis of α-SMA protein. (M) Semiquantitative analysis of collagen Ⅰ protein.
As shown in Figure 8E, the CCK-8 assay results showed that compared with the control group, the proliferation of HSCs in the model group increased, while the proliferation of HSCs further increased after interfering with the expression of WTAP. Then, flow cytometry was used to detect differences in the HSC cell cycle under WTAP interference (Figures 8F,G). The results showed that the number of HSCs in the G0/G1 phase in the model group was significantly lower than that in the control group, while the number of HSCs in S phase and G2/M phase increased significantly. Compared with the model group, the number of HSCs in G0/G1 phase in the si-WTAP group further decreased, while the number of HSCs in the S phase and G2/M phase further increased. Interfering with WTAP promotes the proliferation of HSCs by inducing S phase and G2/M phase arrest.
Moreover, we also detected the expression of the HSC activation markers α-SMA and collagen Ⅰ. As shown in Figure 8H-8M, the mRNA and protein expression levels of α-SMA and collagen Ⅰ were significantly increased in the model group, while the mRNA and protein expression levels of α-SMA and collagen I were further increased after WTAP interference compared with expression in the model group, which also indicated that WTAP interference significantly promoted the activation of HSCs.
Discussion
Modifications through m6A methylation modification, as a kind of RNA modification that exists widely in liver disease, has naturally received extensive attention (Wu et al., 2020; Pan et al., 2021). With regard to the effect of m6A methylation on the biological function of liver cells, existing studies have focused on the regulatory mechanism of genes and pathways (Zhang C. et al., 2020; Cao et al., 2021). A study by Zhu Y. et al. (2020) found that METTL3-mediated m6A methylation could be regulated by ASIC1a, which in turn affects the processing of miR-350, thus inducing the activation of HSCs and promoting the occurrence and development of LF. Unlike their studies, our study compared the difference in m6A methylation between the control and LF liver tissue, and confirmed that the m6A modification level changed significantly in LF.
Herein, we first constructed m6A-seq and RNA-seq libraries and investigated the changes in m6A methylation and the expression levels of genes in the liver of mice with hepatic fibrosis by methylated RNA immunoprecipitation combined with next-generation sequencing, and the results were analyzed by bioinformatics. We found 6,221 m6A modification genes in the control group and 6,982 m6A modification genes in the model group. Further analysis identified 3,315 different m6A methylation genes, of which 2,498 m6A hypermethylated genes and 817 m6A hypomethylated genes were identified, suggesting that there are some differences in the occurrence and development of m6A methylation in LF. Interestingly, although the m6A methylation of the gene was different, the distribution of m6A methylation in the control livers was similar to that in the model livers. We found that m6A methylation of most genes was distributed in CDS, 3′UTR and stop codon regions, accounting for 90% of the total. This is consistent with the report of Dominissini et al. (2012), who found that m6A methylation sites are mainly concentrated in long exons, stop codons and 3′UTR regions, and this distribution pattern is highly conserved between humans and mice. This distribution pattern may be related to the function of m6A methylation modification. Dynamic m6A modification in different regions affects biological functions such as splicing, output, stability and translation of mRNA (Wang et al., 2014; Wang and He, 2014; Maity and Das, 2016). Therefore, m6A modification may play an important role in regulating the expression of genes related to hepatic fibrosis.
The m6A methylation site exists mainly in the RRACH motif, which is caused by the binding of m6A methyltransferase with the corresponding consensus sequence (Liu et al., 2018; Zhang Z. et al., 2019). The RNA binding motifs of METTL3, METTL14 and WTAP have been confirmed to be GGAC, GGAC and GACU, which are highly conserved between humans and animals (Liu et al., 2014). When the RRACH sequence is mutated, the single nucleotide polymorphism of the corresponding site changes, which affects m6A methylation. Kane et al. (Kane and Beemon, 1987) found that the mutation from GAC to GAU in the consensus sequence leads to the reversal of m6A methylation in Rous sarcoma virus mRNA transcripts. In the current study, we found many similar m6A consensus motifs in the control and LF tissues, but there were also some differences in the sequences, which further confirmed the emergence of specific m6A methylation sites in the process of LF. However, the RRACH consensus sequence is critical for m6A methylation, but not all RRACH sites in the body will have m6A modification (Gilbert et al., 2016), which corresponds to our results; that is, there are unmutated sequence sites, showing that m6A methylation modification is also regulated by other molecular mechanisms and needs further study.
To better understand the functions of these differentially expressed m6A methylated genes, GO and KEGG distribution analyses were conducted. We found that differential m6A genes were primarily involved in biological processes associated with the endoplasmic reticulum stress response, such as the unfolded protein response and the protein catabolic process, and were also related to the development and regeneration of liver organs. In addition, they were closely related to the PPAR signaling pathway, TGF-β signaling pathway and PI3K-Akt signaling pathway. Endoplasmic reticulum stress refers to the state of protein folding damage caused by the destruction of endoplasmic reticulum homeostasis, and some studies have confirmed that endoplasmic reticulum stress plays a role in the occurrence and development of various liver diseases (Huang et al., 2019; Wu et al., 2021). Virginia et al. (Hernández-Gea et al., 2013) found that oxidative stress disrupts endoplasmic reticulum homeostasis in stellate cells and causes the endoplasmic reticulum to enter a stressed state. To reduce the stress response, hepatic stellate cells initiate an unfolded protein response by limiting the accumulation of unfolded proteins during transient stress, which promotes cell activation and accelerates the development of LF. Peroxisome proliferation-activated receptor (PPAR) belongs to the nuclear hormone receptor family and plays an important role in many biological processes, such as adipogenesis (Lefterova et al., 2014), cell differentiation (Kim et al., 2019), cell growth regulation (Zhang X. et al., 2019) and inflammation (Bougarne et al., 2018). Previous studies have found that the activation of the PPAR pathway can delay the progression of hepatic fibrosis, and its activation can inhibit the transformation of HSCs from a resting state to an activated state (Guo et al., 2005; Anty and Lemoine, 2011). Liu and others have further found that the activation of PPAR-γ can reduce the expression of α-SMA and collagen I in HSCs (Yang et al., 2006). Both the TGF-β and PI3K-Akt signaling pathways are one of the classical signaling pathways involved in the progression of LF. Abnormalities in TGF-β can stimulate HSCs to secrete excessive ECM, and the activity of the PI3K-Akt signaling pathway is significantly correlated with collagen production, HSC proliferation and apoptosis (Shah et al., 2013; Wu et al., 2017). Interestingly, the fibrogenic effects of TGF-β and PI3K-Akt are synergistic to some extent. Runyan et al. (2004) found that TGF-β can not only induce the activation of PI3K/Akt, but also enhance the transcriptional activity of Smad3, the target downstream of TGF-β signaling, thus enhancing the expression of collagen I.
By combining analyses of m6A-seq and RNA-seq data, we discovered 90 genes with differences in their m6A methylation peaks and synchronously differential mRNA expression in LF. The expression of these genes may be regulated by m6A modification of mRNAs. Among the genes with the highest differences, many have been identified to be closely related to the occurrence and development of LF, such as ApoA4 (apolipoprotein A4). Wang Y. et al. (2021) found that ApoA4 may reduce LF and liver injury by inhibiting LF mediators and inflammatory cytokines and suppressing proinflammatory hepatic M1 cell invasion. Although some genes have not been proven to be related to LF, they are involved in fibrosis in other tissues. For example, Ninj1 has been shown to promote the activation of macrophages by enhancing the interaction with epithelial cells, thus enhancing the inflammatory response of macrophages to participate in the occurrence and development of pulmonary fibrosis (Choi et al., 2018). These genes regulated by m6A modification may play key roles in the occurrence and development of LF and may also become an important target for the treatment of LF. However, the specific molecular mechanism of the effect of m6A methylation of these genes on LF is still unclear and needs further exploration and research in the future.
The most prominent finding in our data is that there is a significant difference in m6A modification between the LF and control tissues. The dot blot results also confirmed this significant difference, and we found that the overall level of m6A methylation in LF decreased significantly, which suggested that the modification of the m6A genes affected the progression of LF. A possible explanation for the global change in this m6A modification pattern may be the unique expression of the key m6A regulator or its own methylation modification. Considering that methylases play very important roles in regulating m6A methylation of liver fibrosis, we selected WTAP, ALKBH5 and YTHDF1 as the representative of methyltransferase, demethylase and m6A binding protein for further study, which verify the differences in mRNA and protein expression levels. Interestingly, not only did the expression of WTAP and YTHDF1 decrease in LF, but the expression of the demethylase ALKBH5 also decreased significantly. Combined with the decrease in the overall level of m6A modification in LF, we speculated that the m6A level in the body involves the regulation of a variety of methylases, and the change in one or several methylation enzymes alone cannot be used as a decisive factor in determining the level of m6A methylation. The decrease in the m6A level in LF was because the overall degree of demethylation was greater than the decrease in the m6A level of methylation.
As an important component of the m6A methyltransferase complex, WTAP, unlike METTL3 and METTL14, does not have N6-methyladenine methyltransferase activity but is necessary for effective RNA methylation in vivo and for the localization of METTL3 and METTL14 in nuclear spots (Śledź and Jinek, 2016). WTAP has been proven to participate in some basic physiological processes, such as mRNA stability (Horiuchi et al., 2006), organ development (Anderson et al., 2014), cell proliferation, apoptosis and cell cycle regulation (Horiuchi et al., 2013). A recent study by Zhu B. et al. (2020) demonstrated that in a rat model of balloon injury-induced hyperplasia of vascular smooth muscle cells (VSMCs), the expression of WTAP decreased significantly. The suggested mechanism is that WTAP regulates p16INK4a through m6A modification and thus causing abnormal proliferation of VSMCs. Nevertheless, contrary to the above findings that WTAP can inhibit cell proliferation, some other studies have shown different results. A study by Chen et al. (2020) confirmed that WTAP could regulate the stability of HMBOX1 mRNA in an m6A methylation-dependent manner, thereby promoting the proliferation and metastasis of osteosarcoma cells. These studies confirmed that as a pivotal enzyme of m6A modification, WTAP can regulate the m6A methylation level in the body, thus fulfilling functionally different roles in different diseases.
Interestingly, in the present study, we found through sequencing that the m6A level of WTAP was significantly upregulated in LF mice, while the expression of mRNA was reduced. Further verification experiments showed that the mRNA and protein expression levels of WTAP decreased significantly, consistent with the sequencing results. Subsequently, we focused on the effect of WTAP interference on HSCs in LF and found that interfering with WTAP promoted the proliferation of HSCs and increased the expression of α-SMA, a marker of HSC activation and collagen I, the main component of extracellular matrix, which indicated that interfering with WTAP could promote the occurrence and development of LF. Therefore, based on the findings of the above study, we speculated that the possible mechanism of WTAP involved in the development of LF was that WTAP acted as a methyltransferase to affect the m6A level on downstream target genes related to cell proliferation and the cell cycle, thus regulating the mRNA expression levels of these genes and ultimately affecting the occurrence and development of LF. These findings may provide new thoughts and insights for other research on WTAP and m6A methylation in LF.
In summary, our findings established a m6A transcriptome map of LF mice, provided a comprehensive investigation of the potential relationship between m6A methylation and mRNA expression in LF, and revealed the key enzymes of m6A modification, especially WTAP, involved in the occurrence and development of LF.
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: BioProject: PRJNA761579, SRA accession: SRP336482.
Ethics Statement
The animal study was reviewed and approved by the Animal Ethics Committee of Anhui University of Chinese Medicine.
Author Contributions
HJ made substantial contributions to the conception and design of the study. CF, YM, SC, and QZ performed the experiments. CF, HJ, JZ, and FW contributed to data acquisition, data analysis and interpretation. HJ revised the manuscript critically for important intellectual content. HJ, JZ and FW confirm the authenticity of all the raw data. All authors agreed to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of the work are appropriately investigated and resolved. All authors read and approved the final manuscript.
Funding
National Natural Science Foundation of China (grant no.81973648). National Natural Science Foundation of China (grant no. 82004139). Leading Talents Introduction and Cultivation Plan Project of Colleges in Anhui Province (gxfxZD2016118).
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 are grateful to Shanghai Bohao Biotechnology Corporation (Shanghai, China) for providing sequencing service. We thank AJE (www.aje.cn) for its linguistic assistance during the preparation of this manuscript.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.767051/full#supplementary-material
References
Anderson, A. M., Weasner, B. P., Weasner, B. M., and Kumar, J. P. (2014). The Drosophila Wilms׳ Tumor 1-Associating Protein (WTAP) Homolog Is Required for Eye Development. Develop. Biol. 390 (2), 170–180. doi:10.1016/j.ydbio.2014.03.012
Anty, R., and Lemoine, M. (2011). Liver Fibrogenesis and Metabolic Factors. Clin. Res. Hepatol. Gastroenterol. 35 (Suppl. 1), S10–S20. doi:10.1016/s2210-7401(11)70003-1
Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., et al. (2000). Gene Ontology: Tool for the Unification of Biology. Nat. Genet. 25 (1), 25–29. doi:10.1038/75556
Bataller, R., and Brenner, D. A. (2005). Liver Fibrosis. J. Clin. Invest. 115 (2), 209–218. doi:10.1172/jci24282
Berulava, T., Buchholz, E., Elerdashvili, V., Pena, T., Islam, M. R., Lbik, D., et al. (2020). Changes in m6A RNA Methylation Contribute to Heart Failure Progression by Modulating Translation. Eur. J. Heart Fail. 22 (1), 54–66. doi:10.1002/ejhf.1672
Bougarne, N., Weyers, B., Desmet, S. J., Deckers, J., Ray, D. W., Staels, B., et al. (2018). Molecular Actions of PPARα in Lipid Metabolism and Inflammation. Endocr. Rev. 39 (5), 760–802. doi:10.1210/er.2018-00064
Bushkin, G. G., Pincus, D., Morgan, J. T., Richardson, K., Lewis, C., Chan, S. H., et al. (2019). m6A Modification of a 3′ UTR Site Reduces RME1 mRNA Levels to Promote Meiosis. Nat. Commun. 10 (1), 3414. doi:10.1038/s41467-019-11232-7
Cao, X., Shu, Y., Chen, Y., Xu, Q., Guo, G., Wu, Z., et al. (2021). Mettl14-Mediated m6A Modification Facilitates Liver Regeneration by Maintaining Endoplasmic Reticulum Homeostasis. Cell Mol. Gastroenterol. Hepatol. 12 (2), 633–651. doi:10.1016/j.jcmgh.2021.04.001
Chen, S., Li, Y., Zhi, S., Ding, Z., Wang, W., Peng, Y., et al. (2020). WTAP Promotes Osteosarcoma Tumorigenesis by Repressing HMBOX1 Expression in an m6A-dependent Manner. Cell Death Dis 11 (8), 659. doi:10.1038/s41419-020-02847-6
Choi, S., Woo, J. K., Jang, Y.-S., Kang, J.-H., Hwang, J.-I., Seong, J. K., et al. (2018). Ninjurin1 Plays a Crucial Role in Pulmonary Fibrosis by Promoting Interaction between Macrophages and Alveolar Epithelial Cells. Sci. Rep. 8 (1), 17542. doi:10.1038/s41598-018-35997-x
Dominissini, D., Moshitch-Moshkovitz, S., Schwartz, S., Salmon-Divon, M., Ungar, L., Osenberg, S., et al. (2012). Topology of the Human and Mouse m6A RNA Methylomes Revealed by m6A-Seq. Nature 485 (7397), 201–206. doi:10.1038/nature11112
Du, Y., Hou, G., Zhang, H., Dou, J., He, J., Guo, Y., et al. (2018). SUMOylation of the m6A-RNA Methyltransferase METTL3 Modulates its Function. Nucleic Acids Res. 46 (10), 5195–5208. doi:10.1093/nar/gky156
Fan, C., Wu, F. R., Zhang, J. F., and Jiang, H. (2020). A Network Pharmacology Approach to Explore the Mechanisms of Shugan Jianpi Formula in Liver Fibrosis. Evidence-Based Complement. Altern. Med. 2020, 4780383. doi:10.1155/2020/4780383
Garsmeur, O., Droc, G., Antonise, R., Grimwood, J., Potier, B., Aitken, K., et al. (2018). A Mosaic Monoploid Reference Sequence for the Highly Complex Genome of Sugarcane. Nat. Commun. 9 (1), 2638. doi:10.1038/s41467-018-05051-5
Gilbert, W. V., Bell, T. A., and Schaening, C. (2016). Messenger RNA Modifications: Form, Distribution, and Function. Science 352 (6292), 1408–1412. doi:10.1126/science.aad8711
Gu, C., Wang, Z., Zhou, N., Li, G., Kou, Y., Luo, Y., et al. (2019). Mettl14 Inhibits Bladder TIC Self-Renewal and Bladder Tumorigenesis through N6-Methyladenosine of Notch1. Mol. Cancer 18 (1), 168. doi:10.1186/s12943-019-1084-1
Guo, Y.-T., Leng, X. S., Li, T., Peng, J. R., Song, S. H., Xiong, L. F., et al. (2005). Effect of Ligand of Peroxisome Proliferator-Activated Receptor γ on the Biological Characters of Hepatic Stellate Cells. World J. Gastroenterol. 11 (30), 4735–4739. doi:10.3748/wjg.v11.i30.4735
Heinz, S., Benner, C., Spann, N., Bertolino, E., Lin, Y. C., Laslo, P., et al. (2010). Simple Combinations of Lineage-Determining Transcription Factors Prime Cis-Regulatory Elements Required for Macrophage and B Cell Identities. Mol. Cel 38 (4), 576–589. doi:10.1016/j.molcel.2010.05.004
Hernández-Gea, V., Hilscher, M., Rozenfeld, R., Lim, M. P., Nieto, N., Werner, S., et al. (2013). Endoplasmic Reticulum Stress Induces Fibrogenic Activity in Hepatic Stellate Cells through Autophagy. J. Hepatol. 59 (1), 98–104. doi:10.1016/j.jhep.2013.02.016
Horiuchi, K., Umetani, M., Minami, T., Okayama, H., Takada, S., Yamamoto, M., et al. (2006). Wilms' Tumor 1-Associating Protein Regulates G2/M Transition through Stabilization of Cyclin A2 mRNA. Proc. Natl. Acad. Sci. 103 (46), 17278–17283. doi:10.1073/pnas.0608357103
Horiuchi, K., Kawamura, T., Iwanari, H., Ohashi, R., Naito, M., Kodama, T., et al. (2013). Identification of Wilms' Tumor 1-Associating Protein Complex and its Role in Alternative Splicing and the Cell Cycle. J. Biol. Chem. 288 (46), 33292–33302. doi:10.1074/jbc.M113.500397
Huang, H.-H., Lee, W.-J., Chen, S.-C., Chen, T.-F., Lee, S.-D., and Chen, C.-Y. (2019). Bile Acid and Fibroblast Growth Factor 19 Regulation in Obese Diabetics, and Non-Alcoholic Fatty Liver Disease after Sleeve Gastrectomy. J. Clin. Med. 8 (6), 815. doi:10.3390/jcm8060815
Jiang, X., Jiang, L., Shan, A., Su, Y., Cheng, Y., Song, D., et al. (2018). Targeting Hepatic miR-221/222 for Therapeutic Intervention of Nonalcoholic Steatohepatitis in Mice. EBioMedicine 37, 307–321. doi:10.1016/j.ebiom.2018.09.051
Kane, S. E., and Beemon, K. (1987). Inhibition of Methylation at Two Internal N6-Methyladenosine Sites Caused by GAC to GAU Mutations. J. Biol. Chem. 262 (7), 3422–3427. doi:10.1016/s0021-9258(18)61520-0
Kim, R. S., Hasegawa, D., Goossens, N., Tsuchida, T., Athwal, V., Sun, X., et al. (2016). The XBP1 Arm of the Unfolded Protein Response Induces Fibrogenic Activity in Hepatic Stellate Cells through Autophagy. Sci. Rep. 6, 39342. doi:10.1038/srep39342
Kim, Y.-H., Jang, W.-G., Oh, S.-H., Kim, J.-W., Lee, M. N., Song, J. H., et al. (2019). Fenofibrate Induces PPARα and BMP2 Expression to Stimulate Osteoblast Differentiation. Biochem. Biophysical Res. Commun. 520 (2), 459–465. doi:10.1016/j.bbrc.2019.10.048
Lefterova, M. I., Haakonsson, A. K., Lazar, M. A., and Mandrup, S. (2014). PPARγ and the Global Map of Adipogenesis and beyond. Trends Endocrinol. Metab. 25 (6), 293–302. doi:10.1016/j.tem.2014.04.001
Lin, Z., Niu, Y., Wan, A., Chen, D., Liang, H., Chen, X., et al. (2020). RNA M 6 A Methylation Regulates Sorafenib Resistance in Liver Cancer through FOXO 3‐mediated Autophagy. Embo j 39 (12), e103181. doi:10.15252/embj.2019103181
Liu, J., Yue, Y., Han, D., Wang, X., Fu, Y., Zhang, L., et al. (2014). A METTL3-METTL14 Complex Mediates Mammalian Nuclear RNA N6-Adenosine Methylation. Nat. Chem. Biol. 10 (2), 93–95. doi:10.1038/nchembio.1432
Liu, N., Dai, Q., Zheng, G., He, C., Parisien, M., and Pan, T. (2015). N6-Methyladenosine-Dependent RNA Structural Switches Regulate RNA-Protein Interactions. Nature 518 (7540), 560–564. doi:10.1038/nature14234
Liu, B., Merriman, D. K., Choi, S. H., Schumacher, M. A., Plangger, R., Kreutz, C., et al. (2018). A Potentially Abundant Junctional RNA Motif Stabilized by m6A and Mg2+. Nat. Commun. 9 (1), 2761. doi:10.1038/s41467-018-05243-z
Liu, T., Wei, Q., Jin, J., Luo, Q., Liu, Y., Yang, Y., et al. (2020). The m6A Reader YTHDF1 Promotes Ovarian Cancer Progression via Augmenting EIF3C Translation. Nucleic Acids Res. 48 (7), 3816–3831. doi:10.1093/nar/gkaa048
Lu, J., Qian, J., Yin, S., Zhou, L., Zheng, S., and Zhang, W. (2020). Mechanisms of RNA N6-Methyladenosine in Hepatocellular Carcinoma: From the Perspectives of Etiology. Front. Oncol. 10, 1105. doi:10.3389/fonc.2020.01105
Ma, J. Z., Yang, F., Zhou, C. C., Liu, F., Yuan, J. H., Wang, F., et al. (2017). METTL14 Suppresses the Metastatic Potential of Hepatocellular Carcinoma by Modulating N 6 ‐Methyladenosine‐Dependent Primary MicroRNA Processing. Hepatology 65 (2), 529–543. doi:10.1002/hep.28885
Maity, A., and Das, B. (2016). N6‐Methyladenosine Modification inmRNA: Machinery, Function and Implications for Health and Diseases. Febs j 283 (9), 1607–1630. doi:10.1111/febs.13614
Mapperley, C., van de Lagemaat, L. N., Lawson, H., Tavosanis, A., Paris, J., Campos, J., et al. (2021). The mRNA m6A Reader YTHDF2 Suppresses Proinflammatory Pathways and Sustains Hematopoietic Stem Cell Function. J. Exp. Med. 218 (3), e20200829. doi:10.1084/jem.20200829
Meyer, K. D., and Jaffrey, S. R. (2014). The Dynamic Epitranscriptome: N6-Methyladenosine and Gene Expression Control. Nat. Rev. Mol. Cel Biol 15 (5), 313–326. doi:10.1038/nrm3785
Nishanth, G., Deckert, M., Wex, K., Massoumi, R., Schweitzer, K., Naumann, M., et al. (2013). CYLD Enhances Severe Listeriosis by Impairing IL-6/STAT3-Dependent Fibrin Production. Plos Pathog. 9 (6), e1003455. doi:10.1371/journal.ppat.1003455
Ondo, K., Isono, M., Nakano, M., Hashiba, S., Fukami, T., and Nakajima, M. (2021). The N6-Methyladenosine Modification Posttranscriptionally Regulates Hepatic UGT2B7 Expression. Biochem. Pharmacol. 189, 114402. doi:10.1016/j.bcp.2020.114402
Pan, X.-Y., Huang, C., and Li, J. (2021). The Emerging Roles of m6A Modification in Liver Carcinogenesis. Int. J. Biol. Sci. 17 (1), 271–284. doi:10.7150/ijbs.50003
Robinson, J. T., Thorvaldsdóttir, H., Winckler, W., Guttman, M., Lander, E. S., Getz, G., et al. (2011). Integrative Genomics Viewer. Nat. Biotechnol. 29 (1), 24–26. doi:10.1038/nbt.1754
Runyan, C. E., Schnaper, H. W., and Poncelet, A.-C. (2004). The Phosphatidylinositol 3-Kinase/Akt Pathway Enhances Smad3-Stimulated Mesangial Cell Collagen I Expression in Response to Transforming Growth Factor-β1. J. Biol. Chem. 279 (4), 2632–2639. doi:10.1074/jbc.M310412200
Shah, R., Reyes-Gordillo, K., Arellanes-Robledo, J., Lechuga, C. G., Hernández-Nazara, Z., Cotty, A., et al. (2013). TGF-β1 Up-Regulates the Expression of PDGF-β Receptor mRNA and Induces a Delayed PI3K-, AKT-, and p70S6K-Dependent Proliferative Response in Activated Hepatic Stellate Cells. Alcohol. Clin. Exp. Res. 37 (11), 1838–1848. doi:10.1111/acer.12167
Śledź, P., and Jinek, M. (2016). Structural Insights into the Molecular Mechanism of the m6A Writer Complex. Elife 5, e18434. doi:10.7554/eLife.18434
Smith-Cortinez, N., van Eunen, K., Heegsma, J., Serna-Salas, S. A., Sydor, S., Bechmann, L. P., et al. (2020). Simultaneous Induction of Glycolysis and Oxidative Phosphorylation during Activation of Hepatic Stellate Cells Reveals Novel Mitochondrial Targets to Treat Liver Fibrosis. Cells 9 (11), 2456. doi:10.3390/cells9112456
Szklarczyk, D., Morris, J. H., Cook, H., Kuhn, M., Wyder, S., Simonovic, M., et al. (2017). The STRING Database in 2017: Quality-Controlled Protein-Protein Association Networks, Made Broadly Accessible. Nucleic Acids Res. 45 (D1), D362–d368. doi:10.1093/nar/gkw937
Vig, S., Talwar, P., Kaur, K., Srivastava, R., Srivastava, A. K., and Datta, M. (2015). Transcriptome Profiling Identifies P53 as a Key Player during Calreticulin Deficiency: Implications in Lipid Accumulation. Cell Cycle 14 (14), 2274–2284. doi:10.1080/15384101.2015.1046654
Wang, X., and He, C. (2014). Dynamic RNA Modifications in Posttranscriptional Regulation. Mol. Cel 56 (1), 5–12. doi:10.1016/j.molcel.2014.09.001
Wang, X., Lu, Z., Gomez, A., Hon, G. C., Yue, Y., Han, D., et al. (2014). N6-Methyladenosine-Dependent Regulation of Messenger RNA Stability. Nature 505 (7481), 117–120. doi:10.1038/nature12730
Wang, X., Zhao, B. S., Roundtree, I. A., Lu, Z., Han, D., Ma, H., et al. (2015). N6-Methyladenosine Modulates Messenger RNA Translation Efficiency. Cell 161 (6), 1388–1399. doi:10.1016/j.cell.2015.05.014
Wang, Q., Wei, S., Zhou, H., Li, L., Zhou, S., Shi, C., et al. (2020a). MicroRNA-98 Inhibits Hepatic Stellate Cell Activation and Attenuates Liver Fibrosis by Regulating HLF Expression. Front. Cel Dev. Biol. 8, 513. doi:10.3389/fcell.2020.00513
Wang, X., Liu, X., Liu, N., and Chen, H. (2020b). Prediction of Crucial Epigenetically-Associated, Differentially Expressed Genes by Integrated Bioinformatics Analysis and the Identification of S100A9 as a Novel Biomarker in Psoriasis. Int. J. Mol. Med. 45 (1), 93–102. doi:10.3892/ijmm.2019.4392
Wang, Y., Yang, Z., Wei, Y., Li, X., and Li, S. (2021a). Apolipoprotein A4 Regulates the Immune Response in Carbon Tetrachloride-Induced Chronic Liver Injury in Mice. Int. Immunopharmacology 90, 107222. doi:10.1016/j.intimp.2020.107222
Wang, Z., Gu, J., Yan, A., and Li, K. (2021b). Downregulation of Circ-RANBP9 in Laryngeal Cancer and its Clinical Significance. Ann. Transl Med. 9 (6), 484. doi:10.21037/atm-21-567
Wu, L., Zhang, Q., Mo, W., Feng, J., Li, S., Li, J., et al. (2017). Quercetin Prevents Hepatic Fibrosis by Inhibiting Hepatic Stellate Cell Activation and Reducing Autophagy via the TGF-β1/Smads and PI3K/Akt Pathways. Sci. Rep. 7 (1), 9289. doi:10.1038/s41598-017-09673-5
Wu, X., Zhang, X., Tao, L., Dai, X., and Chen, P. (2020). Prognostic Value of an m6A RNA Methylation Regulator-Based Signature in Patients with Hepatocellular Carcinoma. Biomed. Res. Int. 2020, 2053902. doi:10.1155/2020/2053902
Wu, S., Ye, S., Lin, X., Chen, Y., Zhang, Y., Jing, Z., et al. (2021). Small Hepatitis B Virus Surface Antigen Promotes Malignant Progression of Hepatocellular Carcinoma via Endoplasmic Reticulum Stress-Induced FGF19/JAK2/STAT3 Signaling. Cancer Lett. 499, 175–187. doi:10.1016/j.canlet.2020.11.032
Yang, L., Chan, C.-C., Kwon, O.-S., Liu, S., McGhee, J., Stimpson, S. A., et al. (2006). Regulation of Peroxisome Proliferator-Activated Receptor-γ in Liver Fibrosis. Am. J. Physiology-Gastrointestinal Liver Physiol. 291 (5), G902–G911. doi:10.1152/ajpgi.00124.2006
Yang, Y., Hsu, P. J., Chen, Y.-S., and Yang, Y.-G. (2018). Dynamic Transcriptomic m6A Decoration: Writers, Erasers, Readers and Functions in RNA Metabolism. Cell Res 28 (6), 616–624. doi:10.1038/s41422-018-0040-8
Yang, W., Tao, Y., Wu, Y., Zhao, X., Ye, W., Zhao, D., et al. (2019). Neutrophils Promote the Development of Reparative Macrophages Mediated by ROS to Orchestrate Liver Repair. Nat. Commun. 10 (1), 1076. doi:10.1038/s41467-019-09046-8
Zhang, W., Sargis, R. M., Volden, P. A., Carmean, C. M., Sun, X. J., and Brady, M. J. (2012). PCB 126 and Other Dioxin-like PCBs Specifically Suppress Hepatic PEPCK Expression via the Aryl Hydrocarbon Receptor. PLoS One 7 (5), e37103. doi:10.1371/journal.pone.0037103
Zhang, K., Han, X., Zhang, Z., Zheng, L., Hu, Z., Yao, Q., et al. (2017). The Liver-Enriched Lnc-LFAR1 Promotes Liver Fibrosis by Activating TGFβ and Notch Pathways. Nat. Commun. 8 (1), 144. doi:10.1038/s41467-017-00204-4
Zhang, X., Yao, J., Shi, H., Gao, B., and Zhang, L. (2019a). LncRNA TINCR/microRNA-107/CD36 Regulates Cell Proliferation and Apoptosis in Colorectal Cancer via PPAR Signaling Pathway Based on Bioinformatics Analysis. Biol. Chem. 400 (5), 663–675. doi:10.1515/hsz-2018-0236
Zhang, Z., Chen, L.-Q., Zhao, Y.-L., Yang, C.-G., Roundtree, I. A., Zhang, Z., et al. (2019b). Single-Base Mapping of M 6 A by an Antibody-Independent Method. Sci. Adv. 5 (7), eaax0250. doi:10.1126/sciadv.aax0250
Zhang, B., Wu, Q., Li, B., Wang, D., Wang, L., and Zhou, Y. L. (2020a). m6A Regulator-Mediated Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization in Gastric Cancer. Mol. Cancer 19 (1), 53. doi:10.1186/s12943-020-01170-0
Zhang, C., Huang, S., Zhuang, H., Ruan, S., Zhou, Z., Huang, K., et al. (2020b). YTHDF2 Promotes the Liver Cancer Stem Cell Phenotype and Cancer Metastasis by Regulating OCT4 Expression via m6A RNA Methylation. Oncogene 39 (23), 4507–4518. doi:10.1038/s41388-020-1303-7
Zhang, Z., Luo, K., Zou, Z., Qiu, M., Tian, J., Sieh, L., et al. (2020c). Genetic Analyses Support the Contribution of mRNA N6-Methyladenosine (m6A) Modification to Human Disease Heritability. Nat. Genet. 52 (9), 939–949. doi:10.1038/s41588-020-0644-z
Zhong, L., Liao, D., Zhang, M., Zeng, C., Li, X., Zhang, R., et al. (2019). YTHDF2 Suppresses Cell Proliferation and Growth via Destabilizing the EGFR mRNA in Hepatocellular Carcinoma. Cancer Lett. 442, 252–261. doi:10.1016/j.canlet.2018.11.006
Zhu, B., Gong, Y., Shen, L., Li, J., Han, J., Song, B., et al. (2020a). Total Panax Notoginseng Saponin Inhibits Vascular Smooth Muscle Cell Proliferation and Migration and Intimal Hyperplasia by Regulating WTAP/p16 Signals via m6A Modulation. Biomed. Pharmacother. 124, 109935. doi:10.1016/j.biopha.2020.109935
Keywords: m6A methylation, m6A-seq, liver fibrosis, HSCs, WTAP
Citation: Fan C, Ma Y, Chen S, Zhou Q, Jiang H, Zhang J and Wu F (2021) Comprehensive Analysis of the Transcriptome-Wide m6A Methylation Modification Difference in Liver Fibrosis Mice by High-Throughput m6A Sequencing. Front. Cell Dev. Biol. 9:767051. doi: 10.3389/fcell.2021.767051
Received: 30 August 2021; Accepted: 01 November 2021;
Published: 16 November 2021.
Edited by:
Chengqi Yi, Peking University, ChinaReviewed by:
Matthew Wong, Children’s Cancer Institute Australia, AustraliaZhiwen Fan, Nanjing Drum Tower Hospital, China
Copyright © 2021 Fan, Ma, Chen, Zhou, Jiang, Zhang and Wu. 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: Hui Jiang, jianghui@ahtcm.edu.cn; Jiafu Zhang, zyfyzjf@163.com; Furong Wu, Fr13866767052@163.com