N⁶‐methyladenine profiling of low‐input multiplex clinical samples on transcriptome reveals RNA modifications implicated in type 2 diabetes and acute myocardial infarction

RNA N6-methyladenines (m6A) not only play an essential role in normal biological processes but also participate in the pathogenesis of diseases including type 2 diabetes (T2D) and comorbidities.1, 2 However, conventional technologies (e.g. m6A-seq)3 require large amounts of input materials often prohibitive for clinically-feasible biospecimens (e.g. a few millilitres of blood or tissue from fine needle aspiration), we, therefore, developed a Strategy of Low-Input Multi-barcode m6A-seq (SLIM-m6A-seq) to facilitate transcriptome-wide m6A profiling with limited RNA. Using this new technology, we implicated novel m6A modifications and pathways in the development of acute myocardial infarction (AMI) in patients with T2D. Technically, SLIM-m6A-seq exploited the barcoding strategy to encode multiple samples with unique barcodes before pooling, library construction, and sequencing (Figure 1A and Table S1).4, 5 The resulting profiles can be split by barcodes (Table S2) and traced back to the original samples. We first tested the feasibility of SLIM-m6A-seq in HeLa cells and calibrated them with spike-in oligonucleotides (Note S1 and Figure S1). Then we compared SLIM-m6A-seq of a pooled sample containing barcoded mRNA varying from 200 to 10 ng with m6A-seq of 800 ng mRNA (Note S1). The SLIM-m6A-seq results from input levels of 200 and 100 ng showed similar distribution (e.g. chromosome 1 in Figure S2A) and complexity with 2.56% relative standard deviation (RSD), compared to the data from m6A-seq in 800 ng input, albeit substantially different starting amounts (Figure S3A). We further validated the capability of SLIM-m6A-seq for low-input samples in two 25-sample sets containing 10 ng (Init10) or 5 ng (Init5) fragmented mRNA from HEK293T cells. Overall, greater than 90% mapping rates on transcriptome were achieved for these samples (Table S3). A similar number of m6A peaks was detected, with an RSD of 34.26% in Init10 and an RSD of 37.43% in Init5 (Figure S4A), as well as a consistent number of m6A-containing genes across different barcodes (Figure S4B). The m6A profiles across the two sets were highly correlated (Figure 1B and Figure S4C,D), thus confirming technical reproducibility in low-input samples. The barcoded samples in Init5 or Init10 sets exhibited high reproducibility of complexity with RSD ≦ 20% in both sets with or without m6A immunoprecipitation (IP), although the complexity would be reduced with smaller input as expected (Figure S4E). The enrichment folds of m6A in individually barcoded samples were comparable between both sets – 6.22-folds in Init5 and 7.51-folds in Init10 (Figure 1C and Figure S5A). The RRACH motifs (Figure 1D and Figure S5B) and consistent m6A distributions (Figure 1E) were discovered across all samples, as shown in an example of chromosome 1 (Figure 1F). We then profiled a cohort of 20 blood samples (Figure S6 and Table S4) from patients with T2D (n = 10) and T2D with AMI (T2D+AMI) (n = 10) matched for age, sex and blood pressure to demonstrate the clinical feasibility of SLIM-m6A-seq. A similar complexity was observed in all samples with the ratio of distinct/total reads ranging from 0.84 to 0.90 and an RSD of 1.74% (Figure 2A). The principal component analysis suggested transcriptome-wide differences between these two groups (Figure 2B). In total, 3223 dysregulated genes (|log2 (fold change)| ≥ 1, adjusted p < 0.05) were detected between these two groups (Figure 2C), with 1845 up-regulated and 1378 down-regulated in T2D+AMI (Figure 2D). A panel of 14 differential genes was further established through the leave-one-out cross-validation and the elastic net regularization (Figure 2E and Table S5). Note that, 7234 m6A peaks were found with typical distributions (Figure 2F) and RRACH motifs (Figure 2G). SLIM-m6A-seq revealed an elevated overall m6A modification level in T2D+AMI patients (Figure 2H), which was confirmed by high-performance liquid chromatography-tandem mass spectrometry (Figure 2I and Table S6) and also reflected by the fact that 806 genes with increased and 260 genes with reduced m6A modification levels in T2D+AMI patients (Figure 2J). Notably, the expression level of the m6A demethylase FTO decreased evidently in T2D+AMI (Figure 3A), supported by real-time quantitative polymerase chain reaction (RT-qPCR) (Figure 3B and Tables S6 and S7). The conjoint analysis uncovered 392 genes exhibiting different combinations of expression and m6A modification patterns (|log2 (fold change)| ≥ 1, adjusted p < 0.05 for expression change, and ≥ 1.68-fold enrichment for m6A). Of these intersecting genes, 311 genes exhibited concordant direction of changes, and 81 showed contrary alterations (Figure 3C). Gene ontology (GO) analysis of the 392 overlapping genes revealed that the biological processes related to T2D+AMI were involved in protein kinase activities, inflammatory response, apoptosis, and angiogenesis (Figure 3D),6 which have been implicated in vascular biology, inflammation, and myocardial injury (Table S5).7 Of the 392 intersecting differential genes, several inflammation- or atherosclerosis-related genes and angiogenesis activators showed correlation (Pearson |R| ≥ 0.5) with one or more clinical cardiac and inflammatory biomarkers (Figure 3E). Of note, TNFSF14, a credible upregulated gene explored in our differential gene analysis and relevant to T2D and angiocardiopathy,8, 9 showed reduced m6A in T2D+AMI. This observation indicated that the effect of m6A in the progression of AMI during T2D could be through regulating the expression of disease-related genes. The GO analysis was further supported by the observed changes in clinical indexes relevant to myocardial injury and inflammation between T2D+AMI and T2D (Figure 3F). Together with mediators in apoptosis and angiogenesis, these differential m6A-containing genes were relevant to inflammation and cardiomyocyte death,10 suggesting the significant role of m6A in the occurrence of T2D+AMI through these mechanisms. In summary, we developed SLIM-m6A-seq to implement transcriptome-wide m6A profiling with limited RNA. Comprehensive analysis of SLIM-m6A-seq data in clinical samples from patients with T2D and T2D+AMI provided novel insights into the epitranscriptome of T2D+AMI regarding the contribution of m6A. Although with certain limitations (e.g. potential confounding due to immune cell composition), the overall m6A levels in T2D+AMI appeared to be elevated compared to T2D, accompanied by down-regulation of FTO. Functional annotation analysis revealed relevant biological processes involved in the development of T2D+AMI. The current study demonstrated the technical robustness of SLIM-m6A-seq and importantly the feasibility of this new method in obtaining transcriptome-wide m6A from clinically convenient samples, thus opening up opportunities for exploiting m6A as an effective clinical approach for disease diagnosis and prognosis. We would like to thank experts in the department of cardiology (Professor Yanggan Wang, Hairong Wang and Jianlei Cao) and the department of endocrinology (Professor Haohua Deng) of Zhongnan Hospital of Wuhan University for disease diagnosis and critical discussion. The numerical calculations in this article were performed using the supercomputing system in the Supercomputing Center of Wuhan University. This work was supported by the National Natural Science Foundation of China (92253202, 91940304 and 22177087 to Xiaocheng Weng; 92153303, 22037004 and 91940000 to Xiang Zhou; 82172359, 81972009 and 81772276 to Song-Mei Liu; and 22277094 to Fang Wang), Hubei Provincial Natural Science Fund for Creative Research Groups (2019CFA018), and Translational Medicine and Interdisciplinary Research Joint Fund of Zhongnan Hospital of Wuhan University (ZNJC201914). The authors declare that they have no conflict of interest. Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.