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

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.

Figure 1. Soft substrates promoted aspect ratio of hFFs. (A) Phase-contrast images of hFFs on TCP and soft substrates with or without TGF-β1 at day 1. Scale bars: 200 μm (up) or 100 μm (down). (B) Aspect ratio of hFFs on TCP and soft substrates with or without TGF-β1 at day 1. ^***P < 0.001 (mean, n = 10).

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 |log₂foldchange|≥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 down-regulated 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).

Figure 2. Differential expression analysis of soft/TCP. (A) Venn diagram of soft and TCP. (B) Scatter plot illustrated DEGs of soft/TCP. (C) Heat-map cluster analysis of DEGs. (D) KEGG pathways of up and down regulated DEGs. The color and horizontal axe were −log₁₀ (p value) and gene number respectively. (E) GO analysis of up- and down- regulated DEGs. The GO analysis categorized mRNA into different groups: BP, CC, and MF. The color and horizontal axe were −log₁₀ (p value) and gene number respectively.

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 down-regulated 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.

Figure 3. Soft substrates inhibited the proliferation of hFFs. (A) GSEA for cell cycle of mRNAs in soft/TCP. Enrichment score of cell cycle genes against their expression profile of soft/TCP. The x axis was the level, and the y axis was the enrichment score of these genes (left). Heat map showed the distribution of mRNAs expression of soft/TCP enriched under the cell cycle KEGG entry in all samples under the gene set, with each row representing one sample. Each line represented a gene, the color from blue to red represented mRNA expression from low to high (right). (B) Immunostaining against Ki67 of soft and TCP with or without TGF-β1 on day 1. Scale bars: 50 μm. (C) The positive rate of Ki67 in B. ^***P < 0.001 (mean, n = 3 randomly selected fields from triplicate samples). (D) Detection of cell cycle by flow cytometry for soft and TCP with or without TGF-β1. (E) The percentage of PI. ^**P < 0.01; ^***P < 0.001 (mean, n = 3). (F) The percentage of G1-, S-, and G2-phase in the cell cycle. ^*P < 0.05 vs. TCP; ^&&P < 0.01, ^&&&P < 0.001 vs. TCP+ TGF-β1 (mean, n = 3).

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 down-regulate 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.

Figure 4. Soft substrates inhibited contractility of hFFs. (A) GSEA for vascular smooth muscle contraction of mRNAs in soft/TCP. Enrichment score of vascular smooth muscle contraction genes against their expression profile of soft/TCP. The x axis was the level, and the y axis was the enrichment score of these genes (left). Heat map showed the distribution of mRNAs expression of soft/TCP enriched under the vascular smooth muscle contraction KEGG entry in all samples under the gene set, with each row representing one sample. Each line represented a gene, the color from blue to red represented mRNA expression from low to high (right). (B) Quantitative reverse transcription PCR (qRT-PCR) analysis comparing α-SMA expression levels in soft and TCP. ^***P < 0.001 (mean, n = 3). (C) Immunostaining against PERIOSTIN and the positive rate on day 1. Scale bars: 50 μm. ^***P < 0.001 (mean, n = 3 randomly selected fields from triplicate samples).

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 log₂foldchange 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.

Figure 5. PPI network analysis of DEGs. (A) The PPI network of 698 DEGs. (B) Cluster1 and cluster 2 sub-networks obtained after MCODE analysis of PPI network. (C) MCC sub-network obtained after CytoHubba analysis of PPI network. (D) 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 log₂ (fold change). The higher and lower the expression is, the redder and bluer it is, respectively.

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

Figure 6. WGCNA and significant module recognition. (A) Sample clustering found no obvious outliers. (B) Analysis of network topology for various soft-thresholding powers. The left panel shows the scale-free fit index (y-axis) as a function of the soft-thresholding power (x-axis). Power 14 was chosen because the fit index curve flattened out upon reaching a high value (> 0.85). The right panel displays the mean connectivity (degree, y-axis) as a function of the soft-thresholding power (x-axis). (C) Clustering dendrogram of all mRNAs dataset based on topological overlap. Each module is given a unique colour and represents a cluster of coexpressed genes. (D) The eigengene adjacency heatmap was used to illustrate the relationship between eigengenes and phenotypic traits. (E) Identification of key modules related to the impact of soft matrices on hFFs. Heatmap displaying the correlations and significant differences between gene modules and samples. Each row corresponds to a module eigengene and each column to a trait. Correlation coefficients and P values are displayed in rectangles. Color-coded by relevance based on a color legend. Blue rectangles represent negative correlations between modules and samples, and red rectangles represent positive correlations between modules and samples.

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 log₂foldchange 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.

Figure 7. Screening and analysis of hub genes. (A) Veen diagram of CytoHubba (MCC), MCODE (cluster 1 and 2) and WGCNA (blue and turquoise modules). (B) The PPI network of 63 hub DEGs. (C) KEGG analysis of 63 hub DEGs. (D) GO analysis of 63 hub DEGs. (E) The PPI network of 12 hub DEGs in cell cycle. (F) The PPI network of 25 hub DEGs in cell division. (G) The legend of networks. The round represents down-regulated DEGs, the squar 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 log₂ (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).

Cell culture

Normal human FBs were isolated from foreskin tissues of healthy donors (provided by Dr. Lin Chen) [57]. Briefly, these tissues were washed with PBS, cut into 1 mm³ pieces. After cells were digested with 0.15% collagenase (Roche Applied Science, USA), the suspension was collected, filtered, centrifuged, resuspended, and then cultured in Dulbecco’s modified Eagle’s medium (Gibco, USA) containing 10% fetal bovine serum (Gibco, USA), 1% NEAA (Gibco, USA), 1% PS, 0.5% GlutaMAX (Gibco, USA) with or without 10 ng/ml TGF-β1 (PeproTech, Inc. USA), at 37°C in a humidified incubator containing 5% CO₂. The centrifuge tubes were obtained from Jet Biofil Co., Ltd, China. The cells of FBs at passages 3 to 5 were used to further experiment.

Preparation of the collagen I

Firstly under the condition of ice bath, the amount of solution as described in the following table separately added to the 1.5 mL EP tube and vortexed for 3 min. Add it dropwise in an amount of 300 μl to one well of a 24-well plate (NEST Biotechnology Co., Ltd, China), rotate it slowly for a circle to make it evenly and flatly attached to the bottom of the well plate, incubate at 37°C for 30 min, lay 200 μl of 0.1% gelatin in each well, and incubate at 37°C for 1 h, it could then be used to culture cells.

RNA-seq

RNA of TCP and soft concentration as well as purity was measured using NanoDrop 2000 (Thermo Fisher Scientific, Wilmington, DE, USA). The RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) was used for assessing RNA integrity.

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 anti-rabbit (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 Trizol volume, vortex for 1 min, centrifuge at 12000 rpm for 15 min at 4°C, add an equal volume of isopropanol to the supernatant, invert and mix, let stand for 10 min, centrifuge at 12000 rpm for 30 min at 4°C, wash the precipitate with 75% ethanol, centrifuge at 7500 rpm for 5 min at 4°C, and dissolve the pellet with DEPC water after drying. cDNA was obtained according to the instructions of Hifair^® II 1^st Strand cDNA Synthesis SuperMix for qPCR (gDNA digester plus) kit, and qPCR was performed according to Hieff^® qPCR SYBR Green Master Mix (High Rox Plus) kit (α-SMA forward: CCCATCTATGAGGGCTATGC, reverse: CACGCTCAGTCAGGATCTTC).

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 (Lisha Li: 43973966@qq.com; lilisha@jlu.edu.cn on reasonable request. Correspondence and requests for materials should be addressed to Ziran Xu (1240300750@qq.com; xuzr18@mails.jlu.edu.cn) and Lisha Li (43973966@qq.com; lilisha@jlu.edu.cn).