Profiling the lncRNA-miRNA-mRNA interaction network in the submandibular gland of diabetic mice

Background Hyposalivation is one of the common symptoms of diabetes. Although long non-coding RNAs (lncRNAs) have recently been reported to play important roles in the pathogenesis of diabetes, the role of lncRNAs in diabetes-induced hyposalivation remains unknown. Methods The present study aimed to explore the function of lncRNA-microRNA-mRNA regulatory network in the submandibular gland (SMGs) under the context of diabetes. LncRNA expression profile of the SMGs was analyzed using microarray technology. Differentially expressed lncRNAs were confirmed using real-time quantitative PCR. Bioinformatics analyses were performed, and Coding-non-coding gene co-expression (CNC) and competing endogenous RNA (ceRNA) networks were constructed to explore the potential mechanisms of diabetes-induced hyposalivation. Results A total of 1273 differentially expressed lncRNAs (536 up-regulated and 737 downregulated) were identified in the SMGs tissues of db/db mice. CNC and ceRNA network analyses were performed based on five differentially expressed lncRNAs validated by real-time quantitative PCR. Gene Ontology analysis of target genes of CNC network revealed that “calcium ion binding” was a highly enriched molecular function. Kyoto Encyclopedia of Genes and Genomes pathway analysis of target genes of ceRNA network revealed that the “mammalian target of rapamycin signaling pathway” was significantly enriched. Conclusions On the whole, the findings of the present study may provide insight into the possible mechanism of diabetes-induced hyposalivation. Supplementary Information The online version contains supplementary material available at 10.1186/s12902-022-01019-1.


Introduction
Diabetes is a group of metabolic diseases characterized by abnormal insulin secretion and/or insulin resistance. In addition to causing long-term damage to different organs, diabetes can also change the function of salivary glands and may cause changes in the composition and volume of salivary secretion, thus affecting the oral homeostasis [1]. Some studies have found that the salivary flow of diabetic patients is significantly decreased compared with control group [2,3]. Clinical studies showed that 92.5% of the elderly patients with type 2 diabetes suffered decreased salivary flow rate [4]. Saliva played a key role in digestion, taste, cleaning, hydration of oral mucosa and tooth protection and was essential for maintaining the dynamic balance of the oral environment [5]. Damage to the salivary glands was typically manifested as a reduction in salivary flow, which can be converted into bad subjective feelings such dry mouth, taste disorders, difficulty with swallowing and chewing, and increased risk of caries. These feelings ultimately have a negative impact on quality of life. Many studies have demonstrated that the incidence of oral fungal and bacterial infection, lichen planus and caries was high in patients with diabetes [2,3,6]. In our previous study, we observed atrophy of the acini, and decreased stimulated salivary flow of SMG in db/db mice [7]. However, the mechanism by which diabetes induces SMG damage is not clear.
LncRNA is a class of RNA molecules that are greater than 200 nucleotides in length. In recent decades, lncRNAs have been identified as important epigenetic factors that regulate various human diseases, such as cancer, neurodegenerative diseases and cardiovascular diseases [8][9][10][11][12]. Recent studies have shown that the incidence and development of type 2 diabetes was associated with lncRNAs. Downregulation of the lncRNA ANRIL improved the cardiac function index and inflammatory factor expression in diabetic rats, enhanced the pathological state of myocardial tissue and myocardial remodeling, and reduced the area of myocardial collagen deposition [13]. The lncRNA NONRATT021972 was upregulated in diabetic rat nervous system cells, suggesting that it may participate in the pathophysiological processes associated with sympathetic neurons in diabetes [14]. Downregulation of the lncRNA SOX2OT regulated the NRF2/HO-1 signaling pathway and protected against high glucoseinduced damage in retinal ganglion cells in vivo and in vitro [15]. Based on these findings, we hypothesized that lncRNAs represented a potential cause of diabetes-induced reductions in salivary secretion. CeRNA can regulate each other at post-transcription level by competing for common miRNAs [16]. LncRNA MEG3 acted as a ceRNA of miR-214 to promote hepatic insulin resistance by facilitating the expression of ATF4 [17]. Protein-coding mRNA DKK1 and PTEN served as ceRNA, showing crosstalk and affecting the expression of each other via competition for miRNAs binding [18]. In studies exploring lncRNAs, the identification of a target gene and its correlation with diabetes-related xerostomia can help us to understand the pathology of diabetes-related xerostomia and provide further information for the treatment of diabetes-related xerostomia.
To better understand the function of lncRNA in the impairment of salivary function in diabetes, this study characterized the expression profiles of lncRNA in the SMG tissues of db/db and db/m mice using microarray. The establishment of CNC and ceRNA networks based on lncRNA may help to explore the potential functions of lncRNAs in diabetes-induced SMG dysfunction.

Animal models
Sixteen-week-old male db/db mice (a model of spontaneous type 2 diabetes mellitus) were purchased from Changzhou Cavens Laboratory Animal (Changzhou, China). Studies have found that estrogen can reduce insulin resistance, improve the ability to deal with glucose, so as to reduce blood glucose. So we used male mice in our research to avoid the effects of estrogen [19]. Age and weight matched db/m mice were used as control group. All experimental procedures were performed in accordance with the Guidelines for the Care and Use of Laboratory Animals (NIH Publication No. 85-23;revised 1996). The study was carried out in compliance with the ARRIVE guidelines as well as the relevant national laws on the protection of animals. All mice were maintained at room temperature with a relative humidity of ~ 60% environment with a 12 h light-dark cycle and fed normal food and water. After anesthesia and sacrifice with CO 2 , the SMG was removed immediately from the animal, frozen in liquid nitrogen for 1 min, and then stored at − 80 °C.

RNA extraction and microarray analysis
Four mice were andomly chosen from each group. The total RNA was extracted from the SMG tissues using the TRIzol LS Reagent (Invitrogen, CA, USA), and was used for lncRNA microarray analysis (KangCheng, Shanghai, China). The concentration and purity of RNA were measured using a NanoDrop ND-1000 (Thermo, MA, USA), and quantity and integrity of RNA was tested using denaturing agarose gel electrophoresis. Microarray hybridization was performed using Quick Amp Labeling Kit, One-Color (Agilent, CA, USA) based on the manufacturer's standard protocols. Agilent Feature Extraction software (Agilent, CA, USA) was used to analyze acquired array images. Quantile normalization and subsequent data processing were performed using the GeneSpring GX v12.1 software package (Agilent, CA, USA). All microarray hybridization and analyses were performed by KangChen Biotech (Shanghai, China). The accession number of microarray data was GSE141411.

Real-time quantitative polymerase chain reaction (RT-qPCR)
We selected the 11 most significantly differentially expressed lncRNAs for RT-qPCR validation. Briefly, total RNA was extracted from two groups of SMG tissues using TRIzol reagent (Invitrogen, CA, USA). Total RNA (1 μg) was reverse-transcribed into cDNA using 5X All-In-One RT MasterMix (ABM, Richmond, BC, Canada) according to the manufacturer's instructions. The RT-qPCR was performed on ViiA ™ 7 Real-Time PCR System (Thermo, MA, USA) using a 2X SYBR Green qPCR Master Mix (Bimake, HOU, USA). The RT-qPCR primers were synthesized by Sangon Biotech (Shanghai, China). ΔΔCq method was used.

Bioinformatics analysis
Volcano plot filtering was used to identify differentially expressed lncRNAs based on the threshold defined as fold-change > 2.0 (Student's t test P < 0.05). The differentially expressed lncRNAs between db/db and db/m mice were displayed by hierarchical clustering. GO analysis (http:// www. geneo ntolo gy. org/) and KEGG analysis (http:// www. genome. jp/ kegg/) were used to explore the roles of the target genes of differentially expressed lncR-NAs [20][21][22].

Coding-noncoding co-expression and competing endogenous RNA networks
The normalized signal intensities of differentially expressed mRNA from whole-genome expression profiling (https:// www. ncbi. nlm. nih. gov/ geo/ query/ acc. cgi? acc= GSE14 1412) and differentially expressed lncRNA validated by RT-qPCR were used to construct CNC networks. In this study, lncRNAs, mRNAs, and miRNAs were collected from the same batch of samples. The lncRNA-mRNA pairs were identified based on a Pearson's correlation coefficient of more than 0.95. The CNC network was constructed using Cytoscape software (The Cytoscape Consortium, CA, USA). Next, we constructed the lncRNA-miRNA-mRNA ceRNA networks based on microarray data. The lncRNA-miRNA interactions and miRNA-mRNA pairs were predicted using Arraystar's home-made miRNA target prediction software and Tar-getScan & miRanda [23][24][25][26].

Statistical analysis
Statistical analysis was performed by GraphPad Prism 5.0 (GraphPad Software, USA). Data of relative expression level of lncRNAs were expressed as the mean ± standard error. Student's t-test was used to test whether the two sets of data were statistically significant, and P < 0.05 was considered to indicate statistical significance of gene expression difference between the two group.

Expression profiling of lncRNAs in the SMG in db/db mice
In our study, lncRNA expression profile of the SMG tissues was carried out using microarray analysis. In comparison with db/m mice, 1273 differentially expressed lncRNAs were identified in the SMG tissues from db/db mice, including 536 upregulated and 737 downregulated lncRNAs (fold-change > 2.0, P < 0.05). Table 1 provided the top ten upregulated and downregulated lncRNAs. A volcano plot of lncRNA expression profiles was displayed in Fig. 1A, the heat maps of the 60 most significantly differentially expressed lncRNAs (30 upregulated and 30 downregulated) was displayed in Fig. 1B.

Validation of lncRNA microarray results
Because there is no previous research that has explored the relationship between lncRNA and saliva secretion in diabetes, there is no literature that can be referred to. Considering that the most differentially expressed lncRNAs should play an important role in the process, we chose the top 10 (up and down regulated) differentially expressed lncRNAs in our study [27]. In these 20 lncRNAs, 9 lncRNAs could not be verified by qRT-PCR because primers could not be designed. Finally, 11 lncRNAs (NR-040589, ENSMUST00000142612, ENSMUST00000139794, ENSMUST00000163495, ENSMUST00000137025, AK021108, NR-045306, ENSMUST00000156336, ENSMUST0000014113, ENSMUST00000120706 and ENSMUST0000069768) were selected from microarray results to carry out RT-qPCR. Of these, five lncRNAs, including 2 downregulated (NR_040589 and ENSMUST00000142612) and 3 upregulated (ENSMUST00000139794, ENS-MUST00000163495, and ENSMUST00000137025), were consistent with microarray analysis results (Fig. 2). The advantages of microarray analysis were fast and efficient, but it also has false positives and false negatives. It is a common phenomenon that qRT-PCR and microarray results are inconsistent. Usually, we selected lncRNA with consistent result for in-depth research [28]. The corresponding primers used were listed in Table 2.

Coding-noncoding co-expression network analysis
Five lncRNAs validated by RT-qPCR and the differentially expressed mRNA from the database (GSE141412) was used to construct the CNC network shown in Fig. 3. There were 615 lncRNA-mRNA connection pairs in the network, including 369 positive connection pairs and 246 negative pairs. The top ten positive and negative interaction pairs according to the Pearson's correlation coefficient were showed in Table 3. These closely interaction pairs may be involved in the regulation of saliva secretion of SMG in db/db mice. GO analysis based on the target genes of the CNC network showed that the mostly enriched biological process was "Cellular response to hormone stimulus", the mostly enriched cellular component was "Extracellular region", while The expression profiling changes of significantly differentially expressed lncRNAs in db/db and db/m mice. A Volcano plots presenting differences in the expression of lncRNAs between db/db mice and db/m mice. B Heat maps showing the top 60-fold change distinct mRNAs expression profiles between the db/db mice and db/m mice mostly enriched molecular function was "Calcium ion binding", as shown in Fig. 4A. KEGG pathway analysis revealed that the mainly enriched pathways were "Cysteine and methionine metabolism", "Phosphatidylinositol signaling system", and "Vitamin digestion and absorption", as shown in Fig. 4B. Cystathionine-gamma-lyase (CTH), glutamate-cysteine ligase modifier subunit (Gclm) and phosphoserine   ENSMUST00000137025  TGC ACA ACT ATG TGA GGA GCA TGG  GGA GGC TCA GGC CAT TCT TCA TTC   ENSMUST00000163495  GCA ATG AGG AGA GTG GAG CAACC  GTC TGC CAA GAG TCC GAA GTA ACC   ENSMUST00000139794  GAC TGG CAG GAA GAA CGG ATGAC  GGC ACA CCT GGT TCC ATC TTA GAC   NR_040589  TCC CGC TTC TCT CGC TTC TCTC  CTG TTG CTC GCC TTC CTG CTG   AK021108  GGT GGC TGT GAA GAC GCT GAAG  TTG ACA ATG TGC TCG TGC TGGAG   ENSMUST00000120706  AGG TGG CTG GTT TCT CTG GA  AGC AGC TAC TCT TGC CTG GA   NR_045306  GCC TGC CTG CCT CCT CTC TC  TTC TGA GCC TGG GTG TCT CTGG   ENSMUST00000142612  TTG AGC ATG TGC CAG ACA GAA GTG  TCT CCA TCG TGT TGC CTT CAA GTG   ENSMUST00000069768  ACG GAC ACC TCT CCT GCC TTG  CTT CCT GAA GCT CGC GGA GAATC   ENSMUST00000156336  ATG AGC GAC TGA GCA GGA CTACC  GGC TGG CGT CCT CAC CTC TAG   ENSMUST00000141103  AAC TGA GCA GAG ACA GCG TGTG  TCA TGG TCT TGT CGT GAG ATG ACG aminotransferase 1 (PSAT1) were enriched in "Cysteine and methionine metabolism". These enriched GO terms and pathways may participate in the DIO-induced hyposalivation.

Competing endogenous RNA network analysis
We constructed ceRNA network based on five lncRNAs validated using RT-qPCR and all differentially expressed mRNAs in the microarray (Fig. 5). We totally found 37  Table S1). GO analysis based on the target genes of ceRNA analysis results showed that the mostly enriched biological process was "Localization", the mostly enriched cellular component was "Cell part", and the mostly enriched molecular function was "Protein binding" (Fig. 6A). KEGG pathway analysis revealed that the main pathways included "Rap1 signaling pathway", "Renal cell carcinoma" and the "mTOR signaling pathway" (Fig. 6B). Signal-induced proliferation-associated 1 like 3 (SIPA1L3), enabled homolog (ENAH), thrombospondin 1 (Thbs1) et al. were enriched in "Rap1 signaling pathway". These results showed that biological processes and regulatory pathways might play vital roles in the secretion of SMG in obesity.

Discussion
In humans, approximately 60% of resting saliva and 40% of stimulated saliva is secreted from the SMG [29]. However, studies focusing on the SMG in the context of diabetes are very limited. We collected SMG of mice for high-throughput sequencing to study the underlying mechanism by which diabetes reduced saliva secretion. Compared with db/m mice, 536 upregulated and 737 downregulated lncRNAs were identified, these changed lncRNAs may play an important role in the diabetesinduced hyposalivation.
We performed RT-qPCR verification on 11 differentially expressed lncRNAs with the most significant dysregulations, and five differentially expressed lncRNAs exhibited consistent with the results of high-throughput sequencing. We previously determined the whole genome expression profile of SMG tissues in db/db mouse and found that 1146 mRNAs exhibited significant dysregulation. Of these, 606 mRNAs were upregulated and 540 mRNAs were downregulated, and this data should be found in GSE141412. We used these five lncR-NAs to perform CNC network and ceRNA network analyses in combination with 1146 differentially expressed mRNAs.
GO analysis of mRNA obtained from the CNC network showed that "Calcium ion binding" is highly enriched in the molecular function. The mobilization of Ca 2+ played an important role in salivary secretion, the activation of muscarinic cholinergic receptors rapidly triggered the release of intracellular Ca 2+ from the endoplasmic reticulum and subsequently the influx of Ca 2+ from the extracellular medium, resulting in a sustained increase in intracellular Ca 2+ [30,31]. The increased intracellular Ca 2+ induced the transport of aquaporin 5, leading to the formation of water pores and thus promoting the rapid increase in transcellular water permeability [32]. In addition, our previous study found that adiponectin can also induce salivary secretion of the rat db/db mouse  by activating adenosine monophosphate-activated protein kinase, and the Ca 2+ signaling pathway played an important role in this process [33]. In human and rabbit SMGs, the activation of muscarinic cholinergic receptors and transient receptor potential vanilloid subtype 1 increased salivary secretion via increased intracellular Ca 2+ [34,35]. In human transplanted epiphora SMG, the elevated intracellular Ca 2+ mobilization induced by muscarinic acetylcholine receptors activation contributed to hypersecretion [34]. In spontaneously hypertensive rats, the damaged Ca 2+ response to carbachol was confirmed in the acinar cells of spontaneously hypertensive rats, which may also be related to the reduced salivary secretion caused by hypertension [36]. These studies indicate that the increased intracellular Ca 2+ derived from extracellular and intracellular Ca 2+ plays an important role in the salivary secretion of SMG. Our results from microarray analysis show that significantly altered "calcium ion binding", suggesting that the salivary secretion of the SMG may be also affected by the Ca 2+ signaling pathway during diabetes. However, the regulatory mechanism of related lncRNA-mRNA requires further research.
KEGG analysis from the CNC network showed that "Cysteine and methionine metabolism" is the highest enriched pathway. An important feature of diabetes is abnormal nutrient homeostasis. One of metabolic abnormalities is to change the metabolism of sulfur-containing amino acid cysteine. Cysteine produced hydrogen sulfide   through enzymatic decomposition, which is a gas transmitter regulating glucose and lipid homeostasis [37]. Studies found that type 2 diabetes was associated with lower plasma cysteine levels [38]. Both cysteine and hydrogen sulfide can inhibit β Cells release insulin after glucose stimulation. The level of low hydrogen sulfide found in plasma of patients with type 2 diabetes may be a compensatory response. In addition to regulating insulin secretion, hydrogen sulfide may regulate other aspects of cell physiology include viability and respiratory function of mitochondria. In the stretptozotocin-induced diabetes model, elevated hydrogen sulfide decrease the viability of cultured β Cells [39]. Our previous research found that high glucose (HG) induced mitochondrial dysfunction and PINK1/Parkin-mediated mitophagy in cultivated rat SMG cell line SMG-C6 cells [40]. We also demonstrated that autophagy plays a crucial role in HG-modulated aquaporin 5 degradation, which contributes to the dysfunction of diabetic SMG [7]. However, the mechanism of diabetes inducing mitophagy and mitochondrial dysfunction in salivary gland is still not well understood. Our results from microarray analysis suggested that the salivary secretion of the SMG may be also affected by "Cysteine and methionine metabolism" during diabetes. From these results, we speculated that these 5 validated lncRNAs might the major regulatory molecules.
LncRNA regulates mRNA via various mechanisms, one of which is ceRNA-mediated changes in the expression of downstream molecules regulated by miRNAs. We performed ceRNA network analysis on five significantly altered lncRNAs and 1146 dysregulated mRNAs to identify the related pathways regulated by the miRNA pathway. We also analyzed the obtained mRNAs via GO and KEGG analyses. KEGG analysis suggested that the "mTOR signaling pathway" was significantly enriched, suggesting that these five lncRNAs might affect the downstream "mTOR signaling pathway" through ceRNA. The PI3K/Akt /mTOR pathway is an intracellular signaling pathway that played a key role in regulating cell cyclemediated processes, including cellular quiescence and cell proliferation [41], and various disease such as epithelial ovarian cancer [42]. Many studies have demonstrated that the "mTOR signaling pathway" plays an important role in the pathophysiological process of diabetes. In MC3T3-E1 cells, high glucose increased the production of reactive oxygen species and induced autophagy by inhibiting the Akt/mTOR pathway [43]. The mTOR/ PI3K/Akt pathway was involved in the regulation of autophagy in diabetic kidney [44]. PI3K/Akt/mTOR pathway was significantly downregulated in the brain of streptozotocin-induced type 2 diabetic rats, this might explain the neurodegeneration commonly observed in diabetes [45]. In addition, the mTOR pathway is also related to the process of salivary secretion. Bleomycin induced epithelial to mesenchymal transformation of human SMG cells via the Akt/mTOR pathway [46]. In our previous studies, we found that autophagy induced aquaporin 5 degradation through the PI3K/Akt /mTOR pathway signaling pathway, resulting in decreased salivary secretion of SMG in db/db mice [7]. These studies suggest that the "mTOR signaling pathway" may play an important role in SMG injury caused by diabetes. The five lncRNAs we verified may involve in the regulation of the "mTOR signaling pathway" through the ceRNA mechanism via miRNA sponging, but further research is needed.
KEGG analysis from the ceRNA network also showed that "Rap1 signaling pathway" is the highest enriched pathway. Ras-associated protein1 (Rap1) is a member of the Ras-like small GTP binding protein family. In endothelium Rap1 is a key positive regulator of angiogenesis and an important regulator of endothelial barrier [47]. Constitutive activation of Rap1 has not only led to T cell anergy, but also inhibited autophagy and supported cancer progression through various oncogenic events [48]. In addition, Rap1 signaling has been shown to play a role in Ca 2+ signaling pathways to control the important cellular functions. Rap1 modulates cardiac Ca 2+ homeostasis through the regulation of Ca 2+ -induced Ca 2+ release in adult ventricular cardiomyocytes [49]. Rap1 may act directly to open Ca 2+ -sensitive K + channels, inducing smooth muscle hyperpolarization and leading to vasorelaxation. Rap1 plays an important role in maintaining normal vascular contractile state and contributes to blood pressure regulation by altering Ca 2+ sensitivity of vascular smooth muscle cells [50]. As intracellular Ca 2+ plays an important role in the salivary secretion of SMG, these 5 validated DE lncRNAs might be a novel treatment strategy for future treatment of diabetesinduced xerostomia in the clinic.
It is worth noting that both KEGG analyses from the ceRNA network and CNC network found "Phosphatidylinositol signaling system" was significantly enriched. Phosphatidylinositol is involved in many signal transduction pathways, such as PI3K/Akt pathway. PI3K/Akt signaling plays a central role in cellular physiology by mediating critical cellular processes, like glucose homeostasis. In db/db mice and high fat diet-induced diabetic mice, activates PI3K/Akt in an insulin-independent manner could attenuates hepatic gluconeogenesis and lipogenesis [51]. Zhang et al. found that artesunate and metformin combination improve salivary gland hypofunction in murine Sjögren's syndrome though regulating the PI3K/Akt pathway [52]. We speculated that these 5 validated DE lncRNAs might be key molecules regulating PI3K/