Identification of Two Long Non-Coding RNAs AC010082.1 and AC011443.1 as Biomarkers of Coronary Heart Disease Based on Logistic Stepwise Regression Prediction Model

Coronary heart disease (CHD) is a global health concern with high morbidity and mortality rates. This study aimed to identify the possible long non-coding RNA (lncRNA) biomarkers of CHD. The lncRNA- and mRNA-related data of patients with CHD were downloaded from the Gene Expression Omnibus database (GSE113079). The limma package was used to identify differentially expressed lncRNAs and mRNAs (DElncRNAs and DEmRNAs, respectively). Then, miRcode, TargetScan, miRDB, and miRTarBase databases were used to form the competing endogenous RNA (ceRNA) network. Furthermore, SPSS Modeler 18.0 was used to construct a logistic stepwise regression prediction model for CHD diagnosis based on DElncRNAs. Of the microarray data, 70% was used as a training set and 30% as a test set. Moreover, a validation cohort including 30 patients with CHD and 30 healthy controls was used to verify the hub lncRNA expression through real-time reverse transcription-quantitative PCR (RT-qPCR). A total of 185 DElncRNAs (114 upregulated and 71 downregulated) and 382 DEmRNAs (162 upregulated and 220 downregulated) between CHD and healthy controls were identified from the microarray data. Furthermore, through bioinformatics prediction, a 38 lncRNA-21miRNA-40 mRNA ceRNA network was constructed. Next, by constructing a logistic stepwise regression prediction model for 38 DElncRNAs, we screened two hub lncRNAs AC010082.1 and AC011443.1 (p < 0.05). The sensitivity, specificity, and area under the curve were 98.41%, 100%, and 0.995, respectively, for the training set and 93.33%, 91.67%, and 0.983, respectively, for the test set. We further verified the significant upregulation of AC010082.1 (p < 0.01) and AC011443.1 (p < 0.05) in patients with CHD using RT-qPCR in the validation cohort. Our results suggest that lncRNA AC010082.1 and AC011443.1 are potential biomarkers of CHD. Their pathological mechanism in CHD requires further validation.


INTRODUCTION
Globally, coronary heart disease (CHD) is a major disease with high morbidity and mortality, which seriously threatens human health (Collet et al., 2021), killing approximately 8 million people annually worldwide (Mortality and Causes of Death, 2015). Reportedly, > 700,000 people die from CHD annually in China, and the overall mortality rate has been increasing each year (Wang et al., 2011). CHD is a complex multifactorial disease with many risk factors, such as age, hypertension, diabetes, dyslipidemia, and smoking status (Knuuti et al., 2020). The definitive diagnosis of CHD primarily depends on invasive coronary angiography, which is relatively costly, timeconsuming, and uncomfortable for the patient (Novak et al., 2021). However, CHD in its early stage is not easily detected by routine examinations, such as electrocardiography and cardiac ultrasound, resulting in a high mortality rate of CHD (Miao et al., 2019). Therefore, it is necessary to identify a novel non-invasive biomarker in the early stages of CHD, to enable early diagnosis and prevention of CHD.
Evidence in recent years has showed that long non-coding RNAs (lncRNAs) have a key role in the gene regulation and widely concerned (Statello et al., 2020;Zhang et al., 2020;Hongkuan et al., 2021). LncRNAs are non-coding RNAs with a transcription length >200 bp and lack protein-coding potential. A variety of cellular functions within the nucleus and cytoplasm could be regulated by lncRNAs, which are important for normal development and disease progression (Lu and Thum, 2019). Due to their tissue-specific and condition-specific expression patterns, lncRNAs can be used as biomarkers and therapeutic targets for multiple diseases in blood, plasma, and urine (Beermann et al., 2016). A variety of lncRNAs, such as lncRNA OTTHUMT00000387022 (Cai et al., 2016), AC100865.1 , and uc010yfd.1 (Li et al., 2018), are abnormally expressed in CHD and therefore, can be used as potential biomarkers of CHD. LncRNAs play key roles in specific physiological and pathological processes of CHD, including the induction of vascular smooth muscle cell proliferation, apoptosis, lipid metabolism, and inflammation (Li et al., 2020;Li et al., 2021a). LncRNAs have been considered promising regulatory genes or biomarkers for CHD (Lu and Thum, 2019).
Recently, the proposed competing endogenous RNAs (ceRNAs) hypothesis has provided the possibility for further study of molecular mechanisms underlying several diseased conditions (Cheng et al., 2015;Ouyang et al., 2020;Wang et al., 2020). CeRNAs can act as miRNA sponges and regulate mRNA expression through their miRNA response elements (MREs) (Salmena et al., 2011). A number of studies have shown that lncRNA can regulate the progression of CHD through the ceRNA mechanism. For instance, lncRNA HCG11 affects the expression of FOXF1 by targeting miR-144-3p, thereby regulating the proliferation and apoptosis of vascular smooth muscle cells in atherosclerosis . The lncRNA HOTAIR can protect myocardial infraction and hypoxia-induced cardiomyocyte apoptosis by interacting with miR-519d-3p . In the present study, we aimed to construct the ceRNA regulatory network of CHD using the Gene Expression Omnibus (GEO) database and screen the possible lncRNAs as biomarkers of CHD.

Screening Differentially Expressed lncRNAs (DElncRNAs) and mRNAs (DEmRNAs)
The identification of DElncRNAs and DEmRNAs between patients with CHD and healthy controls were performed by using the "limma" package of R software (version 4.0.1) (R Core Team, 2018). An independent two-sample t-test and the false discovery rate (FDR) adjusted using the Benjamini-Hochberg method were used to analyze the DElncRNAs and DEmRNAs (Long et al., 2019). p < 0.05, FDR <0.05, and |log 2 fold change (FC)| > 1 were considered as the cut-off criteria.

Target Genes Prediction and Construction of the ceRNA Network
The prediction of lncRNAs to target miRNAs was conducted using the miRcode database (Jeggari et al., 2012). DEmRNAs targeted by differentially expressed miRNAs (DEmiRNAs) were retrieved based on the TargetScan, miRDB, and mirtarbase databases, and only those recognized by these three databases were considered candidate DEmRNAs (Agarwal et al., 2015;Wong and Wang, 2015;Chou et al., 2016). Next, the predicted target DEmRNAs and previous identified DEmRNAs were intersected. Based on the above lncRNA-miRNA and miRNA-mRNA interactions, the lncRNA-miRNA-mRNA network was visualized using Cytoscape 3.8.0 software (Shannon et al., 2003). The R package ''pheatmap'' was used to draw heatmaps for DElncRNAs and DEmRNAs. Figure 1 shows the process used for ceRNA network construction.

Functional and Pathway Enrichment Analysis of DEmRNAs
To clarify the potential biological processes of DEmRNAs associated with the ceRNA network, we used Metascape (http://metascape.org) to analyze the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment of the DEmRNAs (Zhou et al., 2019). p < 0.05, min overlap mRNAs≥2, and enrichment factor >1.5 were considered to be statistically significant (Gao et al., 2020).

Constructing a Logistic Stepwise Regression Prediction Model and Screening Hub lncRNAs
The 141 samples from the GSE113079 dataset were randomly divided into training (70% of the samples, 63 patients with CHD and 36 healthy controls) and test (30%, 30 patients with CHD and 12 healthy controls) sets, and they maintained the similar ratio for patients with CHD and healthy controls (p 0.372). SPSS Modeler 18.0 (IBM Canada Ltd., Markham, Ontario, Canada) was used to construct a logistic stepwise regression prediction model for CHD diagnosis based on DElncRNAs. Further, the sensitivity and specificity of the model were calculated using SPSS Modeler. The receiver operating characteristic (ROC) curve was used to evaluate the effectiveness of the classification model, and the area under the curve (AUC) value was calculated from the ROC curve. ROC curves were calculated using the R package "pROC." Principal component analysis (PCA) using the R package "ggplot2" was performed on the screened hub lncRNAs, and their expression levels in patients with CHD and healthy controls were analyzed. The selection process for screening hub lncRNAs is shown in Figure 2.

Study Population
A total of 30 patients with CHD (17 males) and 30 healthy controls (6 males) were recruited in Guang'anmen Hospital, Beijing, China, as the validation cohort. This study complied with the Declaration of Helsinki and was approved by the Ethics Committee of Guang'anmen Hospital, China Academy of Chinese Medical Sciences (no. 2019-224-KY). All subjects were aged between 42 and 77 years. CHD diagnosis was based on the European Society of Cardiology criteria, which specifies at least one vessel lesion (>50% narrowing of luminal diameter) using coronary angiography (Task

Total RNA Isolation and Reverse Transcription
Peripheral blood samples (4-8 ml) were collected in the early morning from patients with CHD and healthy controls using EDTA vacuum anticoagulant blood vessels. Peripheral blood mononuclear cells (PBMCs) were isolated using the Red Blood Cell Lysis kit (TIANGEN, Beijing, China) and total RNA was extracted from PBMCs using TRIzol reagent (Invitrogen, Carlsbad, CA, United States) according to the manufacturers' instructions. RNA quality and quantity were evaluated using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, MA, United States). LncRNA reverse transcription was performed using the lncRcute lncRNA First-Strand cDNA Kit (TIANGEN) from 1 µg total RNA according to the manufacturer's instructions.

Construction of the ceRNA Network
We used the miRcode database to predict miRNA-targeted lncRNAs.  (Figures 4, 5).

Functional Enrichment Analysis
To determine the biological functions and pathways of the 40 DEmRNAs in the constructed ceRNA network, we used the Metascape database to analyze the GO functional enrichment and KEGG pathway enrichment ( Figure 6). The top terms of biological processes were "BMP signaling pathway," "protein polyubiquitination," and "cellular response to transforming growth factor beta stimulus"; the top cellular component terms were "spindle pole," "mitotic spindle pole," and "catenin complex"; the top molecular function terms were "activin binding," "ubiquitinprotein transferase activity," and "transforming growth factor beta receptor binding." According to KEGG pathway analyses, the most significant pathways included "TGF-beta signaling pathway," "TNF signaling pathway," and "leukocyte transendothelial migration."

Logistic Stepwise Regression Prediction Model Based on DElncRNAs
By setting 38 DElncRNAs as independent variables, we constructed a logistic stepwise regression prediction model for CHD. The analysis reveals that AC010082.1 and AC011443.1 (p < 0.05) as lncRNAs were statistically significant for model construction among the 38 independent variables. The prediction model formula is: Logit(P) −13.718 + AC010082.1*8.099 + AC011443.1*3.825. The sensitivity, specificity, and AUC were 98.41%, 100%, and 0.995, respectively, for the training set, and 93.33%, 91.67%, and 0.983, respectively, for the test set ( Figure 7). These two lncRNAs can be regarded as important signatures for the prediction of CHD. PCA of these two lncRNAs can distinguish patients with CHD from healthy controls (Figure 8). At the same time, the results show that AC010082.1 and AC011443.1 are upregulated in patients with CHD (p < 0.05), and that AC010082.1 and AC011443.1 have a significant positive correlation (r 0.7, p < 0.01, Figure 9).

DISCUSSION
Despite considerable advances in modern medicine, CHD diagnosis remains difficult, particularly in the early stages.
LncRNA is an important regulator of gene expression, and one of its regulatory mechanisms is ceRNA. MiRNAs can bind to mRNA through MREs, resulting in mRNA degradation or translation inhibition. When lncRNAs and mRNAs have the same MREs, lncRNAs can compete with mRNA for binding to miRNA to prevent mRNA degradation and achieve indirect regulation of mRNA expression level . The ceRNA mechanism provides a new method for studying the biological functions of lncRNAs and has garnered extensive attention to date (Li et al., 2021b;Lin et al., 2021;Ma et al., 2021). In this study, a ceRNA network consisting of 38 lncRNAs, 21 miRNAs, and 40 mRNAs was constructed using bioinformatics analysis. KEGG function enrichment in this network was mainly   Frontiers in Genetics | www.frontiersin.org November 2021 | Volume 12 | Article 780431 6 concentrated in the TGF-β signaling pathway, TNF signaling pathway, and leukocyte transendothelial migration. These signaling pathways have been previously shown to be closely related to the pathological mechanisms of CHD (Shen et al., 2020;Guo et al., 2021;Sluiter et al., 2021). The TGF-β signaling pathway can regulate the proliferation and migration of vascular smooth muscle and endothelial cells, and affect the stability of plaque, which plays an important role in CHD pathogenesis (Zeng et al., 2016). The TNF signaling pathway is involved in myocardial I/R injury and cardiomyocyte apoptosis (Tan et al., 2018). Leukocyte transmigration is critical to the inflammatory response and accelerates the progression of atherosclerosis (Muller, 2011;Bhui and Hayenga, 2017). These DElncRNAs can interfere with pathological changes in CHD through various signaling pathways.
However, further lncRNA analysis using the logistic stepwise regression prediction model revealed that AC010082.1 and AC011443.1 might be important biomarkers of CHD, which was confirmed using RT-qPCR. AC010082.1 is an antisense molecule with a length of 546 bp and is located on chr7: 18,429,062-18,430,738, whereas AC011443.1 is also an antisense molecule with a length of 365 bp and is located on chr19:39,134,882-39,136,463. Unfortunately, no relevant studies have been conducted on these two lncRNAs, and their direct mechanism of action remains unclear. Through bioinformatics  analysis, it was found that these two DElncRNAs affect the expression of 25 mRNAs by regulating eight miRNAs. The hsa-miR-17-5p, hsa-miR-20b-5p and hsa-miR-27a-3p were jointly regulated by these two DElncRNAs, which may contain a common MRE and have a regulatory relationship through ceRNA mechanism. We found that AC010082.1 and AC011443.1 were overexpressed in patients with CHD and showed a significant positive correlation. The mechanism may be that when lncRNA AC010082.1 is highly expressed in patients with CHD, it will bind more miRNAs, so as to release the inhibition of miRNA on the other lncRNA AC011443.1, thus leading to its high expression (Salmena et al., 2011). AC010082.1 and AC011443.1 bind to miRNA through MREs, thus affecting mRNA expression and biological functions. Interestingly, we found that hsa-miR-17-5p, hsa-miR-20b-5p and hsa-miR-27a-3p co-regulated by these two DElncRNAs have important biological roles in CHD. The expression of hsa-miR-17-5p is down-regulation in patients with CHD, which could be as a biomarker of CHD (Zhelankin et al., 2021) and a good predictor of the severity of coronary atherosclerosis . In addition, lncRNA SNHG16 could regulate NF-κB signaling pathway by binding hsa-miR-17-5p, and then affect proliferation and inflammatory response in atherosclerosis patients, which provide a potential target for treating AS . Furthermore, miR-20b-5p plays an important role in mediating cardiomyocytes apoptosis via the HIF-1 α/NF-κ B pathway (Zhou et al., 2020). And circHIPK3 acting as miR-20b-5p sponge could accelerate cardiomyocyte autophagy and apoptosis during myocardial I/R injury . Reportedly, miR-27a expression is upregulated in patients with CHD, and its expression level is remarkable correlated with the severity of coronary artery stenosis, which can be used as a biomarker of CHD (Babaee et al., 2020). Meanwhile, miR-27a may lead to the progression of atherosclerotic plaques by downregulating ABCA1 and ABCG1 gene expression (Rafiei et al., 2021).
Besides, these two DElncRNAs affected TGFBR3 expression by co-regulating the miR-27a-3p sponge via bioinformatics prediction. TGFBR3 expression is significantly upregulated in patients after myocardial infarction and may serve as a therapeutic target and biomarker for myocardial damage by activating p38 MAPK to induce cardiomyocyte apoptosis (Chu et al., 2011;Chen et al., 2019). Moreover, it has been found that lncRNA SOX2-OT exacerbates hypoxia-induced cardiomyocyte injury by regulating the miR-27a-3p/TGFβR1 axis (Yang and Lin, 2020). Therefore, we speculated that AC010082.1 and  AC011443.1 play important roles in the pathological progression of CHD through the ceRNA mechanism and are potential biomarkers of CHD. The present study had some limitations. First, we only verified the possibility of using AC010082.1 and AC011443.1 as biomarkers through RT-qPCR, and the clinical samples were obtained from only one research center. In future, larger sample sizes and data from multiple centers countrywide are required to verify our findings. In addition, in vitro and in vivo experiments to further explore the functions of AC010082.1 and AC011443.1, and the related pathological mechanisms, are warranted.
In summary, our study identified AC010082.1 and AC011443.1 as possible biomarkers of CHD through bioinformatics analysis and RT-qPCR. The ceRNA network associated with these lncRNAs possibly plays a key role in the pathological process of CHD.

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 in the article/ Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by The Ethics Committee of Guang'anmen Hospital, China Academy of Chinese Medical Sciences (no. 2019-224-KY). The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
CL and JW contributed to the study design. CL and LL searched for and downloaded the gene expression profiles from the Gene Expression Omnibus database. JG was the clinician involved in enrolling the patients and in diagnosing the coronary heart disease cases. CL, LL, and YL performed the RT-qPCR. CL and YL prepared the draft of the manuscript. All authors revised and approved the final manuscript.