Transcriptome Sequencing Analysis of Peripheral Blood of Type 2 Diabetes Mellitus Patients With Thirst and Fatigue

Purpose: The purpose of this study is to explore the differences in transcriptome expression profiles between healthy subjects and type 2 diabetes mellitus patients with thirst and fatigue (D-T2DM) and, in addition, to investigate the possible role of noncoding ribonucleic acids (RNAs) in the pathogenesis of D-T2DM. Methods: We constructed the expression profiles of RNAs by RNA sequencing in the peripheral blood of D-T2DM patients and healthy subjects and analyzed differentially expressed RNAs. Results: Compared with healthy subjects, a total of 469 mRNAs, 776 long non-coding RNAs (lncRNAs), and 21 circular RNAs (circRNAs) were differentially expressed in D-T2DM patients. Furthermore, several genes associated with insulin resistance, inflammation, and mitochondrial dysfunction were identified within the differentially expressed mRNAs. Differentially expressed lncRNAs were primarily involved in biological processes associated with immune responses. In addition, differentially expressed circRNAs may target miRNAs associated with glucose metabolism and mitochondrial function. Conclusions: Our results may bring a new perspective on differential RNA expression involved in the pathogenesis of D-T2DM and promote the development of novel treatments for this disease.


INTRODUCTION
In today's world, with the improvement of living conditions and the change of lifestyle, more and more diabetic patients are diagnosed. The number of people with diabetes is expected to reach 700 million by 2045 (1). In traditional Chinese medicine, diabetes is divided into different types according to the clinical symptoms such as fever, fatigue, sputum, and urinary frequency (2). The "dual deficiency of Qi and Yin syndrome" is a common type of diabetes, and the main symptoms of this type are dry mouth and throat and weakness with general fatigue. Similarly, a proportion of patients with T2DM are accompanied with thirst and fatigue in modern medicine (3)(4)(5)(6)(7). Therefore, the diagnosis and treatment of this type of T2DM (D-T2DM) have become an important issue in clinical work.
Non-coding ribonucleic acids (ncRNAs) are transcripts that cannot be translated into proteins, but are generally considered to have a role in regulating protein expression (8). ncRNAs are divided into long non-coding RNAs (lncRNAs) and circular loops RNAs (circRNAs) according to their morphology (9). lncRNAs are transcripts longer than 200 nucleotides that have been implicated in diverse biological functions, such as transcriptional and posttranscriptional regulation and chromatin modification (10)(11)(12). CircRNAs are a type of circular ncRNA discovered recently and attracted wide attention (13). In our previous research, ncRNAs have been shown to play a key role in the pathogenesis of T2DM and other endocrine diseases (14,15). However, the role of ncRNAs in D-T2DM development remains unclear. Therefore, comprehensive detection and analysis of ncRNAs in the development of D-T2DM are critical for the prevention and treatment of D-T2DM.
In this study, we used RNA sequencing technology to construct RNA expression profiles in peripheral blood of patients with D-T2DM. Differentially expressed mRNAs, lncRNAs, and circRNAs were detected, and their functions were predicted by bioinformatics analysis, in order to discover new targets and provide assistance for the treatment of D-T2DM.

Ethics Statement and Information of Subjects
This research was approved by the Ethics Committee of Beijing University of Chinese Medicine (BUCM) (2017BZHYLL0105). All subjects, enrolled from the Affiliated Hospital of BUCM, agreed to participate in the study after fully understanding the purpose and procedure of the experiment. Full inclusion and exclusion criteria are listed in Table 1. The subjects were divided  into two groups: D-T2DM group (ID: QYD1, QYD2, QYD3,  QYD4, QYD5, and QYD6) and the healthy subjects as the control  group (ID: LZC001, LZC002, LZC003, LZC004, LZC005, and  LZC006). Subsequently, the fasting venous blood of all subjects were obtained and stored at −80 • C until analysis.

Total RNA Extraction, Library Construction, and Illumina Sequencing
According to the previous experimental methods, we extracted the total RNA and constructed the RNA library (16). Briefly, the total RNA was extracted by TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and tested for purity using the NanoPhotometer spectrophotometer (IMPLEN, CA, USA). In addition, the concentration and integrity of RNA were assessed using the Qubit 2.0 fluorometer (Life Technologies, CA, USA) and the Agilent 2100 bioanalyzer (Agilent Technologies, CA, USA), respectively. Next, Ribo-Zero rRNA Removal Kit (Epicentre, USA) was used to remove ribosomal RNA (rRNA). RNA-sequencing libraries were prepared using the rRNA-depleted RNA by NEBnext ultra RNA library prep kit following the manufacturer's instruction. After that, the quality of libraries was determined using the Agilent

Quality Control and Mapping
Raw data were processed through in-house Perl scripts. In this step, the raw data were cleaned by removing reads containing adaptors, contaminants, and low-quality reads. All clean data have been submitted to Sequence Read Archive with accession number SRP274496. Additionally, the Q20, Q30, and GC content were calculated to estimate the quality of clean reads. Next, the clean reads were aligned to the reference genome (GRCh37/hg19) using STAR (v2.5.1b).

Quantitative Analysis of Genes and ncRNAs
The FPKM of transcripts and ncRNAs were calculated by Cuffdiff (v2.1.1) (17). FPKM refers to the number of fragments per kilobase length from a gene per million fragments mapped. It considers the effect of both sequencing depth and gene length on fragments count.

GO and KEGG Pathway Annotation
In order to understand the biological function of differentially expressed genes, we performed enrichment analysis by topGO software. Based on the newest KEGG database, the pathway analysis was performed to determine the significant pathway of the differential genes. Fisher test was used for enrichment analysis, and those with P < 0.05 were significantly enriched.

Interaction Network Analysis
Cytoscape v2.8.2 software (http://www.cytoscape.org/) was used to construct the lncRNA-mRNA and circRNA-miRNA regulatory network based on the differentially expressed gene data in the blood between D-T2DM and healthy subjects.

Statistical Analysis
This study used GraphPad Prism 7 (GraphPad Software, CA) and SPSS software (version 20.0) for the statistical evaluations. The results are expressed as mean ± SEM. Statistical differences were determined by Student independent t-test, and the significance was accepted at P < 0.05. Each experiment was repeated for three technical replicates.

Clinical Characteristics of the Participants
In this study, six D-T2DM patients and six healthy subjects were enrolled. All patients fulfilled the diagnostic criteria for D-T2DM. There was no significant difference in age between the D-T2DM and control groups. The characteristics of all subjects are shown in Table 2. Compared with the control group, the D-T2DM group exhibited significant increases in FPG and TG levels.

Quality Assessments and Mapping Results
To construct expression profiles of mRNAs and ncRNAs of D-T2DM and healthy control, transcriptome data sets were generated by RNA-seq. Subsequently, quality control and mapping analysis were performed for the sequencing output (Supplementary Tables 1, 2). The base percentage of Q20 was >95.29%, Q30 was >88.92%, and the average mapping rate of the 12 samples was 95.21%. These results could indicate the high quality of transcriptome sequencing data with suitable mapping.

Differentially Expressed mRNAs, lncRNAs, and circRNAs
Sequencing technique was used to detect differentially expressed ncRNAs in the peripheral blood of D-T2DM patients. In Figure 1, the results showed that a total of 469 differentially expressed mRNAs (341 up-and 128 down-regulated), 776 differentially expressed lncRNAs (688 up-and 88 downregulated), and 21 differentially expressed circRNAs (5 up-and 16 down-regulated) were detected in D-T2DM patients compared with healthy control (fold change >1.5, P < 0.05). Tables 3-5 list the top 10 up-and down-regulated mRNAs, lncRNAs, and circRNAs, respectively. Furthermore, Figures 1A-F shows the volcano plot, cluster map of differentially expressed mRNAs, lncRNAs, and circRNAs, respectively. Among these genes, LRRC19, GCNT3, and CKMT2 were associated with the pathogenesis of T2DM. Furthermore, mRNAs involved in regulation of mitochondrial function and lipid metabolism were significantly expressed, such as MT-ND1, MT-ND2, and OSBL6.
lncRNA Target Gene Prediction lncRNA may perform its function by regulating genes. Therefore, we predicted the biological function of lncRNA by its colocated and coexpressed genes. We set the threshold of the colocated genes to 100 kb upstream and downstream of lncRNA; mRNA gene with an absolute value of Pearson correlation coefficient >0.95 was defined as lncRNA coexpressed mRNAs (Supplementary Tables 3, 4).

GO and KEGG Enrichment Analysis of Differentially Expressed RNAs in D-T2DM
We performed functional enrichment analysis of differentially expressed RNAs in bioinformatics databases, including GO and KEGG analysis. As shown in Figure 2A, GO analysis of differentially expressed mRNAs revealed that the most significantly enriched biological process (BP) were disruption of cells of other organism and killing of cells of other organism, and the most significantly enriched cellular component (CC) were actin cytoskeleton, cytosolic ribosome, and secretory granule. The GO analysis of differentially expressed lncRNAs showed that the most significantly enriched BPs were adaptive immune response, positive regulation of lymphocyte activation, and response to virus, and the most noteworthy enriched CCs were MHC class II protein complex, cytosolic part, and MHC protein complex. The most significantly enriched molecular functions (MFs) were MHC protein binding, tumor necrosis factor receptor binding, and antigen binding ( Figure 2C). GO analysis of differentially expressed circRNAs revealed that the most significantly enriched BPs were biological regulation, CC organization or biogenesis, and cellular process, and the most notable enrichment of CCs were cell, cell junction, and cell part. The most significantly enriched MFs were binding, catalytic activity, and channel regulator activity ( Figure 3A). KEGG analysis of differentially expressed mRNAs revealed that the most significantly enriched pathway in D-T2DM pathogenesis was ribosome ( Figure 2B). The colocated mRNAs of differentially expressed lncRNAs were significantly enriched in antigen processing and presentation, osteoclast differentiation, and leishmaniasis ( Figure 2D). Furthermore, differentially expressed circRNAs-derived genes were mostly involved in the ubiquitin-mediated proteolysis ( Figure 3B).

Protein-Protein Interaction Network of lncRNA Colocated mRNA Corresponding Genes
We extracted the interaction relationship of differential gene sets to build a network forms the STRING Protein Interaction Database (http://string-db.org/) and then imported them into   Cytoscape software for visual editing (Figure 4). As shown in Table 6, we listed the top 10 BP terms enriched for the genes involved in the protein-protein interaction network. UBA52 had the highest degree of network connectivity and enriched in the BP terms such as immune system process, regulation of immune response, and regulation of immune system process.

Regulatory Network of circRNA and miRNA
Using differentially expressed circRNA as the center and miRNA as the target, we constructed the circRNA-miRNA regulatory network of up-regulated and down-regulated circRNAs, respectively ( Figure 5)

DISCUSSION
In this study, we used sequencing technology to determine the significantly differentially expressed mRNAs, lncRNAs, and circRNAs in the peripheral blood of patients with D-T2DM. Afterward, we performed GO and KEGG pathway analysis on these differentially expressed RNAs to predict their potential biological functions. In addition, we also constructed proteinprotein interaction network and circRNA-miRNA regulatory network. Our results indicated that ncRNAs may play a role in the pathogenesis of D-T2DM and provide some potential targets for the treatment of D-T2DM.
The sequencing results showed that there were 469 mRNAs, 776 lncRNAs, and 21 circRNAs significantly dysregulated in D-T2DM patients compared with healthy subjects. In upregulated mRNAs, we have found some targets related to the pathogenesis of T2DM, such as LRRC19, GCNT3, and CKMT2. Among them, LRRC19 could activate nuclear factor κB (NF-κB) and induce expression of pro-inflammatory cytokines (18). GCNT3 was also related to inflammation (19)(20)(21). Inflammation can lead to a cluster of chronic metabolic disorders such as insulin resistance, obesity, type 2 diabetes, and cardiovascular disease (22-24). In addition, previous studies have revealed that NF-κB and its target inflammatory factor genes such as IL-1 and IL-6 played a key role in the development of insulin resistance and T2DM (25)(26)(27). CKMT2, creatine kinase, mitochondrial 2, was an effective modulator of ATP synthasecoupled respiration (28). Among down-regulated mRNAs, there were some mitochondrial genes, for instance, MT-ND1, MT-ND2, MT-ND5, and MT-ND4L. These mitochondrial genes were associated with fatty acid metabolism, mitochondrial oxidative phosphorylation, mitochondrial energy transduction, and the diabetes mellitus pathogenesis (29)(30)(31). Mitochondria are the most essential energy production organelles, supplying energy for cell metabolism in the form of ATP (32). In traditional Chinese medicine, fatigue is a major symptom of D-T2DM patients. Meanwhile, the hallmark symptom of mitochondrial dysfunction is fatigue (33). Therefore, mitochondrial dysfunction may play a role in the pathogenesis of D-T2DM. Furthermore, OSBPL6 was positively correlated with plasma levels of high-density FIGURE 2 | GO and KEGG analysis of mRNAs and lncRNAs in peripheral blood sample of D-T2DM patients. The scatterplots of GO enrichment results of differentially expressed mRNAs (A) and lncRNAs (C). The scatterplots of KEGG pathway analysis of differentially expressed mRNAs (B) and lncRNAs (D). The degree of GO and KEGG enrichment is assessed by enrichment of factors, Q-values, and number of genes. When the rich factor is larger, the Q-value is closer to zero, and the more the number of genes, the more significant the enrichment.
lipoprotein cholesterol (34). SLC12A1 may play a role in glucoseinduced insulin secretion (35). In general, the results indicated that the pathogenesis of D-T2DM may be related to the glucose and lipid metabolism disorders, occurrence of inflammation, and mitochondrial dysfunction.
At present, the mechanism of interaction between lncRNA and protein-coding genes is not clear. We predicted the biological function of lncRNA through its colocated and coexpressed protein-coding genes. Among the differentially expressed lncRNAs, lncRNA-SLC12A1, lncRNA-OSBPL6, and lncRNA-GCNT3 were colocated with SLC12A1, OSBPL6, and GCNT3, respectively. Therefore, these lncRNAs may play a role by regulating the genes related to glucose metabolism, lipid metabolism, and inflammation. Afterward, this study used GO and KEGG pathway analyses to analyze the biological function and pathways of lncRNAs-related genes in the peripheral blood of patients with D-T2DM.
The results showed that the most enriched GO term was adaptive immune response. Previous studies have demonstrated that adaptive immune factors was recognized as important etiological components in the development of insulin resistance (36). In our future work, the specific role of lncRNAs predicted by GO and KEGG analysis in the pathogenesis D-T2DM needs to be further studied. Protein-protein interaction network represents the interaction of the protein products of the lncRNA colocated genes and was used to predict the biological function of differentially expressed lncRNAs. The results showed that the proteins were mainly enriched in BPs related to immune response. In addition, UBA52 had the highest degree of connectivity, and previous study suggested that UBA52 may be implicated in the diabetic nephropathy (37). Therefore, differentially expressed lncRNAs appear to function in the pathogenesis of D-T2DM by regulating immune response.  In order to discover the molecular mechanism of circRNAs in D-T2DM, we constructed the circRNA-miRNA regulatory network based on the sequencing data. In our results, MiR-877-3p was combined with down-regulated circRNA novel_circ_0016196, novel_circ_0016198, and novel_circ_0005686. Xie et al. (38) found that MiR-877-3p was deregulated in type 2 diabetic kidney disease. Another study revealed that MiR-877-3p exerts its effects via the Blc-2-mediated mitochondrial apoptotic pathway (39). CircRNAs generally act as an miRNA sponge to regulate gene expression. Therefore, the down-regulation of these circRNAs may improve the function of MiR-877-3p and thus mediate mitochondrial apoptosis. On the other hand, up-regulated circRNA has_circ_0002590 and novel_circ_000372 were associated with MiR-149-5p. The overexpression of MiR-149-5p could ameliorate the high glucose-induced injury in human umbilical vein endothelial cells, whereas the inhibition of MiR-149-5p could suppress cell viability, induce cell apoptosis, and inhibit insulin secretion (40,41). In this research, the up-regulated  circRNAs competitively bound with MiR-149-5p to inhibit its function. In general, circRNAs may play a role in the pathogenesis of D-T2DM by regulating the function of MiR-877-3p and MiR-149-5p, but the specific mechanisms need further research.
In this study, we have performed a comprehensive transcriptome analysis in D-T2DM patients' peripheral blood, revealing the contribution of epigenetic changes to D-T2DM. However, the present study has several limitations. First, in order to better understand the mechanism of occurrence and development of D-T2DM, it is important to compare D-T2DM with other types of T2DM, such as a milder form of T2DM. This will be the main part of our future study. Next, the major targets of insulin actions are skeletal muscle, liver, and adipose tissue rather than the blood. Therefore, our results can explain the pathogenesis of D-T2DM only partially. We anticipate the transcriptome profile of typical insulin-targeted tissues of D-T2DM patients will be investigated in the future. In addition, because peripheral blood also includes multiple cell populations and bioactive substances, further studies are required to characterize a more precise role of ncRNAs in peripheral blood of D-T2DM patients. Our results indicate that ncRNAs may exert their biological functions in the pathogenesis of D-T2DM by interacting with each other or related proteins genes. However, the mechanism of ncRNAs is complicated; thus, we will further study these predicted ncRNAs and their involved signaling pathways. Finally, it will help to reveal the internal mechanism and therapeutic targets of D-T2DM.
In conclusion, we used sequencing analysis to study the expression profile of RNAs in the peripheral blood of D-T2DM and healthy subjects. Differentially expressed mRNAs, lncRNAs, and circRNAs were screened, and their potential biological functions were predicted by bioinformatics analysis. These results may bring new perspectives on the pathogenesis of D-T2DM and promote the development of new therapeutic approaches targeting RNAs.

DATA AVAILABILITY STATEMENT
All clean data have been submitted to Sequence Read Archive with accession number SRP274496.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of Beijing University of Chinese Medicine (BUCM) (2017BZHYLL0105). The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
SG and GJ designed the experiments. BL and XB wrote the manuscript. PL and TA interpreted data and revised the manuscript. JL, YW, JZha, and XY performed experiments. TA, TW, JZhu, and YH performed data collection and analysis. All authors reviewed the manuscript and agreed to the publication of this article.