The GSE58389 gene expression dataset was downloaded from the Gene Expression Omnibus (GEO) database. The GSE58389 platform was GPL13607 (Agilent-028004 SurePrint G3 Human GE 8 × 60K Microarray). GSE58389 contained four groups (5 h control, 5 h treated, 24 h control, and 24 h treated). In each group, primary HDFs from ten donors were cultured on BIOFLEX (Ontario, Canada) culture plates and stretched for 5 or 24 h, or were left untreated as controls. Cyclic stretch was applied using the FX-4000T™ Tension Plus™ System (Flexercell International; McKeesport, PA, USA) with 16% elongation at 0.5 Hz in a half sinus regimen. Forty samples were included in the dataset (7).

The microarray datasets at 5 and 24 h were uploaded to the interactive web tool GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/) to screen the DEGs in HDFs between cyclic mechanical stretching and static samples. Next, we assessed the dataset quality of the microarrays and analyzed the differences in gene expression. To maintain the number of DEGs in the original study, only genes with |log (fold change) (logFC)| >0.5 and p < 0.05 were selected as DEGs. To identify the common DEGs, Calculate and Draw Custom Venn Diagrams (http://bioinformatics.psb.ugent.be/webtools/Venn/) was used to obtain the overlapping DEGs in both the 5 h and 24 h groups.

Based on the web-server Metascape (https://metascape.org/gp/index.html#/main/step1), Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed to identify the functions of DEGs in cyclic mechanical stretching and static samples at different times (8). We established the cutoff criteria as p < 0.01, minimum overlap genes = 3, and minimum enrichment factor >1.5. After identifying all statistically enriched GO and KEGG terms, accumulative hypergeometric p values and enrichment factors were calculated and used for filtering. The remaining significant terms were then hierarchically clustered into a tree based on the kappa statistical similarities among their gene memberships. Next, a kappa value of 0.3 was applied as the threshold to cast the tree into term clusters. We then selected a subset of representative terms from these clusters and converted them into a network layout.

Protein–protein interactions (PPIs) occur in all important biological processes in living organisms, such as catalyzing metabolic reactions, DNA replication, DNA transcription, responding to stimuli, and transporting molecules between locations (9). The STRING database (https://string-preview.org, version 11.5) was used to construct the PPI network of DEGs based on the default parameters (required score = 0.4; false discovery rate = 5%) (10).

Cytoscape software version 3.8.2 (Cytoscape Consortium, Boston, MA, USA) was used to visualize the PPI network, as previously described (11). The top 10 hub genes were selected using the plug-in Cytohubba within Cytoscape using the Maximal Clique Centrality method and were then displayed in the extended subnetwork (12). Further GO analyses of hub genes were performed using Metascape.

The NetworkAnalyst web tool (https://www.networkanalyst.ca) was used to predict miRNA gene and gene transcription factor (TF) interactions (13). Visualization of two regulatory hub gene networks was also conducted using Cytoscape.

Descriptive data are presented as the mean ± SD, median, or frequencies and proportions, as appropriate. Statistical significance was set at p < 0.05. All statistical analyses were performed using SPSS version 25.0 (IBM Corporation, Armonk, NY, USA).

The flow diagram created in this study is shown in Figure 1. The boxplots in Figure 2 show the homogeneity of data quality after GEO2R analyses. According to the selection criteria, 669 DEGs were identified at the early stage (5 h), namely, 491 downregulated and 178 upregulated genes, and 249 DEGs were identified at the late stage (24 h), namely, 149 downregulated and 100 upregulated genes. As shown in Figure 3, the Venn diagrams demonstrate overlap in DEGs in both the early and late stages, indicating that there are common DEGs under continuous mechanical stretch.

At the early stage (5 h), the top enriched GO term was “regulation of autophagy” (GO:0010506), and the top KEGG term was “pyrimidine metabolism” (hsa00240). At the late stage (24 h), the top GO term was “sterol biosynthetic process” (GO:0016126), and the top KEGG term was “synaptic vesicle cycle” (ko04721). Genes within each cluster in networks at different stages are shown in Supplementary Table 1. The enrichments and number of DEGs in each statistically significant term are depicted by bubble plots in Figure 4, and the corresponding data are shown in Supplementary Table 2.

Protein-protein interactions (PPIs) were obtained using STRING, visualized using Cytoscape, and further analyzed using the plug-in Cytohubba. Early (5 h), late (24 h), and common PPI networks are shown in Figures 5A, 6A, 7A, respectively. The hub genes at different stages are listed in Table 1. Subnetworks for the top 10 hub genes and their neighbors are illustrated in Figures 5B, 6B, 7B. GO hub gene analyses indicated entirely different terms at different stages (Table 2). Notably, hub genes related to skin development (PKD1, LCE1C, and LCE1D) commonly changed under mechanical stretch.

Hub gene-related miRNAs or TFs were ranked using the Maximal Clique Centrality method. The top three miRNAs or TFs closest to hub genes are listed in Table 3. Notably, certain miRNAs and TFs were predictably related to hub genes in at least two different stages. With regard to miRNAs, hsa-mir-92a-3p was present at all stages, and hsa-mir-193b-3p appeared at the early and late stages. For TFs, POLR2A belonged to early and common hub genes, SMAD5 was part of the early and late hub genes, and MAZ pertained to late and common hub genes. The gene-miRNA regulatory networks are shown in Figure 8. The TF-gene regulatory networks are illustrated in Figure 9.

Gene Ontology (GO) and KEGG analyses can help us elucidate important biological functions and pathways involving DEGs. At the early stage, the top enriched GO term was determined to be “regulation of autophagy.” Consistently, in our previous study, autophagosomes containing mitochondria and cytoplasm were clearly observed in expanded murine scalp fibroblasts at the early stage (14). This suggests that autophagy of HDFs may occur early in response to mechanical stretch. The top KEGG term was “pyrimidine metabolism,” demonstrating that HDFs initiated active nucleotide synthesis. At the late stage, the top GO term was “sterol biosynthetic process,” representing changes in energy metabolism. The top KEGG term was “synaptic vesicle cycle,” representing substance and signal exchanges through repeated exocytosis and endocytosis. However, the GO and KEGG analyses above only provided a general explanation of molecular changes and did not consider the complex interactions between molecules. Therefore, we constructed PPI networks and identified hub genes at different stages.

Following a 5-h cyclic mechanical stretch, representing the early stage, the top 10 hub genes identified primarily participated in DNA and chromatin alterations. Seven histone-related genes were identified and downregulated at this stage: H4 clustered histone 5 (HIST1H4E, H4C5), H4 clustered histone 11 (HIST1H4J, H4C11), H2A clustered histone 17 (HIST1H2AM, H2AC17), H2A clustered histone 8 (HIST1H2AE, H2AC8), H3 clustered histone 15 (HIST2H3A, H3C15), H2A clustered histone 21 (HIST2H2AB, H2AC21), and cluster member H1.4 linker histone (HIST1H1E, H1-4). Histones are fundamental structural components of chromatin (15). Eukaryotic DNA is wound around an octamer of core histones H2A, H2B, H3, and H4. The binding of the linker histone H1 promotes higher-order chromatin organization (15). The marker of proliferation Ki-67 gene (MKI67) encodes a nuclear protein that is required to maintain individual mitotic chromosomes dispersed in the cytoplasm following nuclear envelope disassembly and may be necessary for cell proliferation (16). DNA topoisomerase II alpha (TOP2A) encodes a key decatenating enzyme that alters DNA topology by binding to two double-stranded DNA molecules (17). TOP2A is generally upregulated in proliferating cells (18). However, in another study, fibroblasts, such as mouse NIH 3T3 and 3T6 cells, did not show high TOP2A expression (19). Even under certain stimuli (radiation or drugs), TOP2A expression is downregulated when fibroblasts maintain proliferation activity (20, 21). This may be related to the negative feedback regulation in fibroblasts, which prevents excessive cell proliferation (22). Another upregulated gene, ribonucleotide reductase regulatory subunit M2 (RRM2), encodes one of two non-identical subunits for ribonucleotide reductase, which is necessary for DNA synthesis. RRM2 also functions as a downstream factor of β-catenin as an inhibitor of Wnt signaling (23), and β-catenin activation can stimulate fibroblast proliferation (24, 25).

Following a 24-h cyclic mechanical stretch, representing the late stage, the top 10 hub genes identified primarily participated in cholesterol metabolism. 3-hydroxy-3-methylglutaryl-CoA synthase 1 catalyzes the condensation of acetyl-CoA with acetoacetyl-CoA to form (3S)-hydroxy-3-methylglutaryl-CoA (HMG-CoA), which is then converted by HMG-CoA reductase into mevalonate, a precursor for cholesterol synthesis (26). Next, mevalonate is converted into lanosterol under the action of various enzymes, namely, isopentenyl-diphosphate delta isomerase 1 and squalene epoxidase (26). Lanosterol can then be diverted into either the Bloch pathway, producing cholesterol via desmosterol, or the Kandutsch–Russell pathway, via 7-dehydrocholesterol. Methylsterol monooxygenase 1, 17-betahydroxysteroid dehydrogenase 7, and NAD(P)-dependent steroid dehydrogenase-like protein are involved in these two pathways (26). Furthermore, INSIG1, STARD4, and C14orf1 assist in controlling sterol biosynthesis (27–29). Cholesterol is a critical regulator of lipid bilayer dynamics and plasma membrane organization in eukaryotes (30). The physical properties of the membranes depend on lipid composition; the stiffness and fluidity of the bilayers are essentially determined by the sterol content (31). Various ion channels are modulated by cellular cholesterol and partitioned into cholesterol-enriched membrane rafts (32). After cholesterol depletion, inhibition of stretch-activated cation channels is mediated via actin remodeling and is initiated by the disruption of lipid rafts (31). Thus, cell membrane cholesterol reduction is closely related to changes in cell morphology after stretching. Alterations in metabolic processes play a role in regulating inflammation and extracellular matrix deposition (33). In summary, changes in cholesterol metabolism are an important biological feature of HDFs under stretch conditions at the late stage and may become a potential target to help HDFs adapt more rapidly to changing environments.

Using overlapping DEGs between 5 and 24 h, the consistent top 10 hub genes (seven downregulated and three upregulated) were determined to be different at both stages. This shows that the biological effects of mechanical stretch varied over time. DDX17, XRN2, and STX3 are involved in transcriptional regulation. DDX17 encodes an important context-dependent transcriptional regulator that promotes cell growth by interacting with estrogen receptors (34). XRN2 (DHP1 in the yeast genus Schizosaccharomyces), an upregulated gene, triggers premature transcription termination and nucleates heterochromatin to promote meiotic gene silencing (35). STX3, another upregulated gene, also acts as a transcriptional regulator; inhibition of endogenous STX3 expression alters cellular genes and promotes cell proliferation (36). EXOC7, CASKIN1, RAVER1, and TCP1 are involved in cytoskeleton rearrangement. EXOC7, in addition to functioning in exocytosis, regulates actin at the leading edges of migrating cells, thereby coordinating cytoskeleton and membrane trafficking during cell migration (37). CASKIN1, a scaffold protein, regulates actin filaments (38). RAVER1 interacts with the cytoskeletal proteins actinin and vinculin (39). TCP1, another upregulated gene, inhibits the transformation of fibroblasts into myofibroblasts, thus adjusting for the morphological changes caused by mechanical stretch (40). Previous studies have also reported that changes in the cytoskeleton are important results of mechanical signals and mediate the synthesis of the extracellular matrix triggered by mechanical stretch (41). Overall, hub genes related to transcriptional regulation and cytoskeleton rearrangement are potential targets for promoting HDF regeneration and alignment.

Importantly, three consistent DEGs (LCE1D, LCE1C, and PKD1) were downregulated and identified as hub genes. These are all closely related to skin development, according to the GO analysis. Mechanical stretch is believed to regulate intracellular calcium homeostasis, resulting in changes in a series of downstream signaling pathways (42). LCE1D and LCE1C may be downregulated in response to extracellular calcium alterations (43). However, previous studies on the function of late cornified envelope proteins were mainly conducted in keratinocytes and not in fibroblasts. Thus, the mechanisms by which LCE1D and LCE1C are involved in fibroblast differentiation and growth may be new targets for future research. PKD1 encodes an integral membrane protein involved in the regulation of mechanotransduction signaling (44). The component of heteromeric calcium-permeable ion channels formed by PKD1 and PKD2 is activated by the interaction between PKD1 and a Wnt family member, such as WNT3A or WNT9B (45). PKD1 induces cell migration by regulating rearrangements and cell-cell mechanical adhesion (46), inhibits cell apoptosis through a PKR-eIF2α pathway (47), and regulates the cell cycle by inhibiting DNA binding (48). Therefore, downregulation of PKD1 may be an important mechanism in HDF-sensing mechanical stretching and controlling cell growth and differentiation.

Predicting gene-miRNA and TF-gene regulatory networks at different stages can provide a reference for subsequent interventions. Two miRNAs (hsa-mir-92a-3p and hsa-mir-193b-3p) and three TFs (POLR2A, SMAD5, and MAZ) were identified in at least two stages. Hsa-mir-92a-3p is extensively involved in the regulation of cellular proliferation and angiogenesis but has cell type-specific effects in vivo (49, 50). Hsa-mir-193b-3p is widely expressed in normal human tissues and displays antiproliferative effects associated with its functions in different cell differentiation processes (51). POLR2A encodes the largest subunit of RNA polymerase II, which is responsible for synthesizing messenger RNA in eukaryotes. Multiple pathways have been shown to regulate cell proliferation by affecting POLR2A expression (52–54). SMAD5 proteins are intracellular signaling molecules that mediate canonical bone morphogenetic protein pathways and are involved in tissue regeneration (55). MAZ is a zinc finger transcription factor that activates the expression of tissue-specific genes and represses the expression of the c-myc proto-oncogene. Previous data indicate that MAZ is a growth suppressor protein that affects the cell cycle in fibroblasts (56). Therefore, MAZ knockdown may be effective in continuously activating HDFs. The above regulatory network provides a reference for further intervention approaches.