Identification of critical genes and pathways involved in intervertebral disc degeneration (IDD) An integrated bioinformatics analysis

doi:10.21203/rs.3.rs-32656/v1

Download PDF

Research article

Identification of critical genes and pathways involved in intervertebral disc degeneration (IDD) An integrated bioinformatics analysis

https://doi.org/10.21203/rs.3.rs-32656/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: In western countries and China, back and neck pain has become a common problem that bothers daily life and severely influences the quality of our daily life. Among all factors that lead to chronic neck and back pain, IDD is the one that couldn’t be easily neglected.

Methods: This study aims to figure out the critical genes and pathways involved in the development of IDD and provide a new aspect of following investigations on the etiology of IDD. We firstly systemically searched the GEO database and identified the differentially expressed genes (DEGs) from the expression profile dataset we selected. We secondly constructed the protein-protein interaction (PPI) network for DEGs, identified the top ten hub genes from the whole PPI network and found two statically and medically significant modules from the network, we then performed the GO and KEGG analysis on the DEGs, top ten hub genes, the PPI network and the two statically and medically modules. In the end, we provided the primers of the mRNAs of all DEGs, which will be useful for the validation experiment of this study.

Results: FN1, MMP2, POSTN, COL3A1, TIMP3, FBN1, GJA1, TGFBI, EFEMP1 and ID1 were top ten hub genes identified from this study, and they may play a vital role in the development of IDD. Angiogenesis and integrin binging are crucial biological process and molecular function defined in this study, which are worthy of being intensely investigated.

Conclusion: More studies on the top ten hub genes, the role of angiogenesis and integrin binding in IDD are urgently needed, which will benefit the prevention, screening, diagnosis and prognosis of IDD.

Epigenetics & Genomics

IDD is a major cause of the neck or back pain in many countries, including the United States and China, which will bring mental and physical suffering to the patient himself and a heavy financial burden to the patient’s family and the whole society (1). And but most of the studies in this field is about the management and diagnosis IDD. Risbud et al. reported in 2014 that intervertebral disc degeneration is related to cytokines such as TNF, IL-1α, IL-1β, etc. (2), But the underlying molecular mechanism, especially genetic mechanism is seldom covered.

In this study, we aimed to reveal the genetic factors and pathways that closely related to the etiology of IDD by the integrated bioinformatics analysis and provide a new idea for following investigations on the etiology of IDD.

We firstly systemically searched the GEO database (3-5) and identified profile dataset GSE34095 for expression data extraction and analysis. We assigned all the samples into two groups, one is the normal group, which contains the tissue obtained from healthy people; the other one is the IDD group, which includes the tissue obtained from patients diagnosed as IDD. After getting the expression level of all genes, we set the threshold value to screen the DEGs for the following. In this study, 61 DEGs were identified, of which 43 DEGs were upregulated, 18 DEGs were downregulated. We then constructed the PPI network (6,7) of all DEGs (Figure1). We then used the Cytoscape (version 3.7.2) (8,9) to modify the network and found two medically and statically significant modules by the applet of the Cytoscape (version 3.7.2) called MCODE (version 1.6). In the end, we applied the cytohubba (version 0.0.1) (10), an applet of the Cytoscape, to detect the top ten hub genes of the whole PPI network. They are FN1, MMP2, POSTN, COL3A1, TIMP3, FBN1, GJA1, TGFBI, EFEMP1 and ID1.

Thirdly we conducted the GO (11) and KEGG (12) pathway enrichment analysis on the DAVID (13), which enables us to investigate the relationship between the genes and biological process, cell component, molecular function and pathways of our bodies. The results are as follows: for all DEGs identified, they are mainly enriched in angiogenesis in the term of biological process, extracellular region in the term of cell component, integrin binding in the term of molecular function and PI3K-Akt signaling pathway in the term of pathway. For genes involved in module 1, they are primarily enriched in angiogenesis in the term of biological process, extracellular matrix in the term of cell component, integrin binding in the term of molecular function and proteoglycans in cancer in the term of pathway. For the genes involved in module 2, they are primarily enriched in ribosomal protein import into nucleus in the term of biological process, nuclear membrane in the term of cell component, poly(A) RNA binding in the term of molecular function and nothing could be detected for the module 2 in the term of pathway. For all top ten hub genes, they are mainly enriched in extracellular matrix organization in the term of biological process, extracellular matrix in the term of cell component, integrin binding in the term of molecular function and proteoglycans in cancer in the term of pathway.

For the sake of precision and credibility of our functional enrichment analysis and the DEGs we identified, we ought to perform an independent validation experiment. But due to the outbreak of the COVID-19 virus, we are unable to get enough biopsy samples for validation. So, we provided the primers of the mRNAs of all DEGs in this study, which could be useful for furthering validation if someone interested in this. All the primers are shown in Table S2.

Microarray platform

The dataset GSE34095 is based on the GPL96 ([HG-U133A] Affymetrix Human Genome U133A Array)

Expression analysis of DEGs

We used the GEO2R function on the website of the GEO database to sort the data into two groups, one is the normal group, which contains the tissue obtained from healthy people; the other one is the IDD group, which includes the tissue obtained from patients diagnosed as IDD. And analyzed dataset with all the relevant sets default and exported to Excel. For all four profile datasets, a gene with |log2‑fold change|>0.5 and P-value<0.05 were considered as a DEG.

GO and KEGG pathway enrichment analysis.

The Database for Annotation Visualization and Integrated Discovery (DAVID) is a website that provides a wide range of functional annotation tools to investigate the effect of genes in the term of biological process, molecular function, cell components, etc. Identified DEGs were investigated further using DAVID (version 6.8), GO, and KEGG pathway enrichment analyses. While analyzing all genes in the whole network, genes involved in the two modules identified from the PPI network, P<0.05 and gene counts of >10 indicates a statistically significant difference in the functional enrichment analysis; For the top 10 hub genes, P<0.05 and gene counts of >5 indicates a statistically significant difference in the functional enrichment analysis.

Integration of the PPI network.

The protein-protein interaction network (PPI) of all proteins expressed by the DEGs was mapped by the STRING (https://string-db.org/; version 11.0) to evaluate the interactions between the DEGs. The Cytoscape applet Molecular Complex Detection (MCODE; version 1.6) was used to detect the medically and statistically significant regions within the whole PPI network, which is called modules. Interactions with low confidence (a combined score >0.150) were defined as statistically significant. Cytoscape software (version 3.7.2) was used to visualize and adjust PPI networks. And based on the degree levels in the Cytoscape plugin cytoHubba (version 0.1), the top 10 ranked genes were defined as hub genes.

Design of primers

We designed primers for the mRNAs of all DEGs. (Table S1) We firstly searched the gene on the nucleotide tool of the NCBI website (14-17). And used the “pick primer” function with the minimum PCR product size is 50, maximum PCR product size is 200, minimum primer melting temperature is 58.0, maximum primer melting temperature is 62.0, maximum primer melting temperature difference is 2.0 and all other settings defaulted. The primer with minimum length, minimum melting temperature difference between forward and reverse primer, minimum GC ratio difference between forward and reverse primer, the GC ratio closest to 45% and minimum self-complementarity will be selected.

The identification of DEGs

In this study, 61 DEGs were identified, of which 43 were upregulated, and 18 were downregulated. (Table S2) Volcano plots (41,42) for the profile dataset GSE34095 was constructed. (Figure 1)The red lines represent the threshold values of the identification of DEGs. And an expression heatmap is drawn for all DEGs. (Figure 2)

Functional Enrichment analysis

To identify the pathways which had the most biologically and statistically significant involvement with the genes identified, we used the DAVID for GO and KEGG pathway analysis. For all DEGs identified, they are mainly enriched in angiogenesis in the term of biological process, extracellular region in the term of cell component, integrin binding in the term of molecular function and PI3K-Akt signaling pathway in the term of pathway, and PI3K-Akt signaling pathway is the only pathway shown in the result. (Figure 3, a-c) For genes involved in module 1, they are primarily enriched in angiogenesis in the term of biological process, extracellular matrix in the term of cell component, integrin binding in the term of molecular function and proteoglycans in cancer in the term of pathway. (Figure 4, a-d) For the genes involved in module 2, they are primarily enriched in ribosomal protein import into nucleus in the term of biological process, nuclear membrane in the term of cell component, poly(A) RNA binding in the term of molecular function and nothing could be detected for the module 2 in the term of pathway. (Figure 5, a-c) For all top ten hub genes, they are mainly enriched in extracellular matrix organization in the term of biological process, extracellular matrix in the term of cell component, integrin binding in the term of molecular function and proteoglycans in cancer in the term of pathway. (Figure 6, a-d)

PPI network construction and module analysis

Interactions between the identified DEGs were elucidated by the PPI network. (Figure 7 a) In total, there were 56 nodes in the network, the average number of neighbors is 10.179, the clustering coefficient is 0.484, and the network heterogenicity is 0.661. The top ten hub genes were FN1, MMP2, POSTN, COL3A1, TIMP3, FBN1, GJA1, TGFBI, EFEMP1 and ID1. And after analyzing with The Cytoscape plugin Molecular Complex Detection (MCODE; version 1.6) of Cytoscape (3.7.2), we got two statically and medically significant modules. (Figure 7 b,c)

FN1 gene, which is also called fibronectin 1, is a gene with at least 20 transcriptional variants. It mainly expressed in the placenta and liver compared to other organs. Many studies have revealed that the difference in the expression of FN1 will lead to gastric cancer (18), prostate cancer (19) and chemo-resistant breast cancer (20). Komori T in 2019, reported the FN1 is also involved in the function of osteoblast (21).

MMP2 is a gene locates on chromosome 16. Some studies have pointed out that this gene can be modulated by the FasL in the osteoblast (22). And Tang H et al. reported in 2019 that the MMP2 derived from mature osteoblasts is involved in the angiogenesis of epithelial cells (23), which is corresponding to the results of GO enrichment analysis.

POSTN encodes a secreted extracellular matrix protein that responsible for tissue development and regeneration, including wound healing and ventricular remodeling following myocardial infarction. Some studies have illustrated that this gene is involved in the carcinogenesis of ovarian cancer (24), colorectal cancer (25) and breast cancer (26), but no available evidence shows the relationship between this gene and the IDD.

COL3A1 is a gene that relatively poorly studied compared to other hub genes. Xu W et al. reported in 2018 that COL3A1 may be involved in the development of osteosarcoma (27).

Han XG et al. reported in 2018 that the TIMP3 would improve the sensitivity of osteosarcoma to cisplatin by reducing the production of IL-6 (28).

FBN1, GJA1, TGFBI, EFEMP1 and ID1 are mostly involved in the carcinogenesis of many types of cancers or other genetic disorders, but almost no studies covered the effects of these five genes on the bone or IDD (29-42).

Judging from the results of GO and KEGG pathway enrichment analysis, we can easily find that angiogenesis and integrin binging are shown in the results of at least two sets of genes in the term of respectively biological process and molecular function. But no study can be searched in PubMed that discusses the relationship between angiogenesis or integrin binding and the formation of bone or IDD. And relevant studies are urgently needed.

However, there are some inevitable limitations in this study. First, the number of available profile datasets in the database was not sufficient to analyze more DEGs, which may cause the imprecision of the results. So, more microarray data that compare the gene expression level between the samples extracted from healthy people and samples obtained from patients diagnosed as IDD is urgently needed. Secondly, The microarray data available now is mostly about people in Taiwan, which means, on the one hand, the expression level difference caused by the difference of race is unable to be identified; on the other hand, the heterogenicity between studies was dramatically huge, which indirectly decreases the credibility of studies included and results.

Due to the limitation of the condition and the outbreak of the COVID-19 virus, we are unable to collect enough biopsy samples and conduct qRT-PCR for validation with our data. We have to use the data from others for validation. We additionally provided the primers of mRNAs of all DEGs, which could be helpful for the furthering validation experiments.

More studies on the top ten hub genes, the role of angiogenesis and integrin binding in IDD are urgently needed, which will benefit the prevention, screening, diagnosis and prognosis of IDD.

PVALUE: -log10(p-value); IDD: intervertebral disc degeneration; DEGs: differentially expressed genes; PPI network: protein-protein interaction network; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Gene and Genomes

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Availability of data and materials

The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.

Competing interests

The authors declare that they have no competing interests.

Funding

This study was supported by The Government of the Xuzhou City. (2018TD011).

Authors’ contributions

JZ and WZ designed the study, JZ and FX collected the data, carried out the data analysis, JZ produced the initial draft of the manuscript, CX contributed to drafting the manuscript, CX and JZ polished all the figures in the manuscript.

Acknowledgments

We would like to acknowledge the helpful comments on this paper received from our reviewers.

Navone SE, Marfia G, Giannoni A, et al. Inflammatory mediators and signalling pathways controlling intervertebral disc degeneration. Histol Histopathol. 2017;32(6):523‐542. doi:10.14670/HH-11-846
Risbud MV, Shapiro IM. Role of cytokines in intervertebral disc degeneration: pain and disc content. Nat Rev Rheumatol. 2014;10(1):44‐56. doi:10.1038/nrrheum.2013.160
Barrett T, Wilhite SE, Ledoux P, et al. NCBI GEO: archive for functional genomics data sets--update. Nucleic Acids Res. 2013;41(Database issue):D991‐D995. doi:10.1093/nar/gks1193
Clough E, Barrett T. The Gene Expression Omnibus Database. Methods Mol Biol. 2016;1418:93‐110. doi:10.1007/978-1-4939-3578-9_5
NCBI Resource Coordinators. Database resources of the National Center for Biotechnology Information. Nucleic Acids Res. 2016;44(D1):D7‐D19. doi:10.1093/nar/gkv1290
Szklarczyk D, Gable AL, Lyon D, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019;47(D1):D607‐D613. doi:10.1093/nar/gky1131
von Mering C, Huynen M, Jaeggi D, Schmidt S, Bork P, Snel B. STRING: a database of predicted functional associations between proteins. Nucleic Acids Res. 2003;31(1):258‐261. doi:10.1093/nar/gkg034
Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498‐2504. doi:10.1101/gr.1239303
Su G, Morris JH, Demchak B, Bader GD. Biological network exploration with Cytoscape 3. Curr Protoc Bioinformatics. 2014;47:8.13.1‐8.13.24. Published 2014 Sep 8. doi:10.1002/0471250953.bi0813s47
Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;8 Suppl 4(Suppl 4):S11. doi:10.1186/1752-0509-8-S4-S11
Gene Ontology Consortium: The Gene Ontology (GO) project in 2006. Nucleic Acids Res 34 (Database Issue): D322‑D326, 2006.
Ogata H, Goto S, Sato K, Fujibuchi W, Bono H and Kanehisa M: KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res 27: 29‑34, 1999.
Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44‐57. doi:10.1038/nprot.2008.211
Ye J, Coulouris G, Zaretskaya I, Cutcutache I, Rozen S, Madden TL. Primer-BLAST: a tool to design target-specific primers for polymerase chain reaction. BMC Bioinformatics. 2012;13:134. Published 2012 Jun 18. doi:10.1186/1471-2105-13-134
Li K, Brownley A. Primer design for RT-PCR. Methods Mol Biol. 2010;630:271‐299. doi:10.1007/978-1-60761-629-0_18
Rodríguez A, Rodríguez M, Córdoba JJ, Andrade MJ. Design of primers and probes for quantitative real-time PCR methods. Methods Mol Biol. 2015;1275:31‐56. doi:10.1007/978-1-4939-2365-6_3
NCBI Resource Coordinators. Database resources of the National Center for Biotechnology Information. Nucleic Acids Res. 2018;46(D1):D8‐D13. doi:10.1093/nar/gkx1095
Jiang K, Liu H, Xie D, Xiao Q. Differentially expressed genes ASPN, COL1A1, FN1, VCAN and MUC5AC are potential prognostic biomarkers for gastric cancer. Oncol Lett. 2019;17(3):3191‐3202. doi:10.3892/ol.2019.9952
Das DK, Ogunwobi OO. A novel microRNA-1207-3p/FNDC1/FN1/AR regulatory pathway in prostate cancer. RNA Dis. 2017;4(1):e1503.
Yang X, Hu Q, Hu LX, et al. miR-200b regulates epithelial-mesenchymal transition of chemo-resistant breast cancer cells by targeting FN1. Discov Med. 2017;24(131):75‐85.
Komori T. Regulation of Proliferation, Differentiation and Functions of Osteoblasts by Runx2. Int J Mol Sci. 2019;20(7):1694. Published 2019 Apr 4. doi:10.3390/ijms20071694
Svandova E, Vesela B, Lesot H, Sadoine J, Poliard A, Matalova E. FasL Modulates Expression of Mmp2 in Osteoblasts. Front Physiol. 2018;9:1314. Published 2018 Sep 19. doi:10.3389/fphys.2018.01314
Tang H, He Y, Li L, et al. Exosomal MMP2 derived from mature osteoblasts promotes angiogenesis of endothelial cells via VEGF/Erk1/2 signaling pathway. Exp Cell Res. 2019;383(2):111541. doi:10.1016/j.yexcr.2019.111541
Karlan BY, Dering J, Walsh C, et al. POSTN/TGFBI-associated stromal signature predicts poor prognosis in serous epithelial ovarian cancer. Gynecol Oncol. 2014;132(2):334‐342. doi:10.1016/j.ygyno.2013.12.021
Oh HJ, Bae JM, Wen XY, Cho NY, Kim JH, Kang GH. Overexpression of POSTN in Tumor Stroma Is a Poor Prognostic Indicator of Colorectal Cancer. J Pathol Transl Med. 2017;51(3):306‐313. doi:10.4132/jptm.2017.01.19
Ratajczak-Wielgomas K, Grzegrzolka J, Piotrowska A, et al. Expression of periostin in breast cancer cells. Int J Oncol. 2017;51(4):1300‐1310. doi:10.3892/ijo.2017.4109
Xu W, Li Z, Zhu X, Xu R, Xu Y. miR-29 Family Inhibits Resistance to Methotrexate and Promotes Cell Apoptosis by Targeting COL3A1 and MCL1 in Osteosarcoma. Med Sci Monit. 2018;24:8812‐8821. Published 2018 Dec 6. doi:10.12659/MSM.911972
Han XG, Mo HM, Liu XQ, et al. TIMP3 Overexpression Improves the Sensitivity of Osteosarcoma to Cisplatin by Reducing IL-6 Production. Front Genet. 2018;9:135. Published 2018 Apr 20. doi:10.3389/fgene.2018.00135
Sakai LY, Keene DR, Renard M, De Backer J. FBN1: The disease-causing gene for Marfan syndrome and other genetic disorders. Gene. 2016;591(1):279‐291. doi:10.1016/j.gene.2016.07.033
Cheng SW, Luk HM, Chu YWY, et al. A report of three families with FBN1-related acromelic dysplasias and review of literature for genotype-phenotype correlation in geleophysic dysplasia. Eur J Med Genet. 2018;61(4):219‐224. doi:10.1016/j.ejmg.2017.11.018
Yao X, Wu J, Chen J. Novel FBN1 Heterozygous Mutations Identified in Chinese Families with Marfan Syndrome. Ann Clin Lab Sci. 2019;49(4):539‐545.
Kajiwara Y, Wang E, Wang M, et al. GJA1 (connexin43) is a key regulator of Alzheimer's disease pathogenesis. Acta Neuropathol Commun. 2018;6(1):144. Published 2018 Dec 21. doi:10.1186/s40478-018-0642-x
Basheer WA, Fu Y, Shimura D, et al. Stress response protein GJA1-20k promotes mitochondrial biogenesis, metabolic quiescence, and cardioprotection against ischemia/reperfusion injury. JCI Insight. 2018;3(20):e121900. Published 2018 Oct 18. doi:10.1172/jci.insight.121900
James CC, Zeitz MJ, Calhoun PJ, Lamouille S, Smyth JW. Altered translation initiation of Gja1 limits gap junction formation during epithelial-mesenchymal transition. Mol Biol Cell. 2018;29(7):797‐808. doi:10.1091/mbc.E17-06-0406
Wang BJ, Chi KP, Shen RL, et al. TGFBI Promotes Tumor Growth and is Associated with Poor Prognosis in Oral Squamous Cell Carcinoma. J Cancer. 2019;10(20):4902‐4912. Published 2019 Aug 27. doi:10.7150/jca.29958
Palomeras S, Diaz-Lagares Á, Viñas G, et al. Epigenetic silencing of TGFBI confers resistance to trastuzumab in human breast cancer. Breast Cancer Res. 2019;21(1):79. Published 2019 Jul 5. doi:10.1186/s13058-019-1160-x
Seok Y, Lee WK, Park JY, Kim DS. TGFBI Promoter Methylation is Associated with Poor Prognosis in Lung Adenocarcinoma Patients. Mol Cells. 2019;42(2):161‐165. doi:10.14348/molcells.2018.0322
Hu J, Duan B, Jiang W, Fu S, Gao H, Lu L. Epidermal growth factor-containing fibulin-like extracellular matrix protein 1 (EFEMP1) suppressed the growth of hepatocellular carcinoma cells by promoting Semaphorin 3B(SEMA3B). Cancer Med. 2019;8(6):3152‐3166. doi:10.1002/cam4.2144
Yang T, Zhang H, Qiu H, et al. EFEMP1 is repressed by estrogen and inhibits the epithelial-mesenchymal transition via Wnt/β-catenin signaling in endometrial carcinoma. Oncotarget. 2016;7(18):25712‐25725. doi:10.18632/oncotarget.8263
Sachdeva R, Wu M, Smiljanic S, et al. ID1 Is Critical for Tumorigenesis and Regulates Chemoresistance in Glioblastoma. Cancer Res. 2019;79(16):4057‐4071. doi:10.1158/0008-5472.CAN-18-1357
Yin X, Tang B, Li JH, et al. ID1 promotes hepatocellular carcinoma proliferation and confers chemoresistance to oxaliplatin by activating pentose phosphate pathway. J Exp Clin Cancer Res. 2017;36(1):166. Published 2017 Nov 23. doi:10.1186/s13046-017-0637-7
Huang YH, Hu J, Chen F, et al. ID1 Mediates Escape from TGFβ Tumor Suppression in Pancreatic Cancer. Cancer Discov. 2020;10(1):142‐157. doi:10.1158/2159-8290.CD-19-0529

supplymentaryfile.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

Identification of critical genes and pathways involved in intervertebral disc degeneration (IDD) An integrated bioinformatics analysis

Status:

Version 1

Abstract

Figures

Background

Methods

Results

Discussion

Conclusion

Abbreviation

Declarations

References

Supplementary Files

Status:

Version 1