Genome-Wide Identification of lncRNA and mRNA for Diagnosing Type 2 Diabetes in Saudi Arabia

Introduction

Taif Governorate is located at an altitude of over 1800 m above sea level and has recently experienced an improvement in the quality of life, reflected in the increase in employment opportunities and tourism activities.¹ At this altitude, oxygen levels are low and atmospheric pressure is decreased.² Living at high altitudes is stressful due to susceptibility to hypoxia, an extreme form of altitude sickness.³ People living at high altitudes often have a strong, long-evolved response to hypoxic conditions; this is evident in indigenous populations that have adapted several molecular, cellular, and systemic responses to tolerate hypoxia at high altitudes.⁴ Various physiological responses, including increased heart and respiratory rates and red blood cell production, exist at the systemic level.⁵ Increased red blood cell mass and hemoglobin content in the blood are thought to be induced by gene regulation.⁶ Metabolic studies have found a shift in expression patterns that can provide an increased energy supply for the cells in the absence of aerobiosis (and exhibit less demand for ATP).⁷ This evidence and other physiological responses constitute examples of altitude adaptation.⁸ Usually, adaptations are considered genetic alterations that cause a particular physiological trait to develop, a phenomenon known as adaptive plasticity.⁹ However, not everyone responds in this way. Some individuals reportedly develop adaptive responses, but others, particularly those with chronic diseases like diabetes, experience complications due to living in such locations.⁸ There has been a significant increase in the prevalence of diabetes in high-altitude populations because of urbanization and rapid changes in diets and lifestyles.^10–13 The global and fast expanding diabetes epidemic is likely to become the primary cause of mortality and disability in the future due to the ageing of the population and lifestyle shifts.¹⁴ The International Diabetes Federation estimates that 450 million people aged 18 and above suffer from diabetes, and this number is expected to increase to approximately 690 million by 2045.¹⁵ Notably, the World Health Organization (WHO) has identified Saudi Arabia as having the second-highest incidence of diabetes in the Middle East and the seventh-highest worldwide. Approximately 10 million people in the country have diabetes or are prediabetic.¹⁶

Type 2 diabetes mellitus (T2DM) is the most prevalent form of diabetes¹⁷ and tends to result from genetic, environmental, immunological, and lifestyle factors.^18,19 T2DM is a progressive, chronic disorder whose symptoms advance over time. T2DM is characterized by low insulin sensitivity and defective insulin secretion. High blood glucose levels may also increase the risk of retinopathy, nephropathy, neuropathy, and cardiomyopathy.¹⁷ Early stages of the illness can go undiagnosed, causing symptoms or complications that are not detected until later stages.²⁰ Approximately half of the people living with diabetes are estimated to be undiagnosed.¹⁵ If individuals can be accurately diagnosed early in the asymptomatic phase of the disease, they may benefit from early interventions, limiting the development of the disease and helping them manage their symptoms more effectively. Thus, there has been an increasing focus on finding reliable, responsive, and easily available diagnostics for diabetes. Family history is a significant risk factor for developing this disease; T2DM has a 4- to 6-fold elevated risk among relatives.²¹ Therefore, collecting the full family history of suspected patients is important. Furthermore, as many changes in insulin-responsive tissues are believed to underlie obesity, insulin resistance, and T2DM, it has become increasingly apparent that genetic and epigenetic markers in the blood can also play crucial roles in their respective pathologies.^22–24 Therefore, new predictive biomarkers that can help diagnose diabetes at an early stage are needed, which may also aid in identifying new therapeutic targets.

Currently, genetic and genomic studies are being conducted for disease prevention and treatment.²⁵ New genetic knowledge must be spread across the wide medical field, and the technical skills needed for disease genetic screening, diagnosis, and prevention should not be confined to research or specialist practice.²⁶ Understanding the genetic basis of diseases requires an understanding of variation across the whole genome to determine overall influence. The current focus of clinical genomics is mainly on protein-coding genes; however, the noncoding genome is far larger than the protein-coding equivalent.²⁷ The noncoding genome encompasses transcriptional, regulatory, and structural information, which needs to be integrated into genome annotations to optimize the use of genomic information in the healthcare system.²⁸ According to genome-wide association studies, most diabetes-related genetic variations do not lie in protein-coding regions, making it difficult to identify functional variants.²⁹ This highlights the importance of identifying and characterizing early noncoding RNA (ncRNA) biomarkers for T2DM management. Over the years, several classes of ncRNAs have been discovered.³⁰ Almost all of these ncRNAs are commonly categorized as small ncRNAs (<200 nucleotides), consisting of microRNAs (miRNAs) and circular RNAs (circRNAs), and large ncRNAs, such as long ncRNAs (lncRNAs).^31–33 lncRNA consists of transcripts with a size range from 200 nucleotides to 100 kilobase pair (kbp).^34,35 lncRNAs are transcribed from either strand and classified as sense exonic lncRNAs, antisense exonic lncRNAs, intronic sense and antisense lncRNAs, and 3′- and 5′-UTR-associated RNAs based on their relationship with the neighboring protein-coding genes.³⁶ lncRNAs generate a complex regulatory network by establishing links with transcription factors, transcriptional co-activators, and repressors, which can influence several aspects of transcription.³⁷ Investigations on the effect of lncRNAs under different clinical and physiological conditions have been conducted.^38–40

lncRNAs are implicated in the regulation of numerous biological reactions associated with health and disease.⁴¹ Research has demonstrated the importance of lncRNAs to inflammation,⁴² and the connection between different mediators of inflammation and T2DM has been determined.^43,44 A cross-sectional cohort study showed that the serum neuregulin-4 level is substantially elevated in patients with T2DM compared to that in healthy controls.⁴⁵ This suggests that neuregulin-4 level may serve as a biomarker for T2DM because euregulin-4 has potential anti-inflammatory properties. Furthermore, several other markers have been studied in T2DM. For example, T2DM complications, such as diabetic renal disease, could be diagnosed based on the uric acid to HDL ratio (UHR) because this ratio is connected to T2DM and inflammation.⁴⁶ In T2DM, the UHR ratio is a robust predictor of metabolic syndrome.⁴⁷ Another study found that uncontrolled hypertension is associated with an increased UHR ratio, which is linked to inflammation⁴⁸ and fatty liver disease.⁴⁹ Although inflammation plays a vital role in the development of T2DM and its related complications, hemogram parameters, including mean platelet volume, were regarded as a new inflammatory biomarker in obese patients with T2DM.⁵⁰

As mentioned above, lnc-RNA is linked to inflammatory conditions and T2DM, as well as its associated conditions such as diabetic kidney disease. Additionally, hypertension, obesity, and fatty liver disease are associated with inflammation, so investigating lnc-RNA in diabetes is rational. However, no research has, to the best of our knowledge, expressly investigated the possible function of certain lncRNAs as diagnostic or management biomarkers for T2DM. In this study, we performed transcriptomic analyses to identify molecular biomarkers that influence patients with T2DM who live in high-altitude areas by analyzing noncoding regions (lncRNA) and protein-coding regions (mRNA) of the genome.

Materials and Methods

Ethics and Study Participants

This study was conducted in accordance with the Declaration of Helsinki. The study procedure was approved by the Taif University Research Ethical Committee, Taif, Saudi Arabia (NO.: 43–220). The aim and nature of the methods to be used in this study were discussed with the participants, and written informed consent was obtained from each participant. Two groups of participants living in the Taif region were enrolled—patients diagnosed with T2DM (five women; age: 27–56 years) and a healthy control group (four women; age: 29–57 years)—between January and March 2022. T2DM diagnoses were based on the 1999 WHO diabetes diagnostic criteria.⁵¹ None of the subjects had received hypoglycemic medication. Exclusion criteria for participants included a history of type 1 diabetes, pregnancy, cancer, and chronic or acute diabetic complications.

Blood Collection and RNA Extraction

Fresh blood (5 mL) was collected from each participant. Thereafter, 1.5 mL of the collected blood sample (with 4000–7000 leukocytes/µL) was processed immediately for total RNA extraction using a QIAamp RNA Blood Mini Kit (Qiagen, Hilden, Germany), following the manufacturer’s protocol. The integrity of the RNA was evaluated with an Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA), and its purity was determined using agarose gel electrophoresis and a NanoDrop2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). RNA samples with an RNA integrity number ≥8.0 were processed further.

Library Construction, Sequencing, and Mapping

Ribosomal RNA (rRNA) was removed from total RNA using an rRNA removal kit (Illumina, San Diego, CA, USA), following the manufacturer’s protocol. A KAPA Stranded RNA-Seq Library Preparation Kit (Illumina) was used to complete the RNA sequencing library following the manufacturer’s protocol. Qubit (Thermo Fisher Scientific) and real-time PCR were used to quantify the constructed library, and a bioanalyzer was used to identify the size distribution. Quantified libraries were sequenced on an Illumina HiSeq 2500 platform (Illumina). The annotation data for the reference genome and gene models were acquired directly from the Ensembl genome browser 106 (https://asia.ensembl.org/index.html). Using hierarchical indexing for spliced alignment of transcripts (HISAT 2; version 2.0.4), clean reads were mapped to the Homo sapiens genome (genome assembly: GRCh38.p13).⁵² Figure 1 illustrates the workflow of this study.

Figure 1 Workflow of lncRNA and mRNAs analysis for patients with T2DM versus healthy controls. Abbreviations: CNCI, coding-noncoding-index; CPC, coding potential calculator; PFAM, Pfam Scan database; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

lncRNA Identification

StringTie software (version 3.3.0) was utilized to assemble each sample’s mapped reads⁵³ and run using the library-type option; all other parameters were left at their default values. Transcripts from all samples were merged using 2/cuffmerge. To find new protein-coding transcripts, the transcripts were examined for signs of protein-coding possibility and conserved sequences. Such transcripts were filtered out, and lncRNA candidates comprised those without coding potential.

Coding Potential Prediction

The coding-noncoding index (CNCI) software tool (version 2) was utilized to profile and differentiate protein-coding and noncoding sequences.⁵⁴ The coding potential calculator algorithm (CPC) was used to assess the quantity and quality of the open reading frame in a transcript and search the sequences against a database of known protein sequences to distinguish between coding and noncoding transcripts. In our study, we gathered functional protein information using the UniProt Knowledgebase (https://www.uniprot.org /UniProtKB) and set the e-value to 1e⁻¹⁰. The Pfam Scan tool (version 1.3) was used to determine the presence of any known protein family domains listed in the Pfam database (release 27; Pfam A and Pfam B).⁵⁵ Transcripts with a Pfam match were excluded in the following step.

lncRNA Function Prediction

A correlation analysis was performed using Pearson’s correlation to assess the possibility of co-expression between lncRNAs and mRNAs. An interaction between a lncRNA and an mRNA was considered significant when Pearson’s correlation value was ≥|0.70| and the P-value was <0.05. Two analyses were conducted on the total correlation matrix to determine and categorize the interactions and potential activities of lncRNAs (cis and trans) regarding their target gene. Cis-regulated genes are protein-coding genes co-expressed with a dysregulated lncRNA and located within 30 kb upstream or downstream of the same gene. Some lncRNAs trans-regulate the central transcription factors to engage specific cellular processes.

Differential Expression Analysis

Ballgown R package (version 2.4.2) was used to identify transcripts differentially expressed between the groups using the data from StringTie.⁵⁶ Among any two groups, transcripts with a P-value <0.05 were classified as differentially expressed transcripts.

Comparative Analysis of mRNA Transcripts Differentially Expressed in T2DM with Genome-Wide Association Studies

To verify the functions of the 84 mRNA transcripts that exhibit differential expression in T2DM, the Type 2 Diabetes Knowledge Portal (https://t2d.hugeamp.org) was utilized. This portal contains a collection of genes that have been linked to T2DM and other glycemic traits, including HOMA-B, HbA1c, and fasting insulin adj BMI through various genome-wide association studies (GWAS).

Enrichment Analysis

GOseq R package (version 1.48.0) was used to implement Gene Ontology [GO; annotates genes to biological processes (BPs), molecular functions (MFs), and cellular components (CCs)] enrichment analysis of the differentially expressed genes (DEGs) or lncRNA target genes. GO terms with a P-value < 0.05 were deemed significantly enriched among DEGs.⁵⁷ The Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/) was used to annotate genes to pathways.⁵⁸

Statistical Analysis

GraphPad Prism (version 10.0.0) was used for statistical analyses. The results are presented as mean ± standard error. For all data, P < 0.05 indicated statistical significance. KOBAS R package (Version 3.0) was used to examine the statistical enrichment of DEGs or lncRNA target genes.⁵⁹

Results

Several metrics, including the total number of reads, number of reads, error rate, number of reads mapped to the genome, and number of spliced and non-spliced reads, were used to evaluate the quality of the transcriptome data. The quality parameter findings between the groups are displayed in Table 1.

Table 1 Quality Parameter Information for Transcriptome Data for Both Patients with T2DM and Healthy Controls

lncRNA Identification and Characterization

To demonstrate the differences in lncRNA profiles between patients with T2DM and healthy controls and to determine diabetes-related lncRNAs, RNA-seq was performed. The CNCI, CPC, and Pfam Scan database (PFAM) were used to exclude protein-coding transcripts and predict lncRNAs. Significantly expressed lncRNAs were identified using the overlapping results of these three approaches. Finally, 1766 lncRNAs were assembled using the three software, of which 582 were novel (Figure 2A). The lncRNAs were categorized based on their genomic location to simplify functional interpretation and undertake extensive analysis; this revealed that 637 (45.57%) of the lncRNAs were sense overlapping, 279 (19.96%) were long intergenic noncoding RNAs (lincRNA), 211 (15.09%) were sense intronic, 208 (14.88%) were antisense, and 63 (4.51%) were others (Figure 2B). The findings indicated that sense overlapping lncRNAs were the most abundant lncRNAs in the T2DM group.

Figure 2 lncRNA transcriptome analysis in the T2DM group compared with the healthy control group. (A) Venn diagram representing predicted lncRNA findings using CNCI, CPC, and PFAM. The sum of the numbers in each large circle reflects the overall number of noncoding transcripts, and the portions of the circle that overlap represent the noncoding transcripts identified by all three methods. (B) A pie chart of lncRNA classification—sense overlapping, lincRNA (long intergenic noncoding RNA), sense intronic, antisense, and other distributions. Abbreviations: CNCI, coding-noncoding-index; CPC, coding potential calculator; PFAM, Pfam Scan database; T2DM, type 2 diabetes.

Co-Expression of lncRNAs and mRNAs in Patients with T2DM vs Healthy Controls

A screen was performed for lncRNAs or mRNAs with significant expression (the default threshold of FPKM score was selected as 1), and the results were analyzed to generate Venn diagrams. The co-expressed lncRNAs and mRNAs were displayed in a Venn diagram separately (Figures 3A and B) to determine the total number of lncRNAs and mRNAs specifically expressed within and between the groups. The co-expression of lncRNAs and mRNAs between T2DM and healthy control groups provides insights into the influence of T2DM on the co-expression pattern. Figure 3A shows that 148 lncRNAs were uniquely expressed in patients with T2DM, 118 in healthy controls, and 191 in both groups. Furthermore, 467 mRNAs were exclusively expressed in patients with T2DM, 654 in healthy controls, and 658 in both groups (Figure 3B).

Figure 3 Venn Diagram of uniquely expressed and co-expressed lncRNAs and mRNAs in the T2DM and healthy control groups. Expression pattern of (A) lncRNAs and (B) mRNAs. Abbreviation: lncRNA, long noncoding RNA.

lncRNA and mRNA Expression Profiles in Patients with T2DM vs Healthy Controls

The relative expression of lncRNAs and mRNAs was analyzed using high-throughput sequencing to explore possible correlations between alterations in lncRNAs and mRNAs and the development of T2DM. The results identified 582 lncRNAs and 2131 mRNAs in the T2DM group. We found that in the T2DM group, 22 lncRNA transcripts were differentially expressed, of which 10 (1.72%) were upregulated, 12 (2.06%) were downregulated, and 560 showed no difference (96.22%). Furthermore, 84 mRNA transcripts were significantly differentially expressed, of which 27 (1.27%) were upregulated, 57 (2.67) were downregulated, and 2048 showed no difference (96.06%). Transcripts were categorized as differentially expressed when the fold change in expression was more than 2.0 and P ˂ 0.05. Volcano plots and pie charts were used to compare the expression profiles of lncRNAs and mRNAs between the T2DM and healthy control groups (Figure 4A–D).

Figure 4 lncRNAs and mRNA expression profiles in T2DM and healthy control groups. Volcano plots clustering analysis of (A) lncRNAs and (B) mRNA. Pie charts represent the percentage of differentially expressed (C) lncRNAs and (D) mRNA. P < 0.05 was considered significant; expression changes are shown in the T2DM group compared with those in the healthy control group. Magenta represents genes whose expression has increased by >2 fold, while green represents genes whose expression has decreased by >2 fold. Abbreviations: T2DM, type 2 diabetes; lncRNA, long noncoding RNA; mRNA, messenger RNA.

We demonstrated that the T2DM group had altered expression of lncRNAs and mRNAs compared with that of the healthy control group. The top 10 (5 upregulated and 5 downregulated) differentially expressed lncRNAs and mRNAs are shown in Table 2 and Table 3. Under varied experimental settings, cluster analysis was performed to identify genes with comparable expression patterns. A hierarchical clustering analysis was performed to identify the expression patterns of differentially expressed lncRNAs (22) and mRNA transcripts (84) in study groups by considering the FPKM. The clustering information from several experiments indicated that genes with the same gene expression patterns might have comparable roles or be involved in the same biological processes (Figure 5A and B). These findings suggest that differentially expressed lncRNAs and mRNAs are associated with T2DM development.

Table 2 Top 10 Differentially Expressed lncRNAs in the T2DM and Healthy Control Groups

Table 3 Top 10 Differentially Expressed mRNAs in the T2DM and Healthy Control Groups

Figure 5 Hierarchical clustering analysis of significant differential expression profiles between T2DM and healthy groups. (A) lncRNAs and (B) mRNAs. Each row is a transcript ID, and each column represents a sample. Upregulation is represented by magenta, whereas downregulation is represented by green. Abbreviations: DM, diabetes mellitus; T2DM, type 2 diabetes; H, healthy; lncRNA, long noncoding RNA; mRNA, messenger RNA.

The 84 mRNA transcripts that exhibited substantial differential expression were further cross-referenced with T2DM GWAS to determine their potential relevance to the genetic underpinnings of the disease. These genes were NUDT22, ATM, IL6R, FMNL1, TANGO2, ACRBP, PTPRJ, SMCHD1, and NUCB2; 3 genes related to HOMA-B-related loci, RNF19B, TNRC18, and KXD; 10 genes related to HbA1c-related loci, ZSWIM1, AKAP13, STK10, ZAP70, LAMTOR4, METRNL, CTAGE5, USP34, MAPKAPK5, and APOBEC3A; and 3 related to fasting insulin adj BMI-related loci, PIK3R5, SETX, and TAF13 (see Table 4).

Table 4 Differentially Expressed mRNAs Cross-Referenced with T2DM GWAS Available at https://t2d.hugeamp.org

lncRNA Function Prediction

Correlation analysis was performed to investigate the possibility of co-expression between lncRNAs and mRNAs and to predict the lncRNA target genes. In predicting cis lncRNAs–mRNA, no differentially expressed lncRNAs could be linked to nearby genes. However, several lncRNAs were identified to regulate their target protein-coding genes in a trans manner. The top 10 differentially expressed lncRNAs identified in this study significantly correlated to 64 nearby genes, as listed in Table 5, with a Pearson’s correlation value ≥|0.70| and P-value <0. 05.

Table 5 Prediction of Top Differentially Expressed lncRNA-Target mRNA Genes via lncRNA–mRNA Co-Expression Trans-Interaction in the T2DM Group Compared with the Healthy Control Group

The lncRNA TCONS_00098523 was linked to 11 genes, namely, RIOK3, ZEB1, PPM1B, ZNF621, LRRFIP1, TCF25, ZNF383, ZNF844, ZNF611, SFPQ, and SIN3B. The lncRNA TCONS_00098587 was linked to six genes, namely, TCF25, ZEB1, LRRFIP1, PPM1B, ZNF844, and TRIM22. The lncRNA TCONS_00060460 was linked to six genes, namely, ZNF621, ZNF383, GTF2H2, RIOK3, SFPQ, and PPM1B. The lncRNA TCONS_00007325 was linked to TCF25. The lncRNA TCONS_00098489 was linked to 13 genes, namely, PPM1B, ZNF844, ZNF383, LRRFIP1, ZEB1, RIOK3, SIN3B, ZNF417, UBE2I, ZNF611, SFPQ, MED6, and TRIM22. The lncRNA TCONS_00059776 was linked to eight genes, namely, KLF2, AKNA, ZNF580, CREBBP, ZNF708, ZNF791, REL, and ZNF841. The lncRNA TCONS_00004761 was linked to ZNF414. The lncRNA TCONS_00098679 was linked to two genes, namely, ZNF101 and ZBTB25. The lncRNA TCONS_00060436 was linked to eight genes, namely, ZNF580, REL, CREBBP, ZNF708, ZNF841, KLF2, AKNA, and ZNF791. The lncRNA TCONS_00029866 was linked to eight genes, namely, AKNA, ZNF708, CREBBP, REL, ZNF580, ZNF791, KLF2, and ZNF841.

Gene Ontology Enrichment Analysis of lncRNA Target Genes and mRNA

GO terms were predicted to determine the function and relationship of differentially expressed lncRNA target genes and mRNAs in the T2DM and healthy groups. The most significant GO analysis results of lncRNA targets and mRNAs are shown in Figure 6. For lncRNA target genes, the enriched MF terms were DNA binding transcription factor activity, DNA binding, and ion binding (Figure 6A). Enriched CC terms were intracellular, organelle, and nucleoplasm (Figure 6B). The most significantly enriched BP terms were cellular nitrogen compound metabolic and biosynthetic processes (Figure 6C). The most significant GO terms of the mRNAs were enriched in MFs (Figure 6D).

Figure 6 Gene Ontology enrichment analysis of differentially expressed lncRNA target genes and mRNAs in the T2DM and healthy groups. (A) Molecular functions (MF), (B) cellular components (CC), and (C) biological processes (BP) of lncRNA target genes. (D) MF, (E) CC, and (F) BP of mRNA. Abbreviations: T2DM, type 2 diabetes; lncRNA, long noncoding RNA; mRNA, messenger RNA.

For mRNA, the most significantly enriched MF term was kinase activity. The other top terms, which were not significant, were ion binding, mRNA binding, and helicase activity (Figure 6D). No significantly enriched CC terms were found, but the gene networks appeared to be involved with the intracellular, lysosome, and organelle terms as the top three terms (Figure 6E). No significantly enriched BP terms were found, but the top three terms were cellular protein modification process, cell motility, and response to stress (Figure 6F).

KEGG Pathway Enrichment Analysis of lncRNA Target Genes and mRNA

Key pathways for lncRNA target genes and mRNA were analyzed through KEGG enrichment. lncRNA target genes were enriched in nine pathways but not significantly (Figure 7A). Notably, we found three lncRNA target genes enriched in six pathways. UBE2I was enriched in the NF-kappa B signaling pathway, ubiquitin-mediated proteolysis, and RNA transport. GTF2H2 was enriched in basal transcription factors and nucleotide excision repair. PPM1B was enriched only in the MAPK signaling pathway.

Figure 7 KEGG pathway analysis of differentially expressed lncRNA targets and mRNAs in T2DM and healthy groups. (A) Upregulated and (B) downregulated KEGG pathways of lncRNA target genes. (C) Upregulated and (D) downregulated KEGG pathways of mRNA. Abbreviations: lncRNA, long noncoding RNA; KEGG, Kyoto Encyclopedia of Genes and Genomes; mRNA, messenger RNA.

Twenty-seven pathways were downregulated, of which only two were significantly downregulated (Figure 7B shows the top 20 pathways). The significantly enriched pathways identified were the FoxO signaling (P = 0.00075) and viral carcinogenesis pathways (P = 0.00172). Based on the results, the affected lncRNA target genes in the FoxO signaling pathway were KLF2 and CREBBP, and those in the viral carcinogenesis pathway were REL and CREBBP. Notably, we found that CREBBP was enriched in the most relevant downregulated pathways, including notch, TGF-beta, glucagon, HIF-1, wnt, and Jak-STAT signaling pathways; long-term potentiation; adherens junction; and cell cycle.

Upregulated mRNA transcripts were enriched in 22 pathways but not significantly (Figure 7C shows the top 20 pathways). The downregulated mRNA transcripts were enriched in 98 pathways, of which 81 were significantly downregulated (Figure 7D shows the top 20 pathways). The related pathways include the Ras signaling pathway (P = 0.0000053), Jak-STAT signaling pathway (P = 0.000028), EGFR tyrosine kinase inhibitor resistance pathway (P = 0.000105), HIF-1 signaling pathway (P = 0.000209), T cell receptor signaling pathway (P = 0.000221), cholinergic synapse (P = 0.000259), natural killer cell-mediated cytotoxicity (P = 0.000454), PI3K-AKT signaling pathway (P = 0.000529), mTOR signaling pathway (P = 0.000661), aldosterone-regulated sodium reabsorption (P = 0.000910), axon guidance (P = 0.000966), chemokine signaling pathway (P = 0.001148), carbohydrate digestion and absorption (P = 0.001245), and type II diabetes mellitus (P = 0.00135).

Figure 8 shows the most likely KEGG pathways linked to downregulated mRNA transcripts involved in T2DM .

Figure 8 Enriched mRNA transcript genes in the KEGG pathways most likely involved in T2DM. The graph was generated using Origin Pro 2023 (OriginLab, Northampton, MA, USA). Abbreviations: T2DM, type 2 diabetes; KEGG, Kyoto Encyclopedia of Genes and Genomes; mRNA, messenger RNA.

Discussion

According to the latest data from WHO, Saudi Arabia is ranked seventh worldwide in the number of individuals diagnosed with diabetes mellitus.⁶⁰ In addition, over the past 3 years, Saudi Arabia has recorded an increase in diabetes mellitus cases, roughly equivalent to a 10-fold increase.⁶¹ The pathogenesis of T2DM is complicated and consists of multiple factors that operate in concert to produce this condition.⁶² Genetic, environmental, immunological, and lifestyle factors typically contribute to developing T2DM.^18,19 Recent research has demonstrated the importance of lncRNA in T2DM.⁶³ The present study utilized genome analysis using RNA sequencing to investigate the expression of lncRNA and mRNA transcripts of female patients with T2DM compared with those of healthy females in Taif City, Saudi Arabia. To our knowledge, this study is the first to be conducted in a high-altitude area, such as Taif City, to evaluate the lncRNA and mRNA expression profiles in females with T2DM to gain a better understanding of the molecular mechanisms behind the etiology of T2DM at high altitudes. In the present study, we identified 1766 lncRNAs in the T2DM group, of which 582 were novel. Additionally, we found that compared with those in the healthy control group, 22 lncRNA transcripts (10 upregulated and 12 downregulated) and 84 mRNA transcripts (27 upregulated and 57 downregulated) were differentially expressed in patients with T2DM, and most of these transcripts were novel. Hierarchical clustering analysis of expression profiles showed significant differences between the T2DM and healthy control groups. The data indicated that this analysis may lead to identifying important target genes implicated in the development of T2DM.

Based on whole-genome sequencing, lncRNA target genes in patients with diabetes were downregulated in two pathways: Forkhead box O (FoxO) signaling and viral carcinogenesis. KLF2 and CREBBP genes were most likely affected in the FoxO signaling pathway, while in the viral carcinogenesis pathway, REL and CREBBP were the most likely affected genes. FOXO is a family of transcription factors, and the FoxO signaling pathway controls many cellular physiological processes, such as glucose metabolism, cell death, cell-cycle regulation, DNA damage repair, resistance to oxidative stress, and adaption to stress stimuli.^64–66 Post-translational modification strictly controls the activity of FOXOs. Patients with diabetes are at an elevated risk of acquiring various severe health complications. Evidence indicates that diabetes-induced activation of FOXO1 is linked to several diabetic problems.⁶⁷ In vivo model knockdown of FOXO1 can help eliminate retinal microvascular endothelial cells that occur in the initial phase of diabetic retinopathy.⁶⁸ In our study, KLF2 (encoding a zinc-finger transcription factor) was the most likely affected gene in the FoxO signaling pathway. According to reports, KLF2 is crucial in preserving endothelial function.⁶⁹ Cell-based investigations have demonstrated that KLF2 directly regulates important endothelial genes, including endothelial nitric oxide synthase (eNOS), thrombomodulin (THBD),^70,71 and genes that encode proteins with anti-thrombotic and anti-inflammatory properties.⁷² KLF2 is inhibited by 30 mM glucose exposure in human umbilical vein endothelial cells.⁷³ KLF2 inhibition by high glucose is a potential diabetic vasculopathy mechanism.⁷⁴ Furthermore, KLF2 is a powerful angiogenesis inhibitor; as shown in an animal angiogenesis model, the overexpression of KLF2 suppresses vascular endothelial growth factor A (VEGFA).⁷⁵ In addition, KLF2 can reduce HIF1-α production and affect its function.⁷⁶ HIF1 is a key transcription factor that regulates metabolic adaptation to hypoxia.⁷⁷ Moreover, HIF1 regulates the promotion of glycolysis and inhibition of mitochondrial respiration, thereby decreasing oxygen uptake and inhibiting the generation of reactive oxygen species.⁷⁸ Under intermittent hypoxic conditions, HIF1 increases the expression of pro-inflammatory and pro-angiogenic genes to induce angiogenesis.⁷⁹ In endothelial cells, the expression of KLF2 was increased under hypoxia, whereas KLF2 knockdown boosted HIF1-α expression.⁸⁰ The results of the present study show that CREBBP most likely plays a role in downregulating the FoxO viral carcinogenesis signaling pathway. CREBBP, a lysine acetyl transferase involved in many signaling pathways, is implicated in controlling the accessibility of chromatin and transcription.⁸¹ Based on our study, CREBBP downregulates the FoxO signaling pathway to reduce diabetes complications. We also found that the viral carcinogenesis pathway is significantly downregulated.⁸² Patients with T2DM are associated with a higher chance of contracting viral infections, as was recently demonstrated during the COVID-19 pandemic.⁸²

We found that the mRNAs significantly downregulated 81 pathways. The most relevant pathways included the Ras, Jak-STAT, PI3K-AKT, mTOR, HIF-1, T cell receptor, and chemokine signaling pathways; cholinergic synapse; natural killer cell-mediated cytotoxicity; aldosterone-regulated sodium reabsorption; axon guidance; carbohydrate digestion and absorption; type II diabetes mellitus pathway; and EGFR tyrosine kinase inhibitor resistance pathway.

The Ras signaling pathway is an essential growth regulator in all eukaryotic organisms.⁸³ The renin–angiotensin system (RAS) is closely associated with the pathogenesis of insulin resistance/diabetes,⁸⁴ and RAS inhibition improves insulin sensitivity in humans.⁸⁵

In our study, PIK3CD and PIK3R5 were enriched in all relevant significantly downregulated pathways. Consistent with our findings, PIK3CD expression was significantly reduced in T2DM in a previous study.⁸⁶ As insulin resistance is frequently identified as the most important contributor to the development of T2DM, insulin resistance might be treated by targeting the PIK3CD gene.⁸⁶ Furthermore, by analyzing the microRNA–mRNA expression patterns and functional network of the submandibular gland in T2DM mice, PIK3CD was surmised to play essential roles in developing diabetes-mediated hyposalivation.⁸⁷ PIK3CB and PIK3CA are among the genes predicted to be predominantly ordered, according to a comprehensive analysis of the functions of highly disordered proteins in T2DM.⁸⁸ These findings elucidated the primary biological functions of these proteins as well as the functional significance of some of their sites, which often play a part in binding between proteins and possess sites for diverse post-translational modifications.⁸⁸ A previous study used high-throughput sequencing to investigate the lncRNA and circular RNA network in T2DM. A protein–protein interaction network was built to identify several hub mRNAs, including PIK3R5, enriched in key pathways such as the mTOR signaling and apoptosis pathways.⁸⁹ In a previous in silico study, bioinformatics analysis was performed to comprehend differential gene expression and patterns and the enriched pathways for obesity and T2DM. Several overexpressed genes that are direct components of the T cell signaling pathway, including PIK3R5, were identified.⁹⁰

In the current study, the IL6R gene was enriched in four relevant pathways, including the Jak-STAT, HIF-1, and PI3K-Akt signaling pathways and EGFR tyrosine kinase inhibitor resistance. Serum levels of the IL6/IL6R are considerably elevated in T2DM;⁹¹ IL6/IL6R has important implications for T2DM. IL6R suppresses pancreatic beta-cell viability and decreases apoptosis-related gene expression to inhibit cell apoptosis by controlling the JAK/STAT signaling pathway via miR22.⁹² IL-6 primarily activates the JAK/STAT signaling pathway but also activates ERK1/2 and PI3K.⁹³ Modifications in JAK/STAT signaling are linked to numerous complications of T2DM.⁹⁴ In the present study, TYK2 was enriched in the Jak-STAT signaling pathway and osteoclast differentiation. Tyk2 is a member of the Janus family kinases (Jaks), which are activated by cytokines, including IL10, IL12, and IL18, and perform important functions in signal transduction.⁹⁵ In mice with gene-targeted knockout of Tyk2 kinase, the function of Tyk2 in the progression of obesity and diabetes was examined. As these animals aged, they developed obesity and T2DM, suggesting that Tyk2 kinase plays a vital role in the progression of these disorders.⁹⁶ Furthermore, a study investigated the association of TYK2 gene polymorphisms with T1DM and T2DM, focusing on the correlation with flu-like syndrome. The results revealed that the variant of the TYK2 promoter has been linked with an increased risk for diabetes, making it an attractive candidate for virus-induced diabetes.⁹⁷

In the current study, ZAP70 was enriched in the Ras and T cell receptor signaling pathways and natural killer cell-mediated cytotoxicity. ZAP70 is a Syk family kinase that plays a key role in triggering the T cell receptor signaling pathway and cell migration and death.⁹⁸ Utilizing gene expression profiles from the Gene Expression Omnibus and a weighted gene correlation network, a comprehensive study was conducted to identify key genes implicated in the development of T2DM-associated cardiovascular disease; the researchers identified 19 genes, including ZAP70, involved in atherosclerosis.⁹⁹ Earlier work combined miRNA and mRNA datasets to identify significant sepsis-related miRNA and mRNA pairings.

In the present study, the LAMTOR4 gene was enriched in the mTOR signaling pathway. mTOR signaling controls development, growth, motility, and protein production, in addition to various cellular and metabolic functions.¹⁰⁰ A study showed that mTOR dysregulation has a significant pathology in the progression of diseases, including T2DM.¹⁰¹ Earlier research emphasized the crucial role of LAMTOR4 as a regulatory element.¹⁰² LAMTOR1 and LAMTOR4 are important in the mTOR signaling pathway. To the best of our knowledge, information on the role of this gene in the development of T2DM at the molecular level is unknown.

In the current study, the SSH2 gene was enriched in axon guidance pathways. These pathways control axon guidance, synaptic development, progenitor movement, and cell migration.¹⁰³ Axon guidance pathways are stimulated in patients with T2DM.^104,105 The profiles and networks of miRNA–mRNA expression in the submandibular gland tissues of an animal model of spontaneous T2DM were described in a previous study, which demonstrated that the 11 differentially expressed microRNAs were related to 820 mRNAs, indicating a link between the miRNAs and mRNAs of their target genes. From these, a network of 11 differentially expressed microRNAs and their target genes was built. According to the network, every miRNA was associated with many mRNAs, and every mRNA was associated with different miRNAs. The mRNA SSH2, for instance, interacts with three miRNAs.⁸⁷ Studies to uncover the correlations between diabetes and sensorineural hearing loss identified two new genes, NOX1 and SSH2.¹⁰⁶

To highlight the origin-specific targets, our results were compared to previously published transcriptomes of T2DM and healthy neutrophils of people of different ethnicity, including 9 Caucasians, 1 Hispanic, and 11 African-Americans, In their investigation, the researchers found a considerable difference in gene expression between individuals with T2DM and those with healthy neutrophils.¹⁰⁷ Specifically, they observed a reduction in gene expression associated with inflammation and lipid metabolism in T2DM, as evidenced by the downregulation of SLC9A4, NECTIN2, and PLPP3. Furthermore, the top KEGG pathways included sphingolipid metabolism, glycerophospholipid metabolism, ether lipid metabolism, Fc gamma R-mediated phagocytosis, and phospholipase D signaling pathway. The top GO terms in the biological processes category included ammonium ion metabolic process and surfactant homeostasis; those associated with molecular functions included sphingosine-1-phosphate-phosphatase activity; and those involved in cellular components included plasma membrane and integral component of plasma membrane.¹⁰⁷

There are some limitations to this study. The small number of samples used for RNA sequencing might have influenced the precision of the results; therefore, it is essential to increase the sample size to validate the results. The results acquired are preliminary and must be verified.

Conclusion

To the best of our knowledge, this comprehensive study is the first to explore the applicability of certain lncRNAs as diagnostic or management biomarkers for T2DM in females in Taif City, Saudi Arabia through the genome-wide identification of lncRNA and mRNA profiling using RNA seq and bioinformatics analysis. This study identified three lncRNA target genes, namely KLF2, CREBBP, and REL. Seven mRNAs, namely PIK3CD, PIK3R5, IL6R, TYK2, ZAP70, LAMTOR4, and SSH2, were significantly enriched in important pathways and perhaps play an important role in the progression of T2DM. Our findings could help in the early diagnosis of T2DM and in designing effective therapeutic targets.