Human m6A-mRNA and lncRNA epitranscriptomic microarray reveal function of RNA methylation in hemoglobin H-constant spring disease

The thalassemia of Hemoglobin H-Constant Spring disease (HbH-CS) is the most common type of Thalassemia in non-transfusion thalassemia. Interestingly, the clinical manifestations of the same genotype of thalassemia can be vastly different, likely due to epigenetic regulation. Here, we used microarray technology to reveal the epigenetic regulation of m6A in modifiable diseases and demonstrated a role of BCL2A1 in disease regulation. In this study, we revealed that methylating enzyme writers including METTL16, WTAP, CBLL1, RBM15B, and ZC3H13 displayed low expression and the demethylating enzyme ALKBH5, along with reader proteins including IGF2BP2 and YTHDF3 exhibited high expression. In addition, BCL2A1 was hypo-methylated and showed low expression. We also revealed that the BCL2A1 methylation level and IGF2BP2 expression were negatively correlated. Additionally, the mRNAs expression between ALKBH5 and IGF2BP2 were positively correlated. In HbH-CS, most genes were hypo-methylated. This included BCL2A1, which may play an important role in the process of red blood cell differentiation and development of HbH-CS. Moreover, the mRNA-M6A methylation status may be regulated by the demethylating enzyme ALKBH5 via IGF2BP2.


Differential expression of genes in T versus N.
Based on our analysis, we selected 8 differentially expressed mRNAs, based on their m 6 A level, and confirmed their expression patterns in T (n = 15) and N (n = 16) samples, using RT-qPCR. The selected mRNAs were Mettl16, WTAP, CBLL1, RBM15B, ZCH3H13, IGF2BP2, YTHDF3, and ALKBH5. As depicted in Fig. 3A, these genes expressed differently in T versus N. However, their expression patterns were consistent with the epitranscriptomic microarray sequences (Fig. 3B), thereby confirming the reliability of the microarray technique. The primers for our RT-qPCR work are presented in Table 2.
Gene ontology enrichment and pathway analysis. To explore the functional correlations between the differentially m 6 A-methylated and differentially expressed mRNA, we enriched select gene ontological functions and GO terms (http:// www. geneo ntolo gy. org) and used the bioconductor top GO package for analysis. Among the top 10 GO were biological process (BP), cellular component (CC), and molecular function (MF). Moreover, we determined significance using Fisher Exact test p-value and established the enrichment score via -log10 (p) formula. Figures 4A illustrates the differentially m 6 A hypo-methylated mRNA. In addition, we conducted metabolic pathway analysis, using KEGG pathways, on the differentially m 6 A hypo-methylated mRNAs, with statistical analysis as described before. As depicted in Fig. 4B, the pathway enrichment analysis showed involvement of essential pathways like the Herpes simplex virus 1 infection, sphingolipid signaling pathway, NF-kappa B axis, Th17 cell differentiation, B cell receptor axis, Viral myocarditis, Yersinia infection, osteoclast differentiation, phospholipase D axis, and the AGE-RAGE axis. Table 1. Clinical information of all participants. N represents healthy volunteers (n = 16) and T refers to the HbH CS thalassemia patients (n = 15). T-test was used for analysis. No discernible difference was observed in age and gender within the N and T cohorts.  Table 3 and 20 differentially hypo-methylated m 6 A non-mRNAS including lncRNA and other small RNAs are listed in Table 4. Here, we evaluated m 6 A levels of BCL2A between T (n = 10) and N (n = 10). The primary methylated sites were in the CDS and 5'UTR regions 21 . Figure 5A illustrates the methylated mRNA and positions of methylation. Based on our data, we demonstrated that the mRNA expression and m 6 A levels of BCL2A1 were markedly down-regulated in T versus N (Fig. 5B,C). The primers used for BCL2A1 identification were 5′AGA ATC TGA AGT CAT GCT TGGA3′and 5′CTC CTT TTC CAT CAC TTG GTTG3. In addition, using KEGG analysis, we showed that the methylation of BCL2A1 was linked to IGF2BP2 mRNA expression (Fig. 5D). Figure 5E showed that ALKBH5 is   www.nature.com/scientificreports/ positively correlated with IGF2BP2. To investigate the the potential roles of IGF2BP2, we examined the effects of IGF2BP2 knockdown in K562 cells and confirmed that IGF2BP2 was knocked down by using siRNA sequences (Fig. 5F,I). We found that knockdown of IGF2BP2 significantly inhibited the ALKBH5 and BCL2A1 expression of k562 cells ( Fig. 5G-I).

Discussion
In recent years, it became evident that people carrying the same Globin genotype can have vastly different phenotypes. Moreover, thalassemia can aggravate in presence of stressors like pregnancy, infection, surgery and so on. Based on these facts, it is possible that thalassemia may be affected by other, non-linked genes 19 . Multiple reports have confirmed that cis-regulatory elements like DNA methylation, histone acetylation, and micro RNAs can regulate pathogenesis and clinical heterogeneity of thalassemia [19][20][21] . In addition, other modifiers were also shown to affect thalassemia phenotype, such as mutations in the molecular chaperone haemoglobin stabilizing protein (AHSP) and transcription regulator GATA1. These emerging reports offer novel view into the therapeutics and prevention of thalassemia, including prevention of the disease in thalassemia heterozygote carriers 22 . Recent studies on m 6 A-mediated RNA modification demonstrated the dynamic reversibility of RNA modification in regulating RNA metabolism, processing, and directional differentiation of stem cells 23 . Moreover, the m 6 A methylation was shown to modulate the coordinated efforts of transcription and post-transcriptional expression 24 . This includes events like mRNA splicing, export, localization, translation, and stability. M 6 A-mediated methylation of mRNA can potentially enable interaction with modulatory proteins which can regulate downstream gene expression. Emerging evidences suggest a crucial role of m 6 A-mediated mRNA modification in driving mRNA translation 25 .  www.nature.com/scientificreports/ To determine whether m 6 A modification is involved in producing the thalassemia phenotype, we assessed immature erythrocytes from the peripheral blood of 15 T and 16 N. We identified that methylating enzymes like METTL16, WTAP, CBLL1, RBM15B, and ZC3H13 were scarcely expressed and the demethylating enzyme ALKBH5, along with reading protein-related genes, such as, IGF2BP2 and YTHDF3, were highly expressed in T versus N. In addition, we detected differential regulation of other m 6 A-methylated RNAs like RBM15, YTHDF3, IGF2BP1, IGF2BP3, and multiple mRNAS, micRNA, and LNC RNA, which were hypo-methylated in T versus N. Given these evidences, we believe that the m 6 A-mediated methylation is crucial for the pathogenesis of HbH-CS. BCL2A1 (B cell lymphoma 2 related A1), from the BCL2 (B cell lymphoma 2) protein family 26 , counters cell apoptosis via prevention of cytoplasmic accumulation of cytochrome and prevents the subsequent stimulation of the internal apoptotic axis. This protein is often highly expressed in advanced cancers and denotes poor prognosis 27 . In this study, we verified that the BCL2A1 mRNA was hypo-m 6 A-methylated and had low expression  www.nature.com/scientificreports/ in T versus N. However, the underlying mechanism is yet to be determined. M 6 A methylation can stabilize target mRNA, thereby affecting downstream gene expression 28,29 . M 6 A methylation is generally common in proteincoding transcripts and are localized in the 3′UTRs 29 and 5′UTR 30 regions. Interestingly, the methylation sites of BCL2A1 are located within the CDS 21 and 5'UTR regions 18 . Based on our data, we propose that alterations in the mRNA methylation status contributes to thalassemia. Due to the severe down-regulation of BCL2A1 in HbH-CS patients, apoptosis remains uninhibited, leading to hemolytic anemia. In thalassemia, erythropoiesis is finely regulated by a complex network of transcription factors, including those involved in erythropoiesis like the erythropoietin receptor EPOR, glycophorin, and Globin peptide chains. In addition, the anti-apoptotic protein BCL-x, induced by the transcription factors STAT5 and GATA1, are activated by erythropoietin EPO 31 . When BCL2A1 levels diminish during erythroid development, further erythrocytes production ceases, thereby aggravating thalassemia. Based on our analysis, we detected BCL2A1 in 3 separate signaling pathways, namely, NF-KappaB, Acute Myeloid Leukemia and transcriptional misregulation in cancer (Supplemental Information). Interestingly, we also discovered that the m 6 A methylation status of BCL2A1 was negatively correlated with IGF2BP2 gene expression. However, m 6 A methylation was not directly affected by the M 6 A-associated enzymes. Alternately, ALKBH5 is positively correlated with IGF2BP2. Collectively, based on our analysis, we propose that ALKBH5 regulates RBC differentiation and development by altering the methylation status of BCL2A1 via IGF2BP2. However, this requires further study and confirmation.

Conclusion
M 6 A plays a significant role in thalassemia. In this study, we have, for the first time, explored the m 6 A RNA methylation status and its regulation of RNA stability in HbH-CS patients. Moreover, using m 6 A-RIP-seq and RNA-seq data, we established a profile of differentially methylated mRNA in T versus N samples. Lastly, we proposed involvement of m 6 A-regulated BCL2A1 in the pathogenesis of thalassemia. To confirm our preliminary data, more investigation is necessary to explore the methylation status of relevant RNAs and modulation of target downstream genes during erythroid differentiation.
Our work had certain deficiencies. Firstly, we discovered that there were two m 6 A methylation sites in BCL2A1. However, we did not elucidate which site holds more significance to thalassemia. Secondly, the underlying mechanism behind RNA methylation of BCL2A1 was not examined, and needs further exploration.

Materials and methods
Participants and samples. We collected samples from 16 HbH-CS thalassemia patients and 15 healthy volunteers, between March and July 2020. The HbH-CS patients had differing degrees of moderate anemia and hepatosplenomegaly and did not receive any blood transfusion in the past 3 months. The HGB range, among the HbH-CS patients, were between 56 and 103 g/l. In addition, Doppler's ultrasound of HbH-CS patients revealed splenomegaly ranging from Grade 1 to Grade 3. HbH-CS thalassemia and healthy volunteers will hereby be referred to as group T and group N, respectively. We conducted the Arraystar Human m 6 A-mRNA and lncRNA epitranscriptomic microarray analysis of 5 pairs of immature erythrocytes, particularly, 5 from HbH-CS thalassemia (T) and 5 from healthy volunteers (N). Our works received ethical approval from the First Affiliated www.nature.com/scientificreports/ www.nature.com/scientificreports/ Hospital of Guangxi Medical University. All participants were recruited from the same university and agreed to sign informed consent forms. If subjects were under 18, the informed consents were signed by their parent and/or legal guardian. All methods were performed in accordance with the Declaration of Helsinki. All HbH-CS patients were confirmed of their diagnosis with blood routine, haemoglobin electrophoresis, and DNA analysis 32,33 . All N groups were without any blood-related diseases. All participants were between the ages of 2 and 50. Patient demographics are provided in Table 1.

Sorting of immature red blood cells (RBCs)
. Immature RBCs were collected by sorting peripheral blood from the T and N cohorts, using a positive CD71 (CD71 Microbeads human, Miltenyi Biotec GmbH, Germany) selection method with a magnetic shelf (MiniMACS Starting Kit, Miltenyi Biotec GmbH, Germany), followed by confirmation with flow cytometry 34,35 . The flow cytometry results are depicted in Fig. 1. The sorted RBC samples were subsequently maintained in TRIZOL at − 80 °C for further analysis.
Total RNA isolation and RT-qPCR. Total  Two-color RNA labeling and array hybridization. The modified RNAs were eluted from the immunoprecipitated magnetic beads as "IP" (immunoprecipitated, Cy5-labelled). The unmodified RNAs were recovered from the supernatant as "Sup" (supernatant, Cy3-labelled). The "IP" and "Sup" RNAs were then labeled with Cy5 and Cy3 respectively, in separate reactions, using Arraystar Super RNA Labeling Kit. The RNAs were combined together and hybridized onto Arraystar Human mRNA & lncRNA Epitranscriptomic Microarray (8 × 60 K, Arraystar). After washing the slides, the arrays were scanned in the two-color channels by an Agilent Scanner G2505C.
Cell culture and transfection. K562  Data analysis. The Agilent Feature Extraction software (version 11.0.1.1) was employed for the analysis of acquired array images. Raw intensities of Modified RNA (Cy5-labelled) and Unmodified RNA (Cy3-labelled) were normalized with the mean of log2-scaled Spike-in RNA intensities. Next, signals with Present (P) or Marginal (M) QC flags in a minimum of 1 in 10 samples were marked as "All Targets Value" in Excel for further "m 6 A methylation level", "m 6 A quantity" and "expression level" calculations. The "m 6 A methylation level" www.nature.com/scientificreports/ was measured via the percentage of modification based on sample normalized intensities. The "m 6 A quantity" was assessed from the Modified RNA (Cy5-labelled) normalized intensities. The "RNA expression level" was analyzed according to the normalized intensities of all RNA. Moreover, we performed additional quartile normalization using the limma package to ascertain array expression, before probe flag screening. Differentially m 6 A-methylated RNAs or differentially expressed RNAs between the HbH-CS and healthy volunteers were recognized by fold change and statistical significance (p-value) values. Lastly, hierarchical clustering was done to show degrees of m 6 A-methylation within samples. P value < 0.05 denotes significance.
Ethics approval and consent to participate. Our works received ethical approval from the First Affiliated Hospital of Guangxi Medical University. All participants agreed to sign informed consent forms.

Consent for publication.
All authors agree to the publication of this article.

Data availability
The following information was supplied regarding data availability. The sequencing data are available in the Supplemental Files. www.nature.com/scientificreports/