Cells adhered to the cotton fibers formed aggregates 48 hours after seeding (Fig. 1A and Fig. 1B). Cells kept their round natural shape in the 3-D environment. Cells seeded in a classical monolayer adhered to the plastic dish (Fig. 1C), they had fibroblast like phenotype in the 2-D environment.
3 − 1 Data Quality
We first run a quality check on the data, to verify gene expression counts and experimental bias between samples. Experimental bias could be due to different techniques applied to the samples, different equipment, different reagents, and different handling processes. Our analysis showed that samples have been treated and analyzed under similar conditions. The expression values box plot shows that conditions are comparable and could be analyzed further (Fig. 2A).
In a principal component analysis (Fig. 2B), MSCs treated with cotton scaffold formed a distinct cluster from the control cells. This shows that both conditions are distinct and could be separated based on gene expression values.
3 − 2 Differently expressed genes
We identified a total of 67 DEG, 31 of which are up-regulated while 36 genes are down-regulated significantly in cotton condition based on Log2 fold change (p < 0.05) and q = 0.8 (Fig. 3 and Fig. 4). Figure 3A is a dot plot showing gene expression in the cells adhering to the scaffold. The red dots show a significant difference compared to the control sample. The genes that are not significantly different from the control are represented.in black (Fig. 3A). The DEG are represented in the bar chart (Fig. 3B), in terms of their expression value in the cotton condition.
In the group of the up-regulated genes, we observed genes involved in cell adhesion: APP, ITGA3 and LAMA4 [28–30], a gene involved in cell proliferation, angiogenesis and cell differentiation, the fibroblast growth factor-2 (FGF2) [31]. Gene with anti-apoptotic roles (TXNDC5) and genes involved in cell cycle and growth (TSG101, ZNF207) were also increased.
On the other hand, the list of down-regulated genes included mainly the cell cycle and growth genes (USP16, NUSAP1, ECT2). Furthermore, genes involved in cell proliferation TRMT1 [32], MAP2K1 and MAP4K4 [33], cell differentiation (SUN1, ECT2, MAP2K1), neurogenesis (ECT2, SMARCE1), cell motility, migration and angiogenesis (SPAG9) [34], apoptosis RABEP1, SRPRB [35], MAP4K4 [33], were also significantly down regulated. The full list of DEG is provided in supplementary data (supplementary Table 1A for the up-regulated genes and supplementary Table 1B for the down-regulated genes).
3–3 Functional enrichment analysis
The biological classification of the identified genes was further analyzed by the gene enrichment analysis
3-3-1 Gene enrichment analysis: Up Regulated Genes
GO term enrichment analysis was used to investigate the biological process and molecular functions of the up-regulated genes. The bar charts in (Fig. 4A and Fig. 4B) show the biological process and molecular function, the number of genes associated with each category and the statistical significance.
Our data revealed that up-regulated genes are principally enriched in the regulation of gene expression, followed by the extracellular matrix organization and extracellular structure organization. Their molecular functions are mainly related to protein heterodimerization activity, fibroblast growth factor receptor and choline transmembrane transporter activity. The top 10 most significant biological processes and molecular functions are represented in Fig. 4A and 4B. The top 20 of these processes together with the cellular localization are shown in supplementary data (supplementary Fig. 1). The top 20 corresponding most significantly increased (p < 0.05) GO term molecular functions are shown in Table 1. The top 20 biological process for the up-regulated genes-are shown in supplementary data (supplementary Table 2).
Table 1
Representation of the top 20 Gene enrichment based on GO Molecular Function of up-regulated genes
Term | Overlap | p value | Adjusted p value | Odds.Ratio | Genes |
choline transmembrane transporter activity (GO:0015220) | 1/5 | 0.008 | 0.145 | 166.38 | SLC44A2 |
aspartic endopeptidase activity, intramembrane cleaving (GO:0042500) | 1/6 | 0.009 | 0.145 | 133.09 | SPPL2A |
protein phosphatase 2B binding (GO:0030346) | 1/6 | 0.009 | 0.145 | 133.09 | ATP2B4 |
nitric-oxide synthase binding (GO:0050998) | 1/7 | 0.011 | 0.145 | 110.91 | ATP2B4 |
P-type calcium transporter activity (GO:0005388) | 1/9 | 0.014 | 0.145 | 83.17 | ATP2B4 |
intramolecular transferase activity, phosphotransferases (GO:0016868) | 1/9 | 0.014 | 0.145 | 83.17 | PMM2 |
nuclear export signal receptor activity (GO:0005049) | 1/11 | 0.017 | 0.145 | 66.53 | CSE1L |
P-type ion transporter activity (GO:0015662) | 1/13 | 0.020 | 0.145 | 55.44 | ATP2B4 |
GTPase inhibitor activity (GO:0005095) | 1/13 | 0.020 | 0.145 | 55.44 | IPO5 |
aspartic-type endopeptidase activity (GO:0004190) | 1/15 | 0.023 | 0.145 | 47.51 | SPPL2A |
nuclear import signal receptor activity (GO:0061608) | 1/16 | 0.025 | 0.145 | 44.34 | IPO5 |
organic cation transmembrane transporter activity (GO:0015101) | 1/16 | 0.025 | 0.145 | 44.34 | SLC44A2 |
histone acetyltransferase binding (GO:0035035) | 1/21 | 0.032 | 0.159 | 33.25 | HIF1A |
protein heterodimerization activity (GO:0046982) | 2/188 | 0.034 | 0.159 | 7.34 | ITGA3; HIF1A |
fibroblast growth factor receptor binding (GO:0005104) | 1/23 | 0.035 | 0.159 | 30.22 | FGF2 |
nuclear localization sequence binding (GO:0008139) | 1/24 | 0.037 | 0.159 | 28.91 | IPO5 |
mRNA 5'-UTR binding (GO:0048027) | 1/25 | 0.038 | 0.159 | 27.70 | IGF2BP2 |
metal ion transmembrane transporter activity (GO:0046873) | 1/27 | 0.041 | 0.159 | 25.57 | ATP2B4 |
receptor antagonist activity (GO:0048019) | 1/28 | 0.043 | 0.159 | 24.62 | FST |
chemokine binding (GO:0019956) | 1/32 | 0.048 | 0.164 | 21.44 | FGF2 |
The table represents the results of the GO Molecular Function of up-regulated genes in cells with cotton. |
3-3-2 Up-regulated Genes: KEGG pathway enrichment analysis
Pathway analysis showed the up-regulated genes were enriched in pathways in cancer, followed by proteoglycans in cancer, regulation of actin cytoskeleton and PI3K-AKt signaling pathway. The up-regulated genes were also enriched in other pathways such as ECM-receptor interaction and focal adhesion (Fig. 4C). The list of the top 20 enriched pathways is provided in the supplementary data (supplementary Fig. 1D). Table 2 represents the top 20 gene enrichment terms based on “KEGG Pathways” showing the genes involved in each pathway and their significance.
Table 2
Representation of the top 20 Gene enrichment terms based on “KEGG Pathways” in up-regulated genes.
Term | Overlap | p value | Adjusted p value | Odds Ratio | Genes |
Fructose and mannose metabolism | 2/33 | 0.001 | 0.056 | 44.36 | PMM2;GMPPA |
Amino sugar and nucleotide sugar metabolism | 2/48 | 0.003 | 0.056 | 29.87 | PMM2;GMPPA |
Proteoglycans in cancer | 3/205 | 0.004 | 0.056 | 10.48 | RDX;HIF1A;FGF2 |
Regulation of actin cytoskeleton | 3/218 | 0.005 | 0.056 | 9.84 | ITGA3;RDX;FGF2 |
Mitophagy | 2/68 | 0.005 | 0.056 | 20.80 | CALCOCO2; HIF1A |
Central carbon metabolism in cancer | 2/70 | 0.005 | 0.056 | 20.18 | G6PD;HIF1A |
ECM-receptor interaction | 2/88 | 0.008 | 0.063 | 15.94 | ITGA3;LAMA4 |
Pathways in cancer | 4/531 | 0.009 | 0.063 | 5.47 | ITGA3;LAMA4;HIF1A;FGF2 |
Small cell lung cancer | 2/92 | 0.009 | 0.063 | 15.23 | ITGA3;LAMA4 |
Choline metabolism in cancer | 2/98 | 0.010 | 0.063 | 14.28 | SLC44A2;HIF1A |
Amoebiasis | 2/102 | 0.011 | 0.063 | 13.70 | RAB5B;LAMA4 |
PI3K-Akt signaling pathway | 3/354 | 0.017 | 0.091 | 5.99 | ITGA3;LAMA4; FGF2 |
Ribosome | 2/158 | 0.025 | 0.122 | 8.76 | RPS10-NUDT3;MRPL34 |
Kaposi sarcoma-associated herpesvirus infection | 2/193 | 0.036 | 0.164 | 7.14 | HIF1A;FGF2 |
Focal adhesion | 2/201 | 0.039 | 0.165 | 6.85 | ITGA3;LAMA4 |
Pentose phosphate pathway | 1/30 | 0.046 | 0.176 | 22.92 | G6PD |
Ras signaling pathway | 2/232 | 0.050 | 0.176 | 5.92 | RAB5B;FGF2 |
Calcium signaling pathway | 2/240 | 0.053 | 0.176 | 5.72 | ATP2B4;FGF2 |
African trypanosomiasis | 1/37 | 0.056 | 0.176 | 18.46 | LAMA4 |
Salmonella infection | 2/249 | 0.057 | 0.176 | 5.51 | RAB5B;CSE1L |
The table shows gene enrichment terms based on “KEGG Pathways” in up-regulated genes in cells with cotton. |
3-3-3 Protein-protein interaction network construct for the up-regulated DEG
The protein-protein interaction network provided the molecular organization of the up-regulated DEG. The up-regulated DEG are connected to a number of Hub genes, involved in cell proliferation and differentiation (TP53, MDM2, MAPK, beta-tubulin, GABARAP [36–41], genes involved in anti-aging, cell protection and anti-inflammatory roles (SUMO2, HGS, Grin1, HSPAS) [42–46]. DEG genes were also connected to genes involved in cancer (SHC1, MEPCE) [47, 48] (Fig. 4D). The top 10 Gene enrichment terms based on protein-protein interaction Hub Genes are represented in supplementary data (supplementary Table 3) showing the genes connected to each hub gene and the significance.
3-3-4 Go enrichment analysis: Down Regulated Genes
GO enrichment analysis of the down-regulated genes was also categorized into three functional groups (biological processes, molecular functions and cellular localization). The results showed that down-regulated genes are enriched in various biological processes mainly in the chromatin disassembly, nucleosome disassembly, protein-DNA complex disassembly, regulation of stress-activated MAPK cascade, positive regulation of transcription, DNA-template (Fig. 5A). Their molecular functions are mainly associated with MAP-Kinase scaffold activity, nucleosomal DNA binding, nuclear receptor coactivator activator, ubiquitin-binding, creatine kinase activity, procollagen-proline 4 dioxygenase activity, JUNK kinase binding, procollagen proline dioxygenase activity, nuclear receptor binding, aldehyde dehydrogenase (NAD+) activity (Fig. 5B). The bar chart of the full list of the top 20 biological processes, molecular functions and cellular localization for the down-regulated DEG is provided in supplementary data (supplementary Fig. 2). Table 3. shows the top 20 gene enrichment terms based on GO molecular function of the down-regulated DEG. The top 20 biological processes in the down-regulated DEG are shown in supplementary data (supplementary Table 4).
Table 3
Representation of top 20 Gene enrichment terms based on “GO Molecular Function” in Down-regulated genes
Term | Overlap | p value | Adjusted p value | Odds Ratio | Genes |
Tat Protein Binding (GO:0030957) | 2/9 | 0.000 | 0.006 | 167.71 | MDFIC;SMARCA4 |
MAP-kinase Scaffold Activity (GO:0005078) | 2/10 | 0.000 | 0.006 | 146.74 | SPAG9;MAP2K1 |
Nucleosome Binding (GO:0031491) | 3/62 | 0.000 | 0.006 | 30.67 | SMARCE1;RNF4;SMARCA4 |
N-methyltransferase Activity (GO:0008170) | 2/27 | 0.001 | 0.024 | 46.92 | TRMT1;METTL9 |
Nucleosomal DNA Binding (GO:0031492) | 2/34 | 0.002 | 0.030 | 36.64 | SMARCE1;SMARCA4 |
Nuclear Receptor Coactivator Activity (GO:0030374) | 2/50 | 0.004 | 0.054 | 24.41 | RNF4;DCAF6 |
RNA Polymerase II-specific DNA-binding Transcription Factor Binding (GO:0061629) | 3/228 | 0.008 | 0.080 | 7.98 | SMARCE1;MDFIC;SMARCA4 |
Histone Deubiquitinase Activity (GO:0140934) | 1/5 | 0.010 | 0.080 | 142.57 | USP16 |
Ubiquitin Binding (GO:0043130) | 2/80 | 0.010 | 0.080 | 15.00 | TSG101;USP16 |
Chromatin DNA Binding (GO:0031490) | 2/83 | 0.010 | 0.080 | 14.44 | SMARCE1;SMARCA4 |
Oxysterol Binding (GO:0008142) | 1/6 | 0.011 | 0.080 | 114.05 | OSBPL5 |
Procollagen-Proline Dioxygenase Activity (GO:0019798) | 1/7 | 0.013 | 0.080 | 95.04 | P4HA1 |
Peptidyl-Proline 4-Dioxygenase Activity (GO:0031545) | 1/7 | 0.013 | 0.080 | 95.04 | P4HA1 |
RNA Polymerase I Core Promoter Sequence-Specific DNA Binding (GO:0001164) | 1/8 | 0.014 | 0.080 | 81.46 | SMARCA4 |
Carbonyl Reductase (NADPH) Activity (GO:0004090) | 1/8 | 0.014 | 0.080 | 81.46 | DHRS7 |
NADP-retinol Dehydrogenase Activity (GO:0052650) | 1/9 | 0.016 | 0.080 | 71.27 | DHRS7 |
RNA Polymerase I Transcription Regulatory Region Sequence-Specific DNA Binding (GO:0001163) | 1/9 | 0.016 | 0.080 | 71.27 | SMARCA4 |
JUN Kinase Binding (GO:0008432) | 1/9 | 0.016 | 0.080 | 71.27 | SPAG9 |
Nuclear Receptor Binding (GO:0016922) | 2/115 | 0.019 | 0.085 | 10.33 | SMARCE1;SMARCA4 |
tRNA (Guanine) Methyltransferase Activity (GO:0016423) | 1/11 | 0.0196 | 0.087 | 57.01 | TRMT1 |
The table shows Gene enrichment terms based on “GO Molecular Function” in Down-regulated genes in cells with cotton. |
KEGG pathway analysis showed the down-regulated genes were mainly enriched in metabolic pathways, arginine and proline metabolism, protein export, tryptophan metabolism and pyruvate metabolism (Fig. 5C). The bar chart for the top 20 pathway analysis of the down-regulated DEG is provided in supplementary data (Fig. 2D). Table 4 represents the top 20 KEGG pathways in the down-regulated DEG.
Table 4
Representation of top 20 Gene enrichment terms based on “KEGG Pathways” in the Down-regulated genes.
Term | Overlap | p value | Adjusted p value | Odds Ratio | Genes |
Hepatocellular carcinoma | 3/168 | 0.003 | 0.216 | 10.91 | SMARCE1;MAP2K1;SMARCA4 |
Arginine and proline metabolism | 2/50 | 0.004 | 0.216 | 24.41 | P4HA1;ALDH7A1 |
Glycerolipid metabolism | 2/61 | 0.005 | 0.216 | 19.85 | AGPAT1;ALDH7A1 |
Phospholipase D signaling pathway | 2/148 | 0.029 | 0.316 | 7.98 | MAP2K1;AGPAT1 |
Histidine metabolism | 1/22 | 0.039 | 0.316 | 27.13 | ALDH7A1 |
Protein export | 1/23 | 0.041 | 0.316 | 25.90 | SRPRB |
Collecting duct acid secretion | 1/27 | 0.048 | 0.316 | 21.91 | ATP6V0A1 |
Ascorbate and aldarate metabolism | 1/30 | 0.053 | 0.316 | 19.64 | ALDH7A1 |
beta-Alanine metabolism | 1/30 | 0.053 | 0.316 | 19.64 | ALDH7A1 |
Citrate cycle (TCA cycle) | 1/30 | 0.053 | 0.316 | 19.64 | ACLY |
SNARE interactions in vesicular transport | 1/33 | 0.058 | 0.316 | 17.80 | STX5 |
Thyroid cancer | 1/37 | 0.065 | 0.316 | 15.82 | MAP2K1 |
Thermogenesis | 2/232 | 0.065 | 0.316 | 5.05 | SMARCE1;SMARCA4 |
Glycine, serine and threonine metabolism | 1/40 | 0.070 | 0.316 | 14.60 | ALDH7A1 |
Bladder cancer | 1/41 | 0.071 | 0.316 | 14.23 | MAP2K1 |
Tryptophan metabolism | 1/42 | 0.073 | 0.316 | 13.88 | ALDH7A1 |
Fat digestion and absorption | 1/43 | 0.075 | 0.316 | 13.55 | AGPAT1 |
Fatty acid degradation | 1/43 | 0.075 | 0.316 | 13.55 | ALDH7A1 |
Endocytosis | 2/252 | 0.075 | 0.316 | 4.64 | RABEP1;TSG101 |
Pyruvate metabolism | 1/47 | 0.081 | 0.316 | 12.37 | ALDH7A1 |
The shows gene enrichment terms based on “KEGG Pathways” in the Down-regulated genes in cells with cotton |
3-3-5 Protein-protein interaction network construct for the down-regulated DEG
The down-regulated DEG were connected to several Hub genes involved in cell proliferation, cell differentiation, inflammatory response (NR3C1 MAPK8, NCK1) [49–51], cell survival and apoptosis (SP1, HSPAS) [42, 52], cancer cell proliferation (AR, CDK1, SHC1) [47, 53, 54] and stem cell proliferation and phenotype (cMYC) [55]. Similarly, to the up-regulated DEG, the down-regulated were also interacting with HGS (Hepatocyte growth factor) (Fig. 5D). Table 5 in supplementary data shows the top 10 Hub genes with the up-regulated DEG connections showing the genes connected to each hub gene and the significance.
3–4 Analysis with a more stringent threshold
The threshold was increased to 0.9 to ensure maximum significant differences in the DEG. The MD plot shows the up and down-regulated genes in the cotton-seeded cells based on the Log2 fold change. Stringent threshold revealed 18 up-regulated genes and 15 down-regulated genes (Fig. 6B). The list of the whole genes is provided in supplementary data supplementary Table 6 Venn diagram (Fig. 7), shows the common DEG between the default threshold (0.8) and the stringent threshold (0.9).
Genes involved in cell adhesion such as APP, ITGA3 and LAMA4 were significantly up-regulated in the MSCs seeded on the cotton fibers. FGF2 as well as the genes with anti-apoptotic role (TXNDC5) and cell cycle and growth role (TSG101) were also significantly up-regulated.
On the other hand, USP16 and NUSAP1 genes involved in cell cycle and growth were down-regulated. Genes involved in cell proliferation, cell motility, cell migration angiogenesis and apoptosis MAP4K4 [33] and gene involved in cell differentiation (SUN1) were significantly down-regulated.