Integrated PPI- and WGCNA-retrieval of hub gene signatures for soft substrates inhibition of human fibroblasts proliferation and differentiation

Fibroblasts (FBs) are the most important functional cells in the process of wound repair, and their functions can be activated by different signals at the pathological site. Although wound repair is associated with microenvironmental stiffness, the effect of matrix stiffness on FBs remains elusive. In this study, TGF-β1 was used to mimic the fibrotic environment under pathological conditions. We found that the soft substrates made FBs slender compared with tissue culture plastic, and the main altered biological function was the inhibition of proliferation and differentiation ability. Through PPI and WGCNA analysis, 63 hub genes were found, including GADD45A, CDKN3, HIST2H3PS2, ACTB, etc., which may be the main targets of soft substrates affecting the proliferation and differentiation of FBs. Our findings not only provide a more detailed report on the effect of matrix stiffness on the function of human skin FBs, but also may provide new intervention ideas for improving scars and other diseases caused by excessive cell proliferation, with potential clinical application prospects.


INTRODUCTION
Pathological scars are pathological products formed during wound healing [1][2][3]. As a typical fibroproliferative disease [4], its pathological changes are mainly caused by massive proliferation of fibroblasts (FBs) and excessive secretion and deposition of collagen [5]. As the core of scar repair, FBs play an important role in skin wound repair [6], and they may be the key to achieving scar-free healing [7]. After trauma, FBs undergo chemotaxis, proliferate and differentiate into myofibroblasts under the control of cytokines [8]. One of the important reasons for the occurrence of pathological scars is the abnormal overproliferation of FBs. For example, the imbalance of microenvironmental homeostasis can lead to the abnormal proliferation of FBs, which leads to the formation of keloids [9]. In conclusion, the fate of FBs determines the end state of the scar. Understanding and controlling the biological behavior of FBs is the basis and key to promote wound healing and prevent scarring.
Wound repair is a complex and orderly biological process (BP), which is the result of the interaction and mutual influence of cells, growth factors and extracellular matrix [10]. Under normal conditions, the matrix is involved in maintaining the quiescent state of FBs, and with the increase of matrix stiffness, the proliferation of FBs is gradually activated, which promotes scarring [11]. Earlier studies have found that the stiffness of scar tissue is significantly higher than that of normal tissue, which is the result of scar formation. More and more studies have proved that the changes in the stiffness of the local microenvironment of injury can directly determine the repair and regeneration ability of the wound, and play an important role in the pathogenesis of pathological scars [12][13][14]. If the wound continues to be stimulated by a large mechanical force, the myofibroblasts continue to proliferate and stimulate the differentiation of myofibroblasts, and the wound will form a pathological scar [15,16]. At present, most of the treatment methods for pathological scars are based on the principle of inhibiting the inflammatory response. The understanding of the effect of mechanical force signals on scar formation is not deep enough. Therefore, it has become a valuable question whether FBs can be modulated by changing mechanical force or modifying mechanical transduction signals to affect scar formation.
The traditional culture of fibroblast is on tissue culture plastic (TCP) (GPa), the stiffness of its cultural environment is much higher than the physiological environment of FBs [17], this will inevitably affect the properties of FBs, therefore, we constructed the culture conditions with similar stiffness to the normal skin microenvironment (700Pa-1120Pa) [18,19], explored the mechanical signal-related mechanism of soft substrates affecting the fate of FBs, and tried to provide new ideas and therapeutic targets for the development of pathological scar treatment strategies.
Transforming growth factor β1 (TGF-β1) is a key factor in promoting tissue and organ fibrosis, it regulates the phenotype and function of fibroblasts [20,21]. And it can also induce the formation of myofibroblasts while inhibit their apoptosis, and promote the deposition of extracellular matrix with collagen as the main component [21,22]. In severely injured tissue and hypertrophic scar tissue, compared with normal skin tissue, the expression level of TGF-β1 is significantly up-regulated [23,24], which is considered to be closely related to the formation of pathological scar. Current research suggests that TGF-β1 should be added to in vitro experiments so that it can simulate the fibrotic environment under pathological conditions to a certain extent [24], and the results are more meaningful to detect the effect of matrix stiffness on the fate of fibroblasts under this condition.

Soft substrates altered the morphology of human foreskin fibroblasts (hFFs)
To explore the effect of soft substrates on FBs, we selected healthy, easily cultured and accessible hFFs, which was developed on either TCP or soft substrates, the elastic modulus of the soft substrates was reported to be 500-1000Pa [25]. At the meantime, the hFFs have been divided into different group based on whether TGF-β1 is added. The cell morphology of the four groups was observed and the aspect ratio was counted, and it was found that after culturing hFFs on TCP and soft substrates with or without TGF-β1, the soft substrates made the boundary of hFFs more distinct, three-dimensional, and slender ( Figure 1A), with significant increased aspect ratio ( Figure 1B), indicating that the soft substrates had a certain effect on hFFs, and which BPs are specifically affected need further analysis.

Differential expression analysis
To explore the BPs and signaling pathways of hFFs that can be altered by soft substrates, we performed mRNA transcriptome sequencing of hFFs cultured on TCP and soft substrates for 24 h, and found that the number of genes co-expressed by soft and TCP was 112,698, the gene expression numbers specific to soft and TCP were 11718 and 10103, respectively ( Figure 2A). The number of differentially expressed genes was large, which was not conducive to subsequent analysis. Therefore, the differential expression was screened. It was considered that |log2foldchange|≥2 was a significantly differentially expressed gene and was used for subsequent analysis. Among them, the number of up-regulated differentially expressed genes was 349, and the number of downregulated differentially expressed genes was 516 ( Figure  2B). The analysis of the inter-sample repeatability of these genes also found that the repeatability within the group was good, meeting the conditions for subsequent analysis ( Figure 2C).
Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of up-regulated and down-regulated differentially expressed genes (DEGs) showed that the pathways promoted by soft substrates included: pathways in cancer, shigellosis, pertussis, chemokine signaling pathway, focal adhesion, bacterial invasion of epithelial cells, regulation of actin cytoskeleton, chagas disease (American trypanosomiasis), etc. The inhibited pathways included: cell cycle, regulation of actin cytoskeleton, proteoglycans in cancer, microRNAs in cancer, focal adhesion, ubiquitin mediated proteolysis, pathogenic Escherichia coli infection, hippo signaling pathway, etc. ( Figure 2D).
Gene Ontology (GO) analysis of up-regulated and downregulated DEGs showed that the BP promoted by soft substrates included: extracellular matrix organization, extracellular structure organization, negative regulation of cellular process, negative regulation of BP, anatomical structure development, multicellular organism development, etc.; Promoted cellular component (CC) included: proteinaceous extracellular matrix, extracellular matrix, extracellular region part, extracellular space, extracellular region, nucleosome, etc.; Promoted molecular function (MF) included: protein binding, protein dimerization activity, protein heterodimerization activity, signaling receptor binding, metalloendopeptidase activity, RNA polymerase II activating transcription factor binding, etc., BP inhibited by soft substrates included developmental process, anatomical structure development, CC organization, CC organization or biogenesis, cell division, cell cycle, cell proliferation, etc.; inhibited CC included cytoplasm, intracellular non−membrane−bounded organelle, non−membrane− bounded organelle, intracellular part, adherens junction, anchoring junction etc., inhibited MF included protein binding, enzyme binding, cell adhesion molecule binding, kinase binding, protein kinase binding, cytoskeletal protein binding, etc. ( Figure 2E). The above results indicated that the effects of soft substrates on hFFs were mainly related to proliferation, differentiation, cytoskeleton and cell-matrix junctions.

Soft substrates inhibited proliferation and differentiation of hFFs
In order to explore the main biological effects of soft matrix on hFFs, we carried out gene set enrichment analysis (GSEA) and found that among the GSEA-enriched downregulated entries, cell cycle was the first one, and the enriched genes involved in this entry were shown in the heatmap ( Figure 3A). To further validate this result, Ki67 immunofluorescence staining was performed on TCP and soft with or without TGF-β1 ( Figure 3B), the results showed that the positive rate of Ki67 in the soft group was much lower than that in the TCP group with or without TGF-β1 ( Figure 3C). The cell cycle of the four groups was examined ( Figure 3D), and it was found that the soft substrates could reduce the PI of hFFs ( Figure 3E), suppress the proportion of cells in G2 phase and increase the proportion of cells in G1 phase ( Figure 3F), indicating that soft substrates inhibited hFFs proliferation.
At the same time, the GSEA analysis also had enriched entries of vascular smooth muscle contraction that interest us, indicating that it had an impact on contractility, including the myofibroblast marker protein α-SMA (also known as ACTA2) [26] ( Figure 4A), it was also considered to be one of the most stable cytoskeletal components of myofibroblasts [27], it was found that the soft substrates did significantly downregulate its expression ( Figure 4B). Periostin regulates myofibroblast differentiation and is persistently overexpressed in abnormal scars and other benign fibrous tissues proliferating with fibroblasts [28], while soft substrates Inhibit its expression, also in the presence of TGF-β1 stimulation ( Figure 4C), it showed that the soft substrates may inhibit the differentiation and contraction of hFFs, and combined with the soft substrates inhibiting the proliferation of hFFs, these phenomena have certain significance for the treatment of skin scarring, sclerosis, fibrosis and other diseases related to the proliferation of FBs.

Protein-protein interaction network (PPI) network analysis
In order to explore the PPI network that can significantly change DEGs and select the major sub-networks, we used the string database to construct a PPI network of 698 DEGs and 9602 relations ( Figure  5A). Interactions between key genes across the network were determined using Cytoscape plugins (MCODE and CytoHubba). Among them, MCODE obtained the top two clusters, there were 56 nodes and 2160 edges in cluster 1, the enrichment score was 39.273, and there were 19 nodes and 192 edges in cluster 2, the enrichment score was 10.667 ( Figure 5B), they were identified from MCODE according to a scoring system. In addition, after importing the PPI network of Figure  5A into another plug-in CytoHubba, 114 key genes were identified by the MCC calculation method, with a total of 1536 edges ( Figure 5C). In the above network, the circles represented down-regulated genes, the squares represented up-regulated genes, the blue to red color gradient of node full color corresponds to the log2foldchange in the range of -21.03 to 21.03, and the node size (26.02 to 61.80) corresponds to the degree (1.57 to 162) in Figure 5A network ( Figure 5D), these genes obtained by MCODE and MCC algorithms were important in DEGs, but the number of genes was large, and further screening was needed.

WGCNA and meaningful module identification
To further screen hub genes, we performed Weighted gene coexpression network analysis (WGCNA) on the transcriptional gene expression profile, sample The legend of networks. The round represents down-regulated DEGs, the square represents up-regulated DEGs, and the size of the node graph represents the degree for (A), which denotes the number of nodes connected to each node. The colors of the nodes indicate the size of log2 (fold change). The higher and lower the expression is, the redder and bluer it is, respectively. AGING clustering was performed, and the clustering dendrogram revealed no obvious outliers ( Figure 6A), therefore, all samples could be analyzed in the next step. Then, the value of the soft threshold power β was calculated before constructing the gene co-expression network, the scale-free topology threshold of the network was 0.85, and when the soft threshold power was 14, the mean connectivity was close to 0, therefore, β = 14 was chosen to build a hierarchical clustering tree ( Figure 6B). Ultimately, 22 modules were identified based on average hierarchical clustering and dynamic tree clipping ( Figure 6C). Eigengene adjacency heatmap was used to illustrate the relationship between eigengenes and phenotypic traits ( Figure 6D).
The direct relationship between each module and the sample in the transcriptome data, including TCP and soft groups, was analyzed, the results showed that among the modules, the blue module had the strongest correlation with soft (r = 0.99, P = 3e-04), including 3215 genes, the turquoise module had the strongest correlation with TCP (r = 0.96, P = 0.002), including 3748 genes ( Figure  6E), indicating that the genes of blue and turquoise modules might play a role in the effect of soft substrates on hFFs important role, therefore, the blue and turquoise modules were considered as key modules for further analysis.

Functional analysis and PPI network construction of hub genes
In order to find more critical hub genes, as well as their relationships and biological functions, we first intersected the 114 genes obtained by the MCC algorithm in CytoHubba, the 75 genes obtained by cluster1 and 2 in MCODE, and the 6963 genes in the blue and turquoise modules of WGCNA, and obtained 63 hub genes ( Figure 7A) (Table 1), a PPI network was constructed for them to show their relationship ( Figure  7B). KEGG and GO analysis of 63 hub genes found that cell cycle and cell division were the most significant entries ( Figure 7C, 7D) (Tables 2, 3), indicating that the most significant effect of soft substrates on hFFs was cell proliferation, and the PPI networks of genes within the entries constructed to represent their relationship ( Figure 7E, 7F). In the above network, circles represented down-regulated genes, squares represented up-regulated genes, the blue to red color gradient of node full color corresponded to log2foldchange in the range of -21.03 to 21.03, and node size (20.24 to 33.33) corresponded to the value of degree (5 to 58) in the network ( Figure 7G), among them, α-SMA (ACTA2) was one of the 63 hub genes, indicating that the main effect of soft substrates on hFFs was to inhibit the proliferation and differentiation of hFFs. At the same time, we further verified the expression of the predicted top-ranked hub genes. Compared with TCP, soft substrates significantly inhibited the expression of CDKN3 and ACTB. In the presence of TGF-β1, soft substrates significantly inhibited the expression of GADD45A, CDKN3 and HIST2H3PS2 ( Figure 7H), which indicated that our predicted hub genes also have certain reliability in the scar hyperplasia model in vitro. After a series of analyses, these hub molecules might become therapeutic targets for hypertrophic scars.

DISCUSSION
After trauma, injury or surgery, FBs in the wound proliferate through mitosis, and their proliferation plays an important role in wound repair in the body. Excessive proliferation of FBs can cause pathological scars and other related diseases, early non-surgical intervention for hypertrophic scarring is the mainstream direction of future development [29].
The increased stiffness of the extracellular matrix is not only the result of fibrosis, but also plays a key role in the processes affecting fibroblast proliferation and matrix synthesis [30]. Studies have shown that the stiffness and other properties of the extracellular matrix have a certain impact on the repair function of FBs, which has potential significance for mesenchymal cell therapy [31]. In fact, many experimental studies have demonstrated that FBs adhere and proliferate more firmly on stiff substrates than soft substrates [30,32,33]. In a study of a mouse model of pulmonary fibrosis, the matrix under conditions of normal stiffness is involved in maintaining the quiescent state of FBs [30,34]. With the increase of matrix stiffness, the proliferation of FBs was gradually activated. In a similar study on mouse cardiac FBs, the experimental results also showed that soft matrix is not conducive to the proliferation and adherence of FBs [35]. In conclusion, matrix stiffness can affect the proliferative activity of FBs [33,36]. However, the specific mechanism by which extracellular matrix stiffness affects fibroblast proliferation is still unclear.
In this study, it was found through transcriptome analysis that the effect of soft substrates on hFFs mainly inhibited proliferation and differentiation, mainly in that it increased the G1 phase and inhibited the G2 phase of the cell cycle, and inhibited the expression of α-SMA and periostin [28]. In the KEGG pathway analysis, we found that there was enrichment of biological functions in the down-regulated pathways of regulation of actin cytoskeleton and cell cycle, and these pathways were closely related to cell proliferation. In addition, we noticed significant enrichment of actin cytoskeleton and focal adhesion related signaling pathways in both up-regulation and down-regulation. Cytoskeleton is the basic structure of intracellular mechanical force signal transduction [37], and focal adhesion is the structural connecting unit between ECM and actin skeleton [38,39]. This result suggests that the regulation of ECM stiffness on the fate of FBs relies on the classic pathway of cytoskeleton as a signal transduction pathway. It was particularly noteworthy that we found the existence of Hippo signaling pathway in the downregulation pathway with the top ranking of KEGG analysis results. Hippo pathway has been identified as the key pathway of mechanical stimulation signal conversion and transmission, which is closely related to the proliferation, differentiation and other functions of cell populations [40,41]. TAP/TAZ is the main effector of Hippo signaling pathway, and it is also the decisive factor for extracellular matrix stiffness to determine the fate and function of fibroblasts [42]. In recent years, studies have proposed that mechanical stimulation signal and Hippo signaling pathway are jointly involved in the regulation of YAP/TAZ activity [43,44]. Our research results confirm this view to a certain extent. The specific mechanism of this regulation is still unclear, and our results may provide key node genes for revealing the relationship between the two. AGING In the GO analysis, BP showed that the enriched upregulated pathways included negative regulation of BP, negative regulation of cellular process, and the down-regulation pathways mainly include cell cycle process, cytoskeleton organization, cell proliferation, cell cycle, and cell division, these results all point to the The legend of networks. The round represents down-regulated DEGs, the square represents up-regulated DEGs, and the size of the node graph represents the degree, which denotes the number of nodes connected to each node. The colors of the nodes indicate the size of log2 (fold change). The higher and lower the expression is, the redder and bluer it is, respectively. (H) Quantitative reverse transcription PCR (qRT-PCR) analysis comparing GADD45A, CDKN3, HIST2H3PS2 and ACTB expression levels in soft and TCP with or without TGF-β1. * P < 0.05; ** P < 0.01; *** P < 0.001 (mean, n = 3).   proliferation function of cells. These results also support that the changes of cytoskeleton may be the reason why soft substrates inhibited FBs proliferation [45,46].
In order to explore the hub genes affected by soft substrates on hFFs, we analyzed the WGCNA and PPI networks, and found 63 hub genes by taking the intersection. The most important functional enrichment of these hub genes was cell cycle and cell division. There were 12 genes in the cell cycle, among which the log2foldchange of MCM2, GADD45A, CCNB1, CHEK1, and PKMYT1 ranked high. The existing literature confirms that the functions of the above genes are all related to cell proliferation [47][48][49][50][51]. There were 25 genes involved in cell division, among which LIG1, CCNB1, AURKA, RCC1 and CDC25A had the highest log2foldchange. These genes have also been proved to be involved in the regulation of cell proliferation [49,[52][53][54][55]. In addition, 63 hub genes contain α-SMA, which is an important marker of FBs differentiation into myofibroblasts [16,56].
Considering the critical role of this differentiation process in wound healing and scar formation, matrix stiffness may have an impact on scar outcome by regulating differentiation function.
In summary, our findings demonstrate the effect of soft substrates on the proliferation and differentiation of FBs even in the presence of TGF-β1, and it is expected to apply soft substrates in the treatment of pathological scars, and genes such as GADD45A, CCNB1, MCM2, CDC25A have the potential to become new therapeutic targets. Our research also provides new clues for the transformation of hypertrophic scars to atrophic scars, and opens up new ways to explore how to reduce scars or no scars after wound healing.

Preparation of the collagen I
Firstly under the condition of ice bath, the amount of solution as described in the following

Immunofluorescence staining
FBs were cultured in a 24-well plate. First, cells were fixed with 4% paraformaldehyde for 10 min at 4°C. Subsequently, cells were washed with PBS for twice (3-5 min/time) and blocked with 10% BSA for 10 min at room temperature. Then aspirating the blocking solution from the specimen without washing the specimen, FBs were incubated with Ki67 primary antibody (Abcam, Cambridge, UK; Rabbit polyclonal antibody) for 1 h. The FBs were washed with PBS and incubated with antirabbit (Invitrogen) away light for 40 min at room temperature. Lastly, Hoechst was used for staining for 10 min at room temperature. Images were obtained on Olympus (FV1000 and FV1200) confocal microscopes.

Cell cycle detection by flow cytometry
After digestion and centrifugation by trypsin, hFFs are washed with PBS for 3 times, 5 min/time and are resuspended in 75% alcohol precooled with 4°C and fixed at 4°C overnight. The next day, the fixed solution is centrifuged and discarded and hFFs are washed with PBS for 3 times, 5 min/time. After that, the PI/RNase mixture (4087S, CST) is added and incubated for 15~20 min in the dark. Lastly, the reaction solution is discarded after centrifugation and PBS is added into hFFs for detection by flow cytometry. The calculation formula of proliferation index (PI) is: PI = (S + G2/M)/ (G0/G1 + S + G2/M) × 100%.

RT-qPCR
After washing the cells twice with PBS, pre-cooled Trizol (1 ml/1 million cells) was added, lysed for 5 min, collected in 1.5 ml RNase-free EP tubes, centrifuged at 12,000 rpm for 5 min at 4°C, and the supernatant was added with chloroform of 1/5

GSEA
GSEA was performed based on our mRNA data using the GSEA software version4.2.1. Gene-lists were derived from the mRNA data by ranking them by their expression levels between two samples. Since then, some genes from c2.cp.kegg.v7.5.1.symbols.gmt was mapped to a list of genes arranged in advance to calculate ES score (enrichment), which nPerm 1000 was arranged to calculate significant sex, and determine 2 set of genes that false discovery rate value (FDR) < 0.25 was considered to have significant sex.

Differential expression analysis
The DESeq2 was used for differential expression analysis between two conditions/groups. Using Benjamini and Hochberg's approach for adjusting P values, the false discovery rate was controlled. Genes with P-value < 0.01 were designated as differentially expressed. FDR < 0.01 and fold change ≥2 were the thresholds to determine whether to identify genes with significant differential expression using DESeq2.

GO enrichment and KEGG pathway enrichment analysis
Based on Wallenius non-central hypergeometric distribution, GOseq R packages analyzed DEG GO enrichment analyses [58], which allowed gene length bias to be adjusted.
KEGG [59] can analyze high-level functions and utilities of cells through molecular datasets (http://www.genome.jp/kegg/). KOBAS [60] software was used to count the enrichment of DEGs in the KEGG pathways.

PPI
Import the DEGs into the STRING database: http://string-db.org/) to obtain the PPI relationship between these DEGs. Afterwards, the PPI relationships of these DEGs were imported into Cytoscape for visualization [61].

WGCNA
A total of 95139 transcripts were detected, the six duplicate samples were all 0 removed, the transcripts with the same gene name were merged (average), and then MAD was used for screening. The genes with the top 10,000 median deviation values were selected for one-step WGCNA network construction. In WGCNA, genes with similar expression patterns are classified into the same module, and the first principal component of the module is called the intrinsic gene of the module.
Here we convert the experimental group into a 0-1 matrix, and we could estimate the correlation between the module and the trait by the expression of the module's eigengenes. The correlation mode was set to Pearson correlation.

Statistical analysis
For statistical analysis, student's t-test was used to compare the differences between two groups, one-way ANOVA Tukey's multiple comparisons test to compare four groups, two-way ANOVA Tukey's multiple comparisons test for G2, S and G1 percentage. Quantitative data are presented as mean ± SD of at least three independent experiments. P-values < 0.05 were considered as statistically significant.

Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author (

AUTHOR CONTRIBUTIONS
Y.L., L.L., and Zhixiang Xu supervised and controlled the research concept, designed and interpreted the data. Ziran Xu controlled research conception, designed the experiments, performed the overall experiment, analysed the data, and wrote the manuscript. T.Z. and Y.W. revised the manuscript and explained relevant data. L.Z. and J.T. analyzed the data with R and WGCNA.