Identification of Common Hub Genes in Human Dermal Fibroblasts Stimulated by Mechanical Stretch at Both the Early and Late Stages.

Background: Mechanical stretch is vital for soft tissue regeneration and development and is utilized by plastic surgeons for tissue expansion. Identifying the common hub genes in human dermal fibroblasts (HDFs) stimulated by mechanical stretch at different stages will help elucidate the mechanisms involved and improve the efficiency of tissue expansion.
Methods: A gene expression dataset (GSE58389) was downloaded from the Gene Expression Omnibus database. Differentially expressed genes (DEGs) in HDFs between cyclic mechanical stretching and static samples were identified at 5 and 24 h. Common DEGs overlapped in both the 5 h and 24 h groups. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed to determine the functions of the DEGs. Protein-protein interaction networks were constructed using the STRING database. The top 10 hub genes were selected using the plug-in Cytohubba within Cytoscape. The regulatory network of hub genes was predicted using NetworkAnalyst.
Results: A total of 669 and 249 DEGs were identified at the early (5 h) and late stages (24 h), respectively. Of these, 152 were present at both stages and were designated as common DEGs. The top enriched GO terms were "regulation of autophagy" at the early stage, and "sterol biosynthetic processes" at the late stage. The top KEGG terms were "pyrimidine metabolism" and "synaptic vesicle cycle" at the early and late stages, respectively. Seven common DEGs [DEAD-box helicase 17 (DDX17), exocyst complex component 7 (EXOC7), CASK interacting protein 1 (CASKIN1), ribonucleoprotein PTB-binding 1 (RAVER1), late cornified envelope 1D (LCE1D), LCE1C, and polycystin 1, transient receptor potential channel interacting (PKD1)] and three common DEGs [5'-3' exoribonuclease 2 (XRN2), T-complex protein 1 (TCP1), and syntaxin 3 (STX3)] were shown to be upregulated and downregulated hub genes, respectively. The GO terms of the common hub genes were "skin development" and "mRNA processing." After constructing the regulatory network, hsa-mir-92a-3p, hsa-mir-193b-3p, RNA polymerase II subunit A (POLR2A), SMAD family member 5 (SMAD5), and MYC-associated zinc finger protein (MAZ) were predicted as potential targets in both stages.
Conclusion: At the early stage, there were clear changes in gene expression related to DNA and chromatin alterations; at late stages, gene expression associated with cholesterol metabolism was suppressed. Common DEGs related to skin development, transcriptional regulation, and cytoskeleton rearrangement identified in both stages were found to be potential targets for promoting HDF growth and alignment under mechanical stretch.

Introduction

Mechanical stretch is a force that is essential for the regeneration and development of skin and other soft tissues (1). Naturally, the growth and plasticity of skin and soft tissues can be observed when other structures grow within the body. For instance, embryo growth causes the abdominal skin to expand during pregnancy, and skull growth causes the scalp to expand in the fetus. Clinically, plastic and reconstructive surgeons utilize the extra skin originating from tissue expansion to reconstruct scars or repair bodily defects. Compared with traditional methods, such as skin graft and flap transfer, tissue expansion can regenerate the area of the donor tissue and substantially reduce the deformities of the donor site. In addition, expanded flaps effectively match the color and texture of the recipient area (2). Therefore, identifying hub genes and determining molecular changes in skin and soft tissue stimulated by mechanical stretch will help elucidate the biological behaviors of the cells and improve the efficiency of tissue expansion.

Recent studies have also investigated the molecular changes in expanded skin induced by mechanical stretch in animal models (3). By conducting a transcription analysis on expanded skin biopsies of pigs, Ledwon et al. (4) revealed changes in immune response activation, cell metabolism, and processes related to muscle contraction and cytoskeleton organization. Using a skin-stretched mouse model, Aragona et al. (5) showed that mechanical stretch created a transient bias in the renewal activity of epidermal stem cells and a second subpopulation of basal progenitors committed to differentiation.

However, the structure and mechanical properties of the skin in humans are different from those in other animals (6). Therefore, comprehensive molecular changes in human skin stimulated by mechanical stretch are not entirely understood. Moreover, it is difficult to collect specific human skin specimens stimulated by mechanical stretch under standard conditions and time periods. Fibroblasts are the main cell type in the dermis. Thus, human dermal fibroblasts (HDFs) are suitable cell models for studying molecular changes under mechanical stretch at distinct time periods. In this study, we focused on identifying common hub genes and their potential regulatory networks in HDFs stimulated by mechanical stretch at different stages.

Materials and Methods

Microarray Data

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).

Identification of Differentially Expressed Genes

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.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment Analyses

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.

Construction of the Protein–Protein Interaction Network

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).

Hub Gene Analysis

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.

Regulatory Network of Hub Genes

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.

Statistical Analysis

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).

Results

Screening of DEGs

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 upregulated and 178 downregulated genes, and 249 DEGs were identified at the late stage (24 h), namely, 149 upregulated and 100 downregulated 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.

Figure 1

Flow diagram.

Figure 2

Data quality boxplot. (A) At 5 h and (B) at 24 h.

Figure 3

Venn diagram. (A) Upregulated and (B) downregulated genes.

GO and KEGG Pathway Enrichment Analyses

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.

Figure 4

GO enrichment and KEGG analyses showed changes in the molecular response and signaling pathway of HDFs under mechanical stretch. (A) GO enrichment analysis of the DEGs at 5 h and (B) at 24 h. (C) KEGG analysis of the DEGs at 5 h and (D) at 24 h. The gene ratio indicates the number of DEGs associated with the GO term divided by the total number of DEGs. The size of the dots represents the number of DEGs associated with the GO term, and the color represents the negative value of log10 of p value. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; HDFs, human dermal fibroblasts; DEGs, differentially expressed genes.

PPI Networks and Hub Genes

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.

Figure 5

PPI network at 5 h. (A) Entire network. Upregulated and downregulated genes are represented by red and blue backgrounds, respectively. (B) Subnetwork for the top 10 hub genes and neighbors. The backgrounds of hub genes are colored from red to yellow by rank, and the backgrounds of neighboring genes are light blue. PPI, protein-protein interactions.

Figure 6

PPI network at 24 h. (A) Entire network. Upregulated and downregulated genes are represented by red and blue backgrounds, respectively. (B) Subnetwork for the top 10 hub genes and neighbors. The backgrounds of hub genes are colored from red to yellow by rank, and the backgrounds of neighboring genes are light blue. PPI, protein-protein interactions.

Figure 7

Common PPI network. (A) Entire network. Upregulated and downregulated genes are represented by red and blue backgrounds, respectively. (B) Subnetwork for the top 10 hub genes and neighbors. The backgrounds of hub genes are colored from red to yellow by rank, and the backgrounds of neighboring genes are light blue. PPI, protein-protein interactions.

Table 1

Early (5 h), late (24 h), and common top 10 hub genes.

Stage	Top 10 hub genes*
	Upregulated	Downregulated
Early	HIST1H4E (H4C5), HIST1H4J (H4C11), HIST1H2AM (H2AC17), HIST1H2AE (H2AC8), HIST2H3A (H3C15), HIST2H2AB (H2AC21), HIST1H1E (H1-4), MKI67	TOP2A, RRM2
Late	None	HMGCR, HMGCS1, SQLE, MSMO1, HSD17B7, IDI1, NSDHL, INSIG1, STARD4, C14orf1 (ERG28)
Common	DDX17, EXOC7, CASKIN1, RAVER1, LCE1D, LCE1C, PKD1	XRN2, TCP1, STX3

*The genes were ranked using the Maximal Clique Centrality method.

Table 2

GO* terms of hub genes.

Stage	Term description	Gene	p value
Early	DNA packaging	HIST1H1E, TOP2A, HIST1H4J, HIST1H4E, HIST2H3A	8.30E−09
	Chromatin silencing	HIST1H1E, HIST1H2AE, HIST1H2AM, HIST2H2AB	1.02E−08
Late	Sterol biosynthetic process	HMGCR, HMGCS1, IDI1, INSIG1, MSMO1, SQLE, C14orf1, NSDHL, HSD17B7, STARD4	1.04E−23
	Steroid biosynthetic process	HMGCR, HMGCS1, IDI1, INSIG1, MSMO1, SQLE, C14orf1, NSDHL, HSD17B7, STARD4	6.24E−23
	Isoprenoid biosynthetic process	HMGCR, HMGCS1, IDI1, INSIG1	1.05E−07
	Cellular ketone metabolic process	HMGCR, INSIG1, STARD4, SQLE, MSMO1	4.59E−05
Common	Skin development	PKD1, LCE1C, LCE1D	9.21E−05
	mRNA processing	DDX17, XRN2, RAVER1	6.08E−04

Gene-miRNA and TF-Gene Regulatory Network

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.

Table 3

Top 3 miRNAs or TFs closest to hub genesa.

Stage	miRNAs	TFs
Early	hsa-mir-34a-5p, hsa-mir-92a-3p, hsa-mir-193b-3pb	SIN3A, POLR2A, SMAD5c
Late	hsa-mir-335-5p, hsa-mir-193b-3p, hsa-mir-92a-3pd	MAZ, CEBPA, SMAD5e
Common	hsa-mir-92a-3p, hsa-mir-615-3p, hsa-mir-331-3p	ETV4, KLF9, POLR2Af

.

Score of hsa-mir-186-5p and hsa-mir-17-5p was equal to the last one (score = 4).

Score of MLLT1 was equal to the last one (score = 7).

Score of hsa-mir-192-5p were equal to the last one (score = 3).

Score of CHD1, MYNN, PPARG, CREB3L1, ZFP37, FOXM1, KDM5B, PHF8, IRF4, SAP30 and EGR1 was equal to the last one (score = 3).

Score of MAZ, HIC1, INSM2, KLF8 and ELK1 was equal to the last one (score = 4).

Bold value means that one miRNA or TF appears in at least two stages.

Figure 8

Gene-miRNA regulatory network. (A) Early stage (5 h); (B) late stage (24 h); (C) common stage. The background colors of genes and miRNAs are red and light blue, respectively. The sole connection of the gene-miRNA is hidden.

Figure 9

TF-gene regulatory network. (A) Early stage (5 h); (B) late stage (24 h); (C) common stage. The background colors of genes and miRNAs are red and light blue, respectively. The sole connection of the TF-gene is hidden. TF, transcriptional factor.

.

Discussion

Mechanical stretch-induced skin regeneration provides sufficient material for wound repair and organ reconstruction. However, the mechanisms by which mechanical stretch promotes skin regeneration remain predominantly unknown. Fibroblasts are one of the main cell types present in the skin; thus, determining the effects of mechanical stretch on HDFs may help to improve our understanding of the mechanisms and assist in promoting expanded skin regeneration. Reichenbach et al. (7) compared stretched HDFs with unstretched HDFs to identify DEGs using mRNA microarray data. However, because of unknown reasons, these authors did not further study or discuss the functions of these DEGs. More crucially, the GO and KEGG pathway enrichment analyses, PPI network construction, hub gene analysis, and regulatory hub gene network were absent in the previous work (7). To solve this problem, we performed this study based on the uploaded mRNA microarray data and previous reports (7). Herein, we provide the DEGs in mechanically stretched HDFs and unstretched HDFs, and the common DEGs present in both the early and late stages. More importantly, we presented novel results of GO and KEGG pathway enrichment analyses, PPI network construction, hub gene analysis, and a regulatory hub gene network. These may assist in identifying genes involved in promoting fibroblast proliferation and skin regeneration.

GO and KEGG Analyses

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.

Early Hub Genes

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 upregulated 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). The simultaneous increase in histones and Ki-67 expression in HDFs under stretching indicates that active DNA packaging was involved in accelerated cell proliferation. 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 this study, it was slightly downregulated (|log[FC]| = −0.511). 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 downregulated 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). Therefore, we hypothesized that RRM2 downregulation is one of the mechanisms of Wnt/β-catenin activation.

Late Hub Genes

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-beta-hydroxysteroid 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). All of the genes mentioned above were downregulated. 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). Ledwon et al. (4) found that, in a 24-h porcine tissue expansion model, a vast majority of DEGs were associated with metabolic processes, such as catabolism of lipids and organic acids, and were downregulated. Similarly, in our study, sterol metabolism-related DEGs were identified after a 24-h cyclic mechanical stretch and were downregulated. This demonstrates the importance of metabolic transitions in the response of HDFs to mechanical stretch in adapting and re-establishing tissue homeostasis (4). 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.

Common Hub Genes

Using overlapping DEGs between 5 and 24 h, the consistent top 10 hub genes (seven upregulated and three downregulated) 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), a downregulated gene, triggers premature transcription termination and nucleates heterochromatin to promote meiotic gene silencing (35). STX3, another downregulated 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 downregulated 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 upregulated 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 upregulated 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, upregulation of PKD1 may be an important mechanism in HDF-sensing mechanical stretching and controlling cell growth and differentiation.

Regulatory Network

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.

Conclusion

At the early stage, DNA and chromatin alterations were clearly observed; at the late stage, cholesterol metabolism was repressed to adapt to the changing environment. In common stages, genes related to transcriptional regulation, cytoskeleton rearrangement, and skin development were found to be potential targets to promote HDF growth and alignment under mechanical stretch. The hub genes and their regulatory networks were different during different periods.

We found that hub genes and their functions were considerably different during distinct periods. Bioinformatics analysis of common hub genes in stretched HDFs present at different stages provides direction for subsequent research. Thus, our findings can provide a reference for the precise regulation of HDF behavior in response to mechanical stretch. In the future, PKD1, among other common hub genes, maybe a new target to promote mechanical stretch-induced skin regeneration.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

CD and WL: designing the study, conducting the study, acquiring and analyzing data, and writing the manuscript. YZ, YS, JD, ZH, and TW: the concept of the study, analyzing data, and editing the manuscript. YZ and XM: concept of the study, designing the study, and writing/editing the manuscript. All authors contributed to the article and approved the submitted version.

Funding

Funding from the National Natural Science Foundation of China (82172229 and 81971851), the Natural Science Basic Research Plan in Shaanxi Province of China (2022JM-600), and the Foundation of Xijing Hospital Grant (XJZT21CM33) was received for this study.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.