Effect of Tensile Frequency on the Osteogenic Differentiation of Periodontal Ligament Stem Cells

Introduction

Mechanical stress enhances bone metabolism and periodontal tissue remodeling.¹ During orthodontic tooth movement (OTM), bone remodeling is initiated via the periodontal ligament.² As the main mesenchymal stem cells (MSCs) in the periodontal ligament, periodontal ligament stem cells (PDLSCs) play an important role in mechanical signal transduction. Currently, a consensus has been reached that cyclic mechanical tension is a strong driver of the differentiation of PDLSCs into the osteoblast lineage.^2–4 Mechanical tension activates calcium channels,⁵ which activate the ERK1/2 and P38 MAPK pathways through integrin-FAK or protein kinase (PKC)-SR signaling⁶ and induce the phosphorylation of Runt-associated transcription factor 2 (Runx2),⁷ promoting osteogenic precursor cell synthesis and the transcription of mineralizable proteins.⁸ At present, the TGF-β, BMP, MAPK, Notch, Wnt, Hedgehog, FGF, and Hippo signaling pathways have been found to be involved in this process.

Force parameters (including magnitude, frequency, and duration) are crucial for well-regulated tissue remodeling. However, numerous in vitro studies performed to date show enormous heterogeneity in tensile force parameters.⁹ In different studies, cyclic tension was applied with magnitudes ranging from 1% to 24%, frequencies ranging from 0.1 Hz to 1.0 Hz, and stimuli duration ranging from 1 hour to 6 days, thereby reducing comparability between different studies.⁹ To establish strategies to optimize tensile force parameters, it is of particular importance to understand how different tensile force parameters affect the osteogenic differentiation of PDLSCs.

The effects of different tensile force magnitudes and durations have been investigated in some studies. Among the magnitudes, a magnitude of 10% generally led to a lower level of inflammation and a higher level of osteogenesis,¹⁰ whereas a magnitude of 12% was found to correlate well with strain conditions at the mid-root under physiological loading conditions^11,12 and to induce optimal effects in both the proliferation and osteogenesis of PDLSCs.¹³ Cyclic tension alone at 3000 μstrain significantly enhanced SATB Homeobox 2 (Satb2) after 3 h of loading and significantly upregulated Runx2 after 6 h.¹⁴ The synthesis of BMP9 increased under 6-h continuously applied cyclic tension.¹⁵ In addition, 12% cyclic tensile force gradually upregulated the expression of Runx2, alkaline phosphatase (ALP), and osteocalcin (OCN) with force durations of 6 h, 12 h, and 24 h, respectively.^16,17 The protein level of osterix increased stepwise following 3 h, 6 h, 12 h, and 24 h of exposure to tensile strain.¹⁴ Recently, temporal gene expression patterns were delineated.¹⁷

Tensile frequency varies largely among different studies. The ROCK‐TAZ pathway and its interaction with Cbfα1 were found to be essential for the cyclic tension (12% elongation, 0.1 Hz)-induced osteogenic differentiation in PDLSCs.¹⁸ Cyclic tension (10% elongation, 0.5 Hz) stimulated the osteogenic differentiation of PDLSCs by inhibiting miR-129-5p expression and activating the BMP2/Smad pathway.¹⁷ LncRNAs-miRNAs-mRNAs networks in PDLSCs were depicted under cyclic tension (10% elongation, 1.0 Hz).¹⁹ However, there have been rare studies examining the impact of different cyclic tensile frequencies on osteogenesis of PDLSCs and the expression of relevant genes thus far. The low-magnitude high-frequency (LMHF) vibration approach was excluded because it is used to simulate a masticatory force, while cyclic tension is used to simulate an orthodontic force, and the two methods of force application are completely different.^2,20 Previous animal studies on long bone distraction osteogenesis have shown that loading frequency affects the osteogenic response of bone tissue.²¹ The mechano-regulation of trabecular bone adaptation is logarithmically dependent on the loading frequency.²² Therefore, we hypothesized that tensile frequency would affect osteogenesis of PDLSCs, in which some tensile frequency-sensitive genes may play an important role. To test our hypotheses, human PDLSCs were subjected to cyclic mechanical tension at different frequencies of 0.1–0.7 Hz to examine the osteoblastic differentiation of PDLSCs, and high-throughput sequencing was performed to characterize the frequency-course expression patterns of mRNA during the osteogenic differentiation of PDLSCs. This study aimed to investigate the effects of tensile frequency on the osteogenic differentiation of PDLSCs as well as the relevant molecular mechanisms.

Materials and Methods

Cultivation of Human PDLSCs

Healthy periodontal ligament tissues were scraped from the middle third of tooth roots, which were extracted for orthodontic reasons, with informed consents. All donors were aged from 14 to 16 years and had no systemic or oral diseases. The periodontal ligament tissues were cut into small pieces and enzymatically digested for 40 min at 37°C with collagenase I (3 mg/mL, Sigma-Aldrich, St. Louis, MO, USA) and dispase II (4 mg/mL; Sigma-Aldrich). The cells were seeded in 25 cm² flasks (Falcon, BD Biosciences, Franklin Lakes, NJ, USA) with α-minimal essential medium (α-MEM; Sigma-Aldrich) supplemented with 10% fetal bovine serum (FBS, Gibco, Life Technologies Co., Grand Island, NY, USA) and antibiotics (100 U/mL penicillin and 100 µg/mL streptomycin, Hyclone, Logan, Utah, USA), and incubated in a humidified atmosphere (37°C, 5% CO₂). The medium was changed every 3 days. After reaching 80% of confluence, the cells were detached with 0.25% trypsin/EDTA (Gibco, Life Technologies Co., Grand Island, NY, USA), and single-cell suspensions were cloned with the limiting-dilution method to purify the stem cells.¹⁷ Cell clusters from the colony were trypsinized and serially sub-cultured.

Immunofluorescent Assay

The third-passage PDLSCs were sub-cultured into six-well plates until confluent. The culture medium was then removed, and cells were fixed with 4% formaldehyde (Zhonghuihecai, Xi’an, CN) for 20 min, permeabilized with 0.3% Triton X-100 (Zhonghuihecai) for 5 min, and incubated with primary antibodies (anti-pan-cytokeratin, 1:300, Abcam, Cambridge, MA, USA; anti-vimentin, 1:500, Abcam; anti-STRO-1, 1:200, Abcam; anti-CD146, 1:200, Abcam) overnight at 4°C. The cells were then washed with PBS and incubated with CY3/FITC-conjugated secondary antibodies (1:500, Zhuangzhi, Xi’an, CN) in darkness for 30 min, and then washed with PBS. Finally, the nuclei were counterstained with 4, 6-diamidino-2-phenylindole (DAPI, Zhuangzhi), and fluorescent images were captured with a fluorescence microscope (Olympus, Japan).

Multilineage Differentiation of PDLSCs

For osteogenic and adipogenic differentiation, PDLSCs were seeded into six-well plates at a density of 2×10⁵ cells/well, and after reaching 80% confluence, the medium was replaced with an osteogenic or adipogenic inductive medium (Osteogenesis or Adipogenesis Differentiation Kit, Cyagen, USA). Seven days after osteogenic induction, the cells were stained with ALP (ALP staining kit, Solarbio, CN), and 14 days after osteogenic incubation, the cells were stained with Alizarin Red (ARS staining kit, Cyagen, USA). After 21 days of adipogenic incubation, the cells were stained with Oil Red O (Cyagen, USA).

Application of Cyclic Tension

Flexcell FX-5000T Tension Plus System (Flexcell International Corporation, Hillsborough, NC, USA) was used to mimic the tensile force exerted on PDLSCs during OTM, according to previous studies.²³ Cyclic tensile loading experiments were performed on the fourth-passage PDLSCs from four different healthy donors in triplicate. PDLSCs were seeded onto six-well type I collagen (COL-I)-coated silicone culture plates (Flexcell International Corporation) at a density of 2×10⁵ cells/well. Upon reaching 80% confluence, the cells were serum-starved overnight, and the medium was changed to osteogenic medium (Cyagen, USA). Cyclic tensile force (12% bottom membrane elongation) was applied to different plates at different frequencies of 0.1 Hz, 0.5 Hz, and 0.7 Hz. Control cells were cultured under identical culture condition but without mechanical stimulation.

Western Blot Analysis

After 10 h of cyclic tensile force, the cells were collected. Proteins were isolated, electrophoretically separated, and immunoblotted as previously described.²³ Briefly, PDLSCs were lysed with RIPA buffer containing 1% phenylmethanesulfonyl fluoride (PMSF, proteinase inhibitor, Zhonghuihecai) and 1% phosphatase inhibitor (Zhonghuihecai). After centrifugation, the supernatant was collected and measured quantitatively using a BCA Protein Assay Kit (Absin, Shanghai, CN). Total protein from cell lysates (20 µg/lane) was separated by SDS-PAGE gels (Beyotime, Hangzhou, CN) and then transferred onto a polyvinylidene difluoride (PVDF) membrane (EMD Millipore, Billerica, MA, USA). After blocking with 5% skimmed milk in tris-buffered saline tween-20 (TBST) for 2 h at room temperature, the membranes were incubated overnight at 4°C with primary antibodies (runt-related transcription factor 2 (Runx2), 1:500, ImmunoWay, USA; COL-I, 1:1000, Proteintech, USA; Glyceraldehyde-3-phosphate dehydrogenase (GAPDH), 1:10000, Proteintech, USA), followed by incubation with horseradish peroxidase (HRP)-conjugated secondary antibody (1:2000, Proteintech, USA) for 2 h at room temperature. The protein expression was visualized using ChemiDocTM XRS+ (Bio-Rad Laboratories, Inc., Hercules, CA, USA) with an enhanced chemiluminescence (ECL) kit (Millipore, Billerica, MA, USA). GAPDH was used as an internal control for normalization.

RNA Sequencing

The fourth-passage PDLSCs from three donors were used for RNA sequencing after tension loading. Total RNA was extracted from the four groups of cells (normal PDLSCs and PDLSCs tensioned at frequencies of 0.1 Hz, 0.5 Hz, and 0.7 Hz for 6 h) using the Trizol (Sigma-Aldrich), according to the manufacturer’s protocols. After digestion with DNase, rRNA were depleted using a Ribo-Zero magnetic kit, and sequencing libraries were constructed as previously described.²⁴ The sequencing of the cDNA library was carried out by Gene Denovo Biotechnology Co. (Guangzhou, China). The gene expression level was evaluated by reads per kilobase transcriptome per million mapped reads (RPKM). Requirements for filtering differentially expressed genes (DEGs) were as follows: (1) |log2 (fold-change)| ≥ 1; (2) p value < 0.05. DESeq2 (differential gene expression analysis based on the negative binomial distribution)²⁵ was used to calculate p values and adjusted p (adj. p) values. Heatmaps and volcano plots analyses were used to visualize these DEGs using the Complex Heatmap package and ggplots2 package of “R” software. The online tool Venny 2.1 (http://bioinfogp.cnb.csic.es/tools/venny/index.html) was applied to identify the common DEGs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were used for annotation visualization and integrated discovery. Raw data of the performed RNA sequencing (RNA-seq) were recorded in the SRA database with the SRA accession: PRJNA665587.

Frequency Series Clustering Analysis and Hub Gene Identification

Gene expression pattern analysis is used to cluster genes of similar expression patterns for multiple samples in a tensile frequency order. To examine the expression pattern of DEGs in different frequencies, the expression data of each sample (in the order of treatment) were normalized to 0, log2 (v1/v0), and log2 (v2/v0), and then clustered using the Short Time-series Expression Miner software (STEM).²⁶ The parameters were set as follows: (1) maximum unit change in model profiles between frequency points was 1; (2) maximum output profiles number was 20 (similar profiles will be merged; and (3) minimum ratio of fold change of DEGs was no less than 2.0. The clustered profiles with p value < 0.05 were considered as significant profiles. Then, the DEGs in all profiles or in each profile were subjected to GO and KEGG pathway enrichment analysis. STRING (https://string-db.org/) was used for PPI network analysis. The DEGs were mapped onto the PPI network with a minimum interaction score of 0.4. Cytoscape v3.7.1 software was used to visualize the PPI network. Gene network clustering analysis was performed to identify the key PPI network modules, using the MCODE and cytoHuba app from the Cytoscape software suite. Adj. p value < 0.05 was set as the significance threshold.

Quantitative Real-Time PCR

Total RNA was extracted using Trizol (Sigma-Aldrich) according to the manufacturer’s protocols. Quantitative real-time PCR (RTqPCR) was performed in triplicate using the Power SYBR Green PCR Mastermix (Applied Biosystems, Foster City, CA, USA). Sequences of the primers used are shown in Table S1. The mean expression values were calculated relative to GAPDH, which was used as an internal control for normalization.

Statistical Analysis

For differential gene expression analysis, a likelihood ratio test was used assuming an underlying zero-inflated negative binomial distribution.^25,27 FDR-corrected p-values were calculated using the Benjamini-Hochberg procedure.²⁸ K-means cluster method²⁹ was used for frequency cluster expression pattern analysis by the STEM software. P and Q values of GO/KEGG functional analysis and frequency cluster analysis were calculated based on a hypergeometric distribution.³⁰ Values of the relative protein and mRNA expression were expressed as mean ± SD within each group. One-way ANOVA followed by SNK post hoc tests was used in the Western blotting assay and RTqPCR assay. The significance threshold was set at 0.05.

Results

Characterization of PDLSCs

Immunofluorescent staining revealed that the isolated cells were positive for CD146, vimentin, and STRO-1 (Figure 1A–C, respectively) but negative for pan-cytokeratin (Figure 1D), confirming that the cells were mesenchymal stem cells of mesodermal origin. The osteogenic potential of PDLSCs was determined by positive ALP staining after 7 days of osteogenic induction (Figure 1E and F) and red mineralized matrix nodules in Alizarin Red staining after 14 days of osteogenic induction (Figure 1G and H). The presence of red lipid droplets in Oil Red staining after 21 days of adipogenic induction (Figure 1I) indicated the adipogenic differentiation of PDLSCs.

Figure 1 Identification of PDLSCs. The Immunofluorescence showed that the cells were positive for CD146 (A, red), vimentin (B, green), and STRO-1 (C, red) but negative for pan-cytokeratin (D). Scale bar = 50 μm. ALP staining was positive both visually (E, black) and under the microscope (F, black, Scale bar = 100 μm), after 7 days of osteogenesis induction culture. After 14 days of osteogenesis induction, ARS staining was observed to be positive by the naked eye (G, red), and mineralized nodules were obvious under the microscope (H, red, Scale bar = 100 μm). After 21 days of adipogenic induction, oil red O-positive lipid clusters were observed microscopically (I, red, Scale bar = 100 μm).

Effect of Tensile Frequency on Osteogenic Differentiation of PDLSCs

PDLSCs were exposed to 12% cyclic tension, which contributed to the cellular reorientation, including an irregular arrangement at the central region and a parallel arrangement at the peripheral region of the plate (Figure 2A). As shown in Figure 2B and C, after consecutive 10-h cyclic tensile force loading, the protein expression levels of Runx2 and COL‑I increased with increasing tensile frequency from 0.1 Hz to 0.7 Hz and were remarkably higher than that in the group without tension application (p < 0.05). The result suggested that mechanical tension (12% deformation) upregulated the osteogenesis of PDLSCs in a frequency-dependent manner. Higher frequencies of cyclic tension were associated with higher osteogenic differentiation of PDLSCs.

Figure 2 Cyclic tension promoted the osteogenic differentiation of PDLSCs. Under cyclic equibiaxial tension, PDLSCs reoriented in parallel alignment at the peripheral region of the plate, while in random orientation at the center of the plate (A, Scale bar = 100 μm). Western blotting was used to detect protein levels of osteogenesis-related genes, COL-I and Runx2, at different tensile frequencies (B and C). **p < 0.01, ***p < 0.001, vs control group.

Characterization of the mRNA Transcriptome of PDLSCs

It has been previously shown that cyclic mechanical tensile stress can improve osteogenesis of PDLSCs, and that consecutive 4–6 h of tension can significantly upregulate the mRNA expression of osteogenesis-related genes.^16,17 To gain insight into the molecular mechanisms by which mechanical tension stimulates osteoblast differentiation of PDLSCs, total RNA was extracted from PDLSCs to conduct RNA-seq after 6 h of cyclic tension. The mRNA expression profiles of the PDLSCs at different tensile frequencies (0.1 Hz, 0.5 Hz, and 0.7 Hz) were detected. Comparative expression analyses were performed according to the different frequencies of the tensile stress (0.1 Hz vs control, 0.5 Hz vs control, and 0.7 Hz vs control). Heat maps of the top 40 DEGs (Figure 3A–C) and volcano plots (Figure 3D–F) were depicted. In total, 50 mRNAs were upregulated, and 261 mRNAs were downregulated at 0.1 Hz. At 0.5 Hz, 656 mRNAs were upregulated, and 1474 mRNAs were downregulated. At 0.7Hz, 139 mRNAs were upregulated, and 194 mRNAs were downregulated. A Venny analysis (Figure 3G and H) showed that 78 genes were simultaneously upregulated and 118 were simultaneously downregulated among the 0.1 Hz, 0.5 Hz, and 0.7 Hz groups. The GO analysis (Figure 4A) demonstrated that changes in biological processes (BPs) were mainly enriched in “metabolic process”, “response to stimulus”, “biological regulation”, “signaling”, and “localization”. Changes in Cellular Components (CCs) were mainly enriched in “organelle”, “membrane”, “macromolecular complex”, and “membrane-enclosed lumen”. Moreover, “binding”, “catalytic activity”, and “nucleic acid binding transcription factor activity” emerged as the highest-ranked Molecular Function (MF) groups. As shown in Figure 4B, DNA replication, cell cycle, and the TNF signaling pathway were significantly enriched in the KEGG pathway. Within the primary category “Environmental Information Processing”, “Signal transduction”, and “Signaling molecules and interaction” were strongly enriched (Figure 4C).

Figure 3 Identification of DEGs among different frequencies. (A–C) Heatmaps of the top 40 DEGs between 0.1 Hz/0.5 Hz/0.7 Hz and static culture, respectively. Red rectangles represent high expression, and blue rectangles represent low expression. (D–F) Volcano plot of DEGs between 0.1 Hz/0.5 Hz/0.7 Hz and static culture, respectively. The red plots represent upregulated genes, the blue plots represent downregulated genes, and the black plots represent nonsignificant genes. (G) Venn diagram of upregulated DEGs among 0.1 Hz, 0.5 Hz, and 0.7 Hz. (H) Venn diagram of downregulated DEGs among all the three frequencies.

Figure 4 Functional enrichment analysis of all DEGs. (A) GO enrichment analysis of all DEGs among different frequencies. (B) Top 20 pathways of the KEGG enrichment analysis of all DEGs among different frequencies, with the KEGG pathway annotation (C). The screening criteria for significance were p value < 0.05.

Gene Expression Patterns of mRNAs at Different Tensile Frequencies

The sequencing data were normalized to the control, and trend analyses of DEGs were identified using STEM. In Figure 5A, within the 20 model profiles, eight mRNA trend profiles were statistically significant. The profile number assigned by STEM was on the top left corner of each profile box, p value was on the bottom left, and the number of the cardinality of each cluster was on the top right corner. As shown in Table 1, among different profiles, the top-ranked KEGG pathways were mainly in the metabolic pathways, PI3K-Akt signaling pathway, cytokine-cytokine receptor interaction, and MAPK signaling pathway. A continuous downregulation pattern was found in profile 0 (Figure 5B), in which the high-ranked BPs, CCs, and MFs in GO enrichment (Figure 5C) were similar to those in Figure 4A, and inflammatory pathways such as arachidonic acid metabolism, peroxisome, and cytokine−cytokine receptor interaction were strongly enriched (Figure 5D).

Table 1 Top 10 of KEGG Enrichment Among Different Profiles

Figure 5 Frequency series clustering analysis on expression profiles of mRNAs by STEM. (A) Within the 20 model profiles, eight mRNA trend profiles were statistically significant. The number at the upper-left corner of each profile box was the profile number assigned by STEM, the number on the bottom left was the p value, and the number on the top-right corner was the number of genes within each cluster. (B) Persistently downregulated genes along frequency were clustered in profile 0. (C) GO enrichment of profile 0. (D) Top 10 pathways of the KEGG enrichment of profile 0.

The interactions of 194 DEGs in profile 0 were analyzed using the STRING online database, and the PPI network was obtained using the Cytoscape software (Figure 6A). The MCODE plugin was then used to investigate the key PPI network modules, and one key module with four genes (EYA1, SIX5, SALL1, FRAS1) was identified (Figure 6B). The cytoHubba plugin was then used to analyze hub genes with maximum correlation criterion (MCC)/Degree, and genes with the top 10 scores were respectively identified. The intersection (EYA1, SALL1) of hub genes according to the above three methods were selected for further RTqPCR validation (Figure 6C). The results of RTqPCR (Figure 6D and E) showed that the mRNA expression of EYA1 and SALL1 decreased with increasing frequency from 0.1 Hz to 0.7 Hz, which were highly consistent with our high-throughput sequencing.

Figure 6 Identification and validation of mechanofreqency-sensitive hub genes. (A) The interaction network between proteins coded by the DEGs in profile 0. The nodes represent genes, and the edges represent links between genes. Blue represents downregulated genes. (B) The highest scoring module was extracted by MCODE. (C) The intersection was obtained among modules measured by MCODE, and the top 10 highly connected genes were identified using MCC and Degree in cytoHubba. (D and E) Validation of the expression of the two intersection genes, EYA1 and SALL1, using RTqPCR. *p < 0.05, vs control group.

Discussion

OTM is based on remodeling processes in the periodontal ligament and the alveolar bone. PDLSCs play an important role in mechano-transduction and in promoting periodontal tissue regeneration in OTM.^31,32 It is evident that cyclic tensile force regulates the osteogenic differentiation of PDLSCs.³³ A complex network of signaling molecules regulates the osteoblastic differentiation of PDLSCs under cyclic tension.^5,17,34 The heterogeneity of mechanical force parameters (duration, magnitude, frequency, and others)⁹ led to the heterogeneity of the osteogenic phenotype and gene regulation in different studies. This study focused on the effect of tensile frequency on the osteogenic differentiation of PDLSCs and attempted to elucidate a potential mechanism.

In the current study, we successfully isolated and characterized human PDSLCs. The cells at the 4th to 6th passage were used, whose phenotype was generally believed to be maintained.³⁵ Mechanical tension was applied using the Flexcell tension system, which has been widely used in PDLSCs studies.^2,20 We observed that the cellar reorientation after force loading was similar to that of the previous study,²³ and this could be attributed to the mechano-responsive stress fibers–focal adhesion system.³⁶ The present study showed that cyclic mechanical tension (magnitude: 12% deformation, duration: 10 h) in the range of 0.1–0.7 Hz promoted osteogenic marker genes including Runx2 and COL-I in PDLSCs, and their protein expressions increased with the increasing tensile frequency. Runx2 is an osteogenic lineage commitment specific transcription factor, which binds to the specific cis-acting elements of osteoblasts to promote the transcription and translation of OCN, osteopontin (OPN), bone sialoprotein (BSP), and COL-I.⁸ COL-I, which acts as a template onto which minerals are deposited to form bone matrix,³⁷ is the major constituent of extracellular matrix in the periodontal ligament and bone, and is confirmed to be essential for osteogenesis in response to tension during OTM.³⁸ Consistently, upregulation of RUNX2 and COL-I in response to tension was reported in most studies.⁹ To the best of our knowledge, the present study is the first to reveal the frequency dependence during cyclic tension in enhancing the expression of osteogenic markers within the first 10 h of cyclic tension application, which may inform a new method of accelerating OTM.

Using RNA-seq, we observed that mRNAs in strained PDLSCs were mainly enriched in response to the stimulus process, signal transduction, and relative pathways such as mismatch repair, TNF signaling pathway, and FOXO signaling pathway, which were associated with cell survival and differentiation as well as immune and inflammatory responses.^39,40 The STEM platform was also used to investigate how gene expression profiles change with tensile frequency during the osteoblast differentiation of PDLSCs under cyclic tension. Eight trend profiles were noted as significant. Genes in these profiles were mainly enriched in the metabolic pathways, PI3K-Akt signaling pathway, cytokine-cytokine receptor interaction, and MAPK signaling pathway. The PI3K-Akt signaling pathway has been reported to be involved in the mechanical force-induced osteoblast differentiation of PDLSCs.^41,42 The MAPK signaling pathway also have been found to participate in the mechano-transduction of PDLSCs.^43,44

Genes in profile 0 showed a continuous downward trend from 0.1 Hz to 0.7 Hz, an inverse trend of osteogenic genes, and were mainly enriched for pathways related to an inflammatory response, such as arachidonic acid metabolism and cytokine−cytokine receptor interaction. Cyclooxygenase-2 (COX2) and prostaglandin E synthase (PTGES) participate in the arachidonic acid metabolism pathway. TNF Receptor Superfamily Member 14 (TNFRSF14) and Interleukin 9 (IL9) are involved in the cytokine−cytokine receptor interaction. PTGES is induced by inflammatory mediators.⁴⁵ COX2 is involved in the synthesis of prostaglandin E2 (PGE2), which is a potent pro-inflammatory cytokine, and participates in bone resorption.² The expression of COX2 and PGE2 after mechanical stimulation was previously reported to be correlated with force duration and force magnitude,⁹ and it showed a negative correlation with tensile frequency in the present study. TNFRSF14 is a membrane-bound receptor leading to the induction of proinflammatory genes by activating the NF-κB pathway.⁴⁶ IL-9 plays a role in regulating inflammatory immunity, and it has demonstrated pro-inflammatory activity in several mouse models of inflammation.⁴⁷ The response of periodontal ligament to mechanical stress generated by OTM is known as an aseptic transitory inflammatory process, which is regulated by various cytokines and chemokines.⁴⁸ Proinflammatory cytokines activate matrix metalloproteinases (MMPs), degrade the ECM, and inhibit the expression of COL-I.⁴⁹ Increased osteogenesis is usually accompanied by lower levels of inflammatory cytokines and chemokines.^10,50,51 Accordingly, in the present study, with increasing frequency, the osteogenic commitment increased, and the suppression of pro-inflammatory genes and the relative inflammatory response pathway were observed.

Furthermore, through PPI network screening, we identified two candidate genes, EYA1 and SALL1, which were specifically sensitized to tensile frequency. The result was validated by RTqPCR, which confirmed the decreased expression of EYA1 and SALL1 with increasing frequency of tension stimulation. EYA1 is a conserved critical regulator of organ-specific stem cells.⁵² SALL1 is also considered a stem cell marker.⁵³ The osteoblastic differentiation of PDLSCs increased with increasing tensile frequency; thus, stemness and related genes correspondingly reduced. In view of the high-throughput sequencing and validation after 6 h of tensile force exposure, further studies over a longer period are needed. Whether overexpression of EYA1 and SALL1 would reverse the frequency-dependent trend of the osteogenic differentiation of PDLSCs also deserves further study.

In the present study, the role of tensile frequency on osteogenic commitment of PDLSCs were identified, and the mRNA transcriptomes of PDLSCs during the osteogenic differentiation under cyclic tension with different frequencies were delineated. Frequency series clustering were defined using STEM, and tensile frequency-sensitive genes were identified. This study extends the knowledge about the role of tensile frequency in cyclic tension induced PDLSCs osteogenesis.

Conclusion

The osteoblastic differentiation of PDLSCs under mechanical tensile force is frequency dependent. EYA1 and SALL1 were identified as potential important tensile frequency-sensitive genes, which may contribute to the cyclic tension-induced osteogenic differentiation of PDLSCs in a frequency-dependent manner.