Transcriptome analysis reveals the mechanism of stromal cell-derived factor-1 and exendin-4 synergistically promoted periodontal ligament stem cells osteogenic differentiation

Stromal cell-derived factor-1 (SDF-1) and Exendin-4 (EX-4) play beneficial roles in promoting periodontal ligament stem cells (PDLSCs) osteogenic differentiation, while the detailed mechanism has not been clarified. In this study, we aimed to evaluate the biological mechanism of SDF-1 and EX-4 alone or synergistic application in regulating PDLSCs differentiation by RNA-sequencing (RNA-seq). A total of 110, 116 and 109 differentially expressed genes (DEGs) were generated in osteogenic medium induced PDLSCs treated by SDF-1, EX-4, and SDF-1+EX-4, respectively. The DEGs in SDF-1 group were enriched in signal transduction related signaling pathways; the DEGs in EX-4 group were enriched in metabolism and biosynthesis-related pathways; and the DEGs generated in SDF-1+EX-4 group were mainly enriched in RNA polymerase II transcription, cell differentiation, chromatin organization, protein phosphorylation pathways. Based on Venn analysis, a total of 37 specific DEGs were identified in SDF-1+EX-4 group, which were mainly enriched in negative regulation of autophagy and cellular component disassembly signaling pathways. Short time-series expression miner (STEM) analysis grouped all expressed genes of PDLSCs into 49 clusters according to the dynamic expression patterns and 25 genes, including NRSN2, CHD9, TUBA1A, distributed in 10 gene clusters in SDF-1+EX-4 treated PDLSCs were significantly up-regulated compared with the SDF-1 and EX-4 alone groups. The gene set enrichment analysis indicated that SDF-1 could amplify the role of EX-4 in regulating varied signaling pathways, such as type II diabetes mellitus and insulin signaling pathways; while EX-4 could aggravate the effect of SDF-1 on PDLSCs biological roles via regulating primary immunodeficiency, tight junction signaling pathways. In summary, our study confirmed that SDF-1 and EX-4 combined application could enhance PDLSCs biological activity and promote PDLSCs osteogenic differentiation by regulating the metabolism, biosynthesis and immune-related signaling pathways.


INTRODUCTION
Periodontal disease-induced tooth loss has become a global public health challenge that greatly affects people's quality of life (Peres et al., 2019). Mesenchymal stem cell (MSC) based periodontal tissue regeneration has aroused great attention in the field of regenerative medicine (Hu, Liu & Wang, 2018). Among all MSCs, periodontal ligament stem cells (PDLSCs) are the main candidate cells for periodontal regeneration. Several studies have demonstrated that transplanting autologous and allogeneic PDLSCs directly into periodontal defect areas or surgically created periodontal defect areas could result in periodontal tissue regeneration, which highlights that PDLSC-mediated tissue engineering is a promising therapy for periodontitis (Bartold, Shi & Gronthos, 2006;Ding et al., 2010;Liu et al., 2008;Liu et al., 2019a;Liu et al., 2019b).
An increasing number of researchers have focused on the recruitment of endogenous PDLSCs to the injury site to enhance healing by harnessing the innate regenerative potential of the body (Lee et al., 2017a). Cytokines, chemokines, and adhesion molecules have been used to enhance cell migration, maintain tissue homeostasis, regulate immune responses, promote wound healing and facilitate periodontal tissue regeneration (Lee et al., 2017b;Onizuka & Iwata, 2019;Wang et al., 2013). Stromal cell-derived factor-1 (SDF-1), a member of the chemokine family, can promote the proliferation and migration of various MSCs by activating the G protein-coupled receptor C-X-C chemokine receptor type 4 (CXCR4) (Du, Feng & Ge, 2016;Kimura et al., 2014;Zhu, Dissanayaka & Zhang, 2019). Our previous study also demonstrated that topical application of SDF-1 could significantly recruit MSCs to the wound area and promote local vascular regeneration in a rat model (Wang, Du & Ge, 2016). SDF-1 possesses great potential in promoting MSC migration and growth; however, the compromised osteogenic differentiation of these cells could not be induced by SDF-1. Therefore, the application of SDF-1 alone is insufficient for favorable bone regeneration, and the optimal method for potentiating periodontal bone regeneration is to combine SDF-1 with other osteogenic factors.
Exendin-4 (EX-4), a full agonist of glucagon-like peptide-1 receptor (GLP-1R), has been widely used in the clinical treatment of type 2 diabetes mellitus (T2DM) (Yap & Misuan, 2019). In addition, EX-4 plays key roles in promoting MSC proliferation and migration Zhou et al., 2015a). Recently, EX-4 has been confirmed to present the potential to promote osteogenic differentiation and bone formation in a variety of stem/precursor cells (Feng et al., 2016;Luciani et al., 2018;Meng et al., 2016). Moreover, in addition to enhancing the MSC osteogenic differentiation capability, EX-4 could promote the recruitment effect of SDF-1 (Zhou et al., 2015b). Our previous study also confirmed that SDF-1 and EX-4 cotherapy synergistically promoted PDLSCs proliferation, migration and osteogenic differentiation (Liang et al., 2021). However, the mechanism of SDF-1 and EX-4 alone or synergetic application for PDLSCs osteogenic differentiation is not fully understood.
High-throughput RNA sequencing (RNA-Seq) has been widely applied to analyze the whole transcriptomic changes of eukaryotes, which can provide progressively greater knowledge of both the quantitative and qualitative aspects of transcript biology (Ozsolak & Milos, 2011). RNA-seq has been successfully applied to identify the potential transcriptional mechanisms of various diseases, such as cancers, metabolic diseases and retinal diseases (Bakhtiarizadeh et al., 2019;Demircioğlu et al., 2019;Farkas et al., 2015). In the present study, RNA-seq transcriptomic analysis was applied to identify the core dynamic differentially expressed genes (DEGs) signature affected by EX-4 and SDF-1 alone or synergistically application in osteogenic medium-induced PDLSCs. Additionally, an integrated network containing specific DEGs generated in EX-4+SDF-1-treated PDLSCs was constructed. The results revealed the whole alteration of gene expression in PDLSCs undergoing EX-4 and SDF-1 application during the osteogenic differentiation process, which establishes a foundation for further research investigating the synergistic application of SDF-1 and EX-4 to promote PDLSCs osteogenic differentiation.

Human subjects and ethics statements
This study was approved by the Medical Ethical Committee of the Stomatology School, Shandong University (NO. 20170801). Five healthy individuals without systemic diseases aged from 18-30 who underwent premolar extraction at the Department of Oral and Maxillofacial Surgery were recruited for this project. All individuals agreed to participate in the research project and signed the informed consent forms according to the Helsinki Declaration of 1975.

Cell isolation and culture
The extracted teeth were stored in Dulbecco's modified Eagle's medium (DMEM, HyClone, Logan, UT, USA) with 5% antibiotics (100 U/mL penicillin, 100 mg/mL streptomycin, Sigma Aldrich, St Louis, MO, USA) and quickly transported from the clinic to the laboratory. Then, single PDLSCs were acquired as previously described in our previous study (Du, Feng & Ge, 2016). Specifically, Primary PDLSCs were cultured with DMEM containing 20% fetal bovine serum (FBS, BioInd, Kibbutz, Israel) at 37 • C in a humidified atmosphere of 5% CO 2 , and cells were trypsinized and passaged at a dilution ratio of 1:3 to expand the culture in 10% FBS medium upon the cell monolayer reached 80-90% confluence. Fourth passage cells were used in all experiments.

RNA-seq analysis
PDLSCs from 5 different individuals were cultured in osteogenic medium (OM) and treated with SDF-1, EX-4 or SDF-1+EX-4 at 21 d, and normal PDLSCs treated with OM served as a negative control (NC). Totally, we collected 20 samples and RNA-seq was used to analyze the whole genome expression at LC Sciences through the Illumina X10 platform (Hangzhou, Zhejiang, China). Firstly, the total RNA were performed quality control based on previous study (Wang, Wang & Li, 2012), and then the clean reads were mapped to the reference genome (GRCh38) via hierarchical indexing for spliced alignment of transcripts (HISAT) (v2.0.4) (Kim, Langmead & Salzberg, 2015). The mapped reads of each sample were assembled using StringTie (v1.3.0, Pertea et al., 2015). Furthermore, all transcriptomes samples were merged to reconstruct a comprehensive transcriptome using Perl scripts. After the final transcriptome was generated, StringTie and edgeR were used to estimate the expression levels of all transcripts. StringTie was used to determine the expression level of mRNAs by calculating fragments per kilobase of exon model per million mapped fragments (FPKM). The DEGs were selected with statistical significance (p value <0.05) by R package (v 3.2.5).

Volcano analysis of DEGs in PDLSCs
Volcano analysis was used to identify the DEGs between each pair of groups (Li, 2012). The up-and down-regulated genes were identified, and the total number of each pair of groups was visualized by the histogram.

Gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis
Based on the DEGs generated by SDF-1, EX-4 and SDF-1+EX-4 compared with NC, the overrepresented GO categories and the significant KEGG pathways were identified (anonymous, 2019;Kanehisa & Goto, 2000). A q value <0.05 was used as the cut-off criterion for the selection of significant GO terms and KEGG pathways.

Venn and Upset analysis
To identify the specific DEGs generated by every two different compared groups, overlapping analysis was performed according to the jvenn website (http://jvenn.toulouse. inra.fr/app/example.html), and an intersection UpSet diagram based on the UpSet R package was drawn (Conway, Lex & Gehlenborg, 2017). The specific genes generated by each group were identified and the gene functions were analyzed according to Metascape website (http://metascape.org/gp/index.html#/main/step1) (Zhou et al., 2019).

Short time-series expression miner (STEM) analysis
STEM software (version 1.3.8) (Douglas et al., 2019) was applied to identify the specific gene expression clusters in PDLSCs treated with SDF-1, EX-4 and SDF-1+EX-4. The genes in the upregulated clusters in SDF-1+EX4-treated PDLSCs were selected, and the expression of all genes in these clusters was shown by a heatmap. The interaction relationship of these genes was analyzed according to the GeneMANIA database (http://genemania.org/search) (Franz et al., 2018).

In-depth mechanism analysis and functional network construction
To identify the function of specific DEGs generated by SDF-1, EX-4 and SDF-1+EX-4treated PDLSCs, a functional network was constructed according to STRING database (https://string-db.org/) and GeneMANIA database (http://genemania.org/search) (Franz et al., 2018). Functional enrichment analysis of genes in the functional network was further performed by the Metascape database (Zhou et al., 2019).

Gene set enrichment analysis (GSEA)
GSEA is one of the functional class scoring analysis methods (Subramanian et al., 2005). To select the genes that were not significantly differentially expressed but were important for the function of biological pathways, GSEA was performed method according to all genes in SDF-1 vs NC, EX-4 vs NC, SDF-1+EX-4 vs NC, SDF-1+EX-4 vs EX-4, SDF-1+EX-4 vs SDF-1, and EX-4 vs SDF-1 groups according to the clusterProfiler and enrichplot R package.

GO and KEGG enrichment analysis
The GO (biological process) enrichment results showed that the DEGs generated by SDF-1 stimulation were mainly enriched in signal transduction regulation of transcription, cell differentiation, RNA polymerase II transcription, and multicellular organism development signaling pathways ( Fig. 2A and Table S1). The KEGG analysis results also indicated that SDF-1-induced DEGs were mainly enriched in the glucagon signaling pathway, NF-kappa B signaling pathway, and TNF signaling pathway ( Fig. 2B and Table S1), and mostly involved in human diseases and organismal systems, such as viral infectious diseases, the immune system, and the endocrine system ( Fig. S2B and Table S1). In addition, the DEGs generated by the EX-4 stimulation were mainly enriched in the oxidationreduction process, metabolic process, regulation of transcription, DNA-templated, protein phosphorylation, and lipid metabolic process pathways ( Fig. 2C and Table S1). KEGG analysis also indicated that EX-4-induced DEGs were mainly enriched in metabolism and biosynthesis-related pathways, such as nicotinate and nicotinamide metabolism, phenylpropanoid biosynthesis, and phenylpropanoid biosynthesis ( Fig. 2D and Table S1), and all these DEGs were significantly involved in energy metabolism, lipid metabolism and biosynthesis of other secondary metabolites (Fig. S2D). Moreover, we identified that the DEGs generated by SDF-1+EX-4 was also mainly enriched in RNA polymerase II transcription, cell differentiation, chromatin organization, and protein phosphorylation pathways according to the GO (biological process) analysis results ( Fig. 2E and Table S1). The KEGG enrichment analysis results indicated that these costimulated DEGs were mainly enriched in pathogenic Escherichia coli infection, mitophagy in animals, bacterial invasion , and EX-4 vs SDF-1 (F). Each point represents individual genes. Black dots represent the genes that were not significantly differentially expressed, red dots indicate the genes that were significantly upregulated and blue dots indicate the genes that were downregulated (|FC| > 1 and p-adjusted value < 0.05). The volcano plots were analyzed by DESeq2.
Full-size DOI: 10.7717/peerj.12091/ fig-1 of epithelial cells and cysteine and methionine metabolism signaling pathways ( Fig. 2F and Table S1). Additionally, we found that all the DEGs generated by the SDF-1, EX-4, SDF-1+EX-4 groups were significantly enriched in the nucleus, membrane, nucleosome and cytoplasm pathways according to the GO cellular component analysis, and based on the GO molecular function analysis, all these DEGs were involved in the DNA binding, protein binding, ATP binding and protein heterodimerization activity pathways ( Figs

Screening of specific DEGs and functional analysis
To characterize the DEGs specifically generated by SDF-1 vs NC, EX-4 vs NC, SDF-1+EX-4 vs NC, SDF-1+EX-4 vs EX-4, SDF-1+EX-4 vs SDF-1, and EX-4 vs SDF-1, we performed an overlapped analysis of the six compared DEG groups, and the results indicated that 34, 30, 37, 39, 35 and 28 DEGs were specially generated (Fig. 3A, Table S2, Fig. S3 and Table S3). In addition, the Metascape website was referenced to analyze the specific gene functions in different groups and the results indicated that the 34 specific DEGs generated in SDF-1 vs NC was mainly enriched in response to the wounding, and macromolecule methylation signaling pathways; the 30 specific DEGs generated in EX-4 vs NC were mainly enriched in the cellular response to external stimulus, regulation of protein complex assembly, and regulation of chromosome organization signaling pathways; the 37 specific DEGs  Table S4).

Screening of core DEGs generated by the SDF-1 and EX-4 combined application
To identify the core DEGs generated by the SDF-1 and EX-4 combined application, STEM software was applied to perform a pattern analysis (mock infection was designated NC), and the results revealed 49 gene clusters among all DEGs generated by SDF-1, EX-4 and SDF-1+EX-4-treated PDLSCs ( Fig. S4 and Table S5). In all gene clusters, we focused on genes in clusters 14,39,47,36,7,28,45,25,18, and 41, which were significantly upregulated by the SDF-1 and EX-4 combined application compared with SDF-1 or EX-4 alone (Fig. 4A). Among the 10 clusters, a total of 25 genes, including NRSN2, CHD9, TUBA1A, AKAP13, VAMP7, NPIPA2 etc., were identified, and all DEG expression in different groups is shown by a Heatmap (Fig. 4B and Table S6). The GeneMANIA analysis network showed that these genes possessed an intricate tangle of connections through genetic interactions and co-expression pattern ( Fig. 4C and Table S7). Functional enrichment analysis indicated that our core DEGs and their related genes constructed in the network were mainly enriched in the metallothioneins binding metals, asparagine N-liked glycosylation, response to stimulus, detoxification, biological regulation and growth signaling pathways ( Fig. 4D and Table S8).

Network construction and the mechanism analysis of the SDF-1 and EX-4 synergistic effect in PDLSCs
To further clarify the mechanism of the SDF-1 and EX-4 combined application to promote PDLSC osteogenic differentiation, a detailed complex network analysis was performed based on the DEGs generated by the different groups through the STRING and GeneMANIA. A network analysis based on the STRING database showed that the DEGs generated by EX-4 and SDF-1 alone were centralized and converged into a network, while the DEGs in the SDF-1+EX-4 group were scattered in the network. The DEGs generated in SDF-1+EX-4 vs SDF-1, and SDF-1+EX-4 vs EX-4 were also centralized (Fig. S5). In addition, compared with the NC group, EX-4 significantly elevated the expression of the core gene MAPK27, which plays central roles by activating MAPK-related signaling pathways; and SDF-1 mainly changed the expression of the key genes CREB1, MMP13, RHOQ and BIRC2, which exert their roles by activating ATP signaling pathways (Figs. 5A and 5B). More importantly, we identified that the gene expression levels of NEDD8, CHCHD1, LMO7, and ATP5L were significantly varied the SDF-1+EX-4 vs EX-4 groups, which played crucial roles in activating ATP-related signaling pathways and indicated that SDF-1 could amplify the EX-4 effect in elevating ATPase activity and promoting PDLSC osteogenic differentiation. Additionally, we found that the DEGs of SOX15, UBC, VAMP7 and ARPC5 were significantly changed in SDF-1+EX-4 vs SDF-1 groups, which played vital roles in promoting cytoskeletal protein formation and degradation (Figs. 5C and 5D). Although the DEGs generated in the SDF-1+EX-4 vs NC group were dispersed, most of these genes could be centralized by the genes HAP1, TP53, TAL1, RRPF40A and TALDO1 (Fig. S6).

Gene set enrichment analysis (GSEA) based on all genes
All of the above analyses were based on the DEGs selected according to our threshold value and statistical analysis technique; however, genes that are not significantly differentially expressed but are important for the function of the SDF-1 and EX-4 biological pathways may be omitted, and thus, we performed a GSEA according to all genes in SDF-1 vs NC, EX-4 vs NC, SDF-1+EX-4 vs NC, SDF-1+EX-4 vs EX-4, SDF-1+EX-4 vs SDF-1, and EX-4 vs SDF-1. The results indicated that SDF-1 played its roles in regulating the metabolism of xenobiotics by cytochrome P450, regulation of autophagy, neuroactive ligand receptor interaction signaling pathways ( Fig. 6A and Table S9); EX-4 played roles mainly through affecting the metabolism related signaling pathways, such as, starch and sucrose metabolism, arginine and proline metabolism, and type II diabetes mellitus signaling pathways (Fig. 6B and Table S9); and the combination of SDF-1 and EX-4 mainly regulating PLDSCs biological process through activating metabolism and immunity related signaling pathways, such as pantothenate and COA biosynthesis, vascular smooth muscle contraction, and maturity onset diabetes of the young (Fig. 6C and Table S9). In addition, we confirmed that SDF-1 could amplify the role of EX-4 in regulating various signaling pathways, such as type II diabetes mellitus, the insulin signaling pathway, and the allograft rejection pathway (Fig. 6D and Table S9), while EX-4 could aggravate the effect of SDF-1 on PDLSC biological roles by regulating signaling pathways, including primary immunodeficiency, tight junction, and basal transcription factors ( Fig. 6E and Table S9). Interestingly, we found that the signaling enrichment analysis based on the DEGs and GSEA shared similarity, which indicated that SDF-1 and EX-4 may regulate PDLSC biological activity by enhancing each other's biological functions.

DISCUSSION
The effects of SDF-1 and EX-4 alone on bone regeneration have been widely reported; however, the synergetic effects of SDF-1 and EX-4 on PDLSC osteogenic differentiation and the potential mechanism have not been reported. In our previous study, we confirmed that the combined application of SDF-1 and EX-4 could significantly promote osteogenic differentiation of PDLSCs (Liang et al., 2021). In this study, we confirmed that SDF-1 could amplify the role of EX-4 by significantly activating metabolism-related signaling pathways, such as type II diabetes mellitus and insulin signaling pathways; and EX-4 could aggravate the effect of SDF-1 on PDLSCs biological roles by regulating primary immunodeficiency and tight junction signaling pathways. Briefly, our study confirmed that the SDF-1 and EX-4 combined application could enhance PDLSCs biological activity and promote PDLSCs osteogenic differentiation by regulating the metabolism, biosynthesis and immune-related signaling pathways. EX-4, a common glucagon-like peptide-1 receptor agonist, has been confirmed to possess excellent effects on treating patients with T2DM by significantly reducing HbA1c content compared to basal insulins (Singh et al., 2017). Our current study confirmed that EX-4-generated DEGs were enriched in the type II diabetes mellitus signaling, starch and sucrose metabolism, arginine and proline metabolism, and alanine aspartate and glutamate metabolism signaling pathways, while the combination of SDF-1 and EX-4 could significantly activate more metabolism-related signaling pathways, such as the valine leucine and isoleucine degradation, insulin signaling, phenylalanine metabolism, and pyruvate metabolism signaling pathways. In the plasma of type II diabetes, a pronounced postprandial rise in amino acids (such as leucine, isoleucine, valine, lysine, and threonine) and glucose-dependent insulinotrophic polypeptide was observed, which finally resulted in glycemic and insulinemic responses (Nilsson, Holst & Björck, 2007). In our study, we confirmed that SDF-1 and EX-4 combination therapy could significantly increase the gene expression of EHHADH, HMGCL, IL4I1, PKLR, PIK3CG, PYGM, SLC2A4, RHOQ, PRKCZ, FBP1, and SH2B2, which play key roles in promoting amino acid degradation and insulin secretion; thus, SDF-1 assists EX-4 in controlling the blood glucose level of diabetes patients. EX-4 could promote the osteogenic differentiation of osteoblasts, adipose-derived stem cells, and PDLSCs by activating the Hedgehog/Gli1, Wnt and NF-κB signaling pathways; however, the interactive relationship of EX-4 regulating PDLSCs metabolism-related pathways and osteogenic differentiation pathways and the mechanism by which SDF-1 amplifies the effect of EX-4 on PDLSCs osteogenic differentiation need to be further clarified (Deng et al., 2019;Gao, Li & Li, 2018;Liu et al., 2019a;Liu et al., 2019b).
SDF-1 could promote PDLSCs proliferation, migration and osteogenic differentiation in vitro, play key roles in recruiting endogenic PDLSCs into the periodontal defect area and then contribute to angiogenesis in vivo (Bae et al., 2017;Du, Feng & Ge, 2016). In the current study, we confirmed that SDF-1 significantly activated metabolism-related signaling pathway of PDLSCs cultured in OM, such as the metabolism of xenobiotics by cytochrome P450, alanine aspartate and glutamate metabolism; however, the effects of metabolism on PDLSCs osteogenic differentiation have not been reported. In addition, SDF-1 activated the renin angiotensin system, basal transcription factor, and neuroactive ligand receptor interaction signaling pathways, which is consistent with previous studies (Chu et al., 2009;Wang et al., 2015). Moreover, we found that the EX-4 and SDF-1 combined stimulation significantly activated the immunodeficiency, tight junction, complement and coagulation cascade signaling pathways compared to SDF-1 stimulation alone in osteogenic medium cultured PDLSCs. The adaptive immune system plays a prominent role in the development of heterotopic ossification (Ranganathan et al., 2016), and tight junctions (TJs) play a pivotal role in the modulation of paracellular permeability. For example, Cldn11 recombinant protein, a well-established component of TJs, could significantly decrease the resorption of lipopolysaccharide-induced calvarial bone and increase the osteogenic activity of calvarial bone formation (Baek et al., 2018). Moreover, stress-induced hematopoietic stem cell mobilization is enhanced by the fibrinolytic and complement cascades (Nguyen, Lapidot & Ruf, 2018). In the current study, we confirmed that EX-4 could enhance the role of SDF-1 in PDLSCs osteogenic differentiation through RNA-seq analysis, although this process was limited due to the lack of biological experimental validation. In the future, we will design related biological experiments to validate the mechanism of the SDF-1 and EX-4 combined application in promoting PDLSCs osteogenic differentiation.

CONCLUSION
Our study confirmed that SDF-1 could amplify the role of EX-4 by activating more metabolism-related signaling pathways, such as type II diabetes mellitus and the insulin signaling pathways, and EX-4 could aggravate the effect of SDF-1 on PDLSCs biological roles by regulating primary immunodeficiency and tight junction signaling pathways. In addition, we confirmed that the SDF-1 and EX-4 combined application could enhance PDLSCs biological activity and promote PDLSCs osteogenic differentiation by regulating metabolism, biosynthesis and immune-related signaling pathways. Our current study lays a solid foundation for exploring the effects of SDF-1 and EX-4 synergistic application on periodontal tissue regeneration.

Grant Disclosures
The following grant information was disclosed by the authors: