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 downregulated and upregulated 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 increased. 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.

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

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 TM Tension Plus TM 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).

Screening of DEGs
The flow diagram created in this study is shown in Figure 1. 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.

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.

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.

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.

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 downregulated at this stage: H4  (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).

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, 17betahydroxysteroid 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)(28)(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.

Common Hub Genes
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.

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)(53)(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 strengthened 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.