Expression and network analysis of YBX1 interactors for identification of new drug targets in lung adenocarcinoma.

Y-Box Binding protein 1 (YBX-1) is known to be involved in various types of cancers. It's interactors also play major role in various cellular functions. Present work aimed to study the expression profile of the YBX-1 interactors during lung adenocarcinoma (LUAD). The differential expression analysis involved 57 genes from 95 lung adenocarcinoma samples, construction of gene network and topology analysis. A Total of 43 genes were found to be differentially expressed from which 17 genes were found to be down regulated and 26 genes were up-regulated. We observed that Polyadenylate-binding protein 1 (PABPC1), a protein involved in YBX1 translation, is highly correlated with YBX1. The interaction network analysis for a differentially expressed non-coding RNA Growth Arrest Specific 5 (GAS5) suggests that two proteins namely, Growth Arrest Specific 2 (GAS2) and Peripheral myelin protein 22 (PMP22) are potentially involved in LUAD progression. The network analysis and differential expression suggests that Collagen type 1 alpha 2 (COL1A2) can be potential biomarker and target for LUAD.

Introduction

Lung Adenocarcinoma, a subtype of non-small cell lung cancer is the most pervasive among lung cancers leading to the death of millions of people each year 1, 2. Though the recent therapy for the lung adenocarcinoma with mutated EGFR and rearranged ALK have been significant 3, 4, but the other driving force for the progress of lung adenocarcinoma have not been deciphered much. The type of interaction among the gene and its environment decides their role in disease progression 5, 6. Understanding the gene interactions and its importance in the regulation mechanism is necessary to identify a potential target for therapeutic application.

In this study, we analyzed the expression of Y- box binding protein 1 (YBX1) and its interactors in lung adenocarcinoma (LUAD). YBX1 belongs to the Y-box binding protein family with a highly conserved cold shock domain and known to be involved in various eukaryotic cellular mechanisms 7-9. YBX1 is usually involved in RNA splicing and translational mechanism in cytoplasm but, it translocates to nucleus and gets involved in transcriptional regulation during stress condition 10-13. It is associated with drug resistance, cell proliferation and cell death 14. The transcription factor of YBX1 binds to the E-box for regulating it 15.

YBX1 is over expressed in various cancers and it has been suggested as biomarker in prognosis of various cancer types 16, 17. Multi drug resistance and metastasis are major reason for it's over expression 18, 19. YBX1 interactors are known to be involved in various cellular mechanisms ranging from cell signaling, DNA and protein repair mechanism to transcriptional regulation 20-26. Genes which are known to be involved in regulation of YBX1 as well as those regulated by YBX1 has been listed in (Table 1).

Table 1

Genes regulates and regulated by YBX1

Gene- Id	Gene	Reference
Genes Up- Regulates YBX1 | 4904
TP73 | 7161	Tumor Protein 73	17
MYC | 4609	v-myc avian myelocytomatosis viral oncogene homolog
MAX | 4149	MYC Associated Factor X
TWIST1 | 7291	Twist basic helix-loop-helix transcription factor 1	14
PABPC1 | 26986	Poly(A) binding protein, cytoplasmic 1	18
GATA1 | 2623	GATA binding protein 1	19
GATA2 | 2624	GATA binding protein 2
PTGER1 | 5731	Prostaglandin E receptor 1 (subtype EP1)	20
SHH | 6469	Sonic Hedgehog	21
Genes Down- Regulates YBX1 | 4094
FOXO3 | 2309	Forkhead Box 03	15
ILK | 3611	Integrin linked kinase	16
TGFB1 | 7040	Transforming growth factor, beta 1	22
C1QBP | 708	Complement component 1, q subcomponent binding protein	23
GAS5 | 60674	Long non coding RNA growth arrest specific transcript 5	24
KAT2B | 8850	K (lysine) acetyltransferase 2B	25
Genes Activated by YBX1 | 4904
CCL5 | 6352	Chemokine (C-C motif) ligand 5	26,27
CD44 | 960	CD44 molecule	28
ITGA6 | 3655	Integrin, alpha 6
MMP2 | 4313	Matrix metallopeptidase 2	29,30
POLA1 | 5422	Polymerase (DNA directed), alpha 1, catalytic subunit	31
EGFR | 1956	Epidermal growth factor receptor	32
ERBB2 | 2064	V-erb-b2 avian erythroblastic leukemia viral oncogene homolog 2
MET | 4233	Met proto-oncogene	33
ABCB1 | 5243	ATP-binding cassette, sub-family B (MDR/TAP), member 1	34,35
MVP | 9961	Major vault protein	36
PDGFB | 5155	Platelet-derived growth factor beta polypeptide	37
PTPN1 | 5770	Protein tyrosine phosphatase, non-receptor type 1	38
SMAD7 | 4092	SMAD family member 7	39
CCNA1 | 8900	Cyclin A1	40
CCNA2 | 890	Cyclin A2
CCNB1 | 891	Cyclin B1
PIK3CA | 5290	Phosphatidylinositol-4,5-bisphosphate 3-kinase, catalytic subunit alpha	41
Genes Repressed by YBX1 | 4904
ACTA1 | 58	Actin, alpha 1, skeletal muscle	42,43
COL1A1 | 1277	Collagen, type I, alpha 1	44,45
COL1A2 | 1278	Collagen, type I, alpha 2	46
CPS1 | 1373	Carbamoyl-phosphate synthase 1, mitochondrial	47
FAS | 355	Fas cell surface death receptor	5
CSF2 | 1437	Colony stimulating factor 2 (granulocyte-macrophage)	48,49
HSPA5 | 3309	Heat shock 70kda protein 5	50
MMP12 | 4321	Matrix metallopeptidase 12	51
MMP13 | 4322	Matrix metallopeptidase 13	52
HLA-A | 3105	Major histocompatibility complex, class I, A	53
HLA-B | 3106	Major histocompatibility complex, class I, B
HLA-C | 3107	Major histocompatibility complex, class I, C
HLA-E | 3133	Major histocompatibility complex, class I, E
HLA-F | 3134	Major histocompatibility complex, class I, F
HLA-G | 3135	Major histocompatibility complex, class I, G
HLA-DRA | 3122	Major histocompatibility complex, class II, DR alpha	54-56
B2M | 567	Beta-2-microglobulin
HLA-DQB1 | 3119	Major histocompatibility complex, class II, DQ beta 1
ABCC2 | 1244	ATP-binding cassette, sub-family C (CFTR/MRP), member 2	57
CDKN1A | 1026	Cyclin-dependent kinase inhibitor 1A (p21, Cip1)	58
TP53 | 7157	Tumor protein p53	59
TSHR | 7253	Thyroid stimulating hormone receptor	60
VEGFA | 7422	Vascular endothelial growth factor A	61,62
SOX2 | 6657	SRY (sex determining region Y)-box 2	63

In this study we investigated the role of YBX1 and its interactors in Lung adenocarcinoma by looking at their expression profile and constructed gene regulatory network in order to decipher their importance in network formation by a system biological approach.

Materials and Methods

RNA Seq data and Sample Quality Analysis

The Cancer Genome Atlas (TCGA) RNA-Seq level 3 data for YBX1 and its interactors belonging to 116 normal and normal matched tumor samples in LUAD were downloaded from Broad Genome Data Analysis Center (GDAC) Firehose site (https://gdac.broadinstitute.org/). RNA Sequence by Expectation Maximization (RSEM) 69 counts for 57 genes of our interest including YBX-1 were used for differential expression analysis.

Principle Component Analysis (PCA) and hierarchical clustering with log transformed datasets of the samples were performed using R-Bioconductor 70,71 package DESeq2 72 to estimate sample dispersion and to filter out the outlier.

Differential Expression (DE) Analysis and Correlation between genes

R-Bioconductor package DESeq2 was used to carry out differential expression analysis with nbionom Wald test for calculating logarithmic fold change and Benjamini- Hochberg method 73 for estimating adjusted-p value. Genes with adjusted-p value less than 0.05 were considered as differentially expressed. Correlation between differentially expressed genes was calculated with Pearson Correlation method using R package Hmisc 74 and the plot for the correlation coefficient was constructed using R package Corrplot 75.

Gene Regulatory Network Construction and Analysis

Gene regulatory network was constructed for the 56 genes omitting GAS5 using the GeneMANIA 76 with maximum resultant gene and attributes to enrich were 30 and 215 respectively. And the constructed network was analyzed for its network topology by looking at its parameters like closeness centrality, betweenness centrality and degree with network analysis tool in Cytoscape 3.4 77. RAIN (RNA - protein Association and Interaction Network) 78 was used to identify interactors for the long non-coding RNA (lncRNA) GAS5. And Search Tool for Retrieval of Interacting Genes (STRING) database 79 was used to analyze the network and the interactions with confidence level of 0.7 were considered significant. Network with not more than 10 interactors for the GAS5 was constructed.

Gene Set Enrichment Analysis

Functional enrichment for the genes common to both differentially expressed gene sets and hub genes from network was done by using Database for Annotation, Visualization and Integrated Discovery (DAVID) 80. Gene Ontology (GO) database 81 provides the annotation for the gene set and Reactome pathway 82 illustrate the pathways in which the genes were involved. With DAVID tool we analyzed the gene set in both GO and Reactome pathway with p-value < 0.05 and gene-count > 3 as condition for enrichment analysis.

Results

Assessment of Sample Quality and Filtering

In this study, we initially took 116 LUAD samples with TCGA level-3 data. With Principle Component Analysis (PCA) and hierarchical clustering removed samples whichever showed variation within its group (Figure S1). And we plotted again PCA and hierarchical clustering for the final 95 samples with 57 genes. After filtering, two methods now able to distinguish the samples based on their conditions- normal or tumor (Figure 1).

Figure 1

Quality of the samples after filtering. (1A) Heatmap shows the unsupervised hierarchical clustering of the normal and tumor samples in LUAD after filtration. The row represents the genes and column represent the samples. Normal and tumor samples are clustered within their group based on their Euclidean distance. (1B) Principle Component Analysis for the filtered samples shows that first principle component (PC1) separates the normal and tumor samples. In both analysis, samples are found to be grouped within their type.

Differential Expression Analysis and Correlation between Genes

The 43 genes which have adjusted-p value less than 0.05 were selected as differentially expressed with GATA1 having a high negative fold change (-1.91476) and MMP13 having a high positive fold change (5.81790). And correlation matrix between the differentially expressed genes was constructed and the correlation coefficient with p-value less than 0.05 were considered significant (Figure 2). YBX1 showed a maximum positive correlation with C1QBP (0.592) and negative correlation with HLA-E (-0.281) (Table S1).

Figure 2

Correlation plot between the differentially expressed genes in LUAD. Correlation with the significance p- value (<0.05) are shown. Insignificant correlation are left blank without colors.

Gene Network Construction and Analysis

Gene network was constructed and analyzed (Figure 3), it has 86 nodes and 3611 interactions (Table 2). Top 10 genes involved in the network formation were selected based on the topology parameter-degree, betweenness centrality and closeness centrality (Table 3). The attributes from GeneMANIA showed cancer pathways hold 15.34% in network with involvement of 17 genes (Table 4).

Figure 3

Gene network constructed and analyzed with Cytoscape. Directed root network built from the gene interaction data obtained from GeneMANIA. Network has 86 nodes with 3611 interactions. Figure , 3B and 3C highlights the top 10 genes with closeness centrality, betweenness centrality and degree respectively. CDKN1A is the common gene between closeness centrality and betweenness centrality. Five genes namely, HLA-DRA, HLA-A, HLA-DRB1, HLA-DMB and B2M showed highest betweenness centrality and degree.

Table 2

Interaction type and number of interactions in the network constructed by GeneMANIA

Type of Interactions	Number of Interactions
Co-expression	1761
Co-localization	116
Genetic Interaction	98
Pathways	182
Physical Interaction	462
Predicted	13
Shared Protein Domain	979

Table 3

Central genes in Gene Regulatory network constructed using GeneMANIA and analyzed by Cytoscape- Network Analysis Tool. Top 10 gene with high Closeness Centrality, Betweenness Centrality and Degree.

Genes	Closeness Centrality	Genes	Betweenness Centrality	Genes	Degree
MMP13	1.000	HLA-DRA	0.018	HLA-DRA	317
YBX-1	1.000	B2M	0.017	HLA-DMA	274
TSHR	1.000	HLA-A	0.015	HLA-A	266
ABCC2	0.800	HLA-C	0.014	HLA-DPB1	266
KATA2B	0.750	CDKN1A	0.014	HLA-DPA1	257
ABCB1	0.750	ILK	0.014	B2M	255
COL1A2	0.700	HLA-B	0.010	HLA-G	252
CDKN1A	0.689	HLA-DMB	0.010	HLA-DMB	246
MMP 12	0.667	VEGFA	0.009	HLA-F	241
MYC	0.667	HLA-DRB1	0.008	HLA-DRB1	226

Table 4

Genes involved in Cancer Pathways, known from GENEMANIA network attributes.

Gene Id	Gene
CCNAI | 8900	Cyclin A
VEGFA | 7422	vascular endothelial growth factor A
MMP2 | 4313	Matrix Metallopeptidase 2
MET | 4233	Met proto-oncogene
MAX | 4149	MYC associated factor X
SHH | 6469	Sonic hedgehog
TGFB1 | 7040	Transforming Growth Factor, Beta 1
ITGA6 | 3655	Integrin, alpha 6
FAS | 355	Fas cell surface death receptor
TP53 | 7157	Tumor Protein p53
PDGFB | 5155	Platelet-derived growth factor beta polypeptide
CDKN1A | 1026	Cyclin-dependent kinase inhibitor 1A
ERBB2 | 2064	v-erb-b2 avian erythroblastic leukemia viral oncogene homolog 2
PIK3CA | 5290	Phosphatidylinositol-4,5-bisphosphate 3-kinase, catalytic subunit alpha
EGFR | 1956	Epidermal growth factor receptor
CDKN1B |	Cyclin-dependent kinase inhibitor 1B

As GeneMANIA construct network and interactions for the coding gene, it was not able to predict for GAS5. So RAIN database which identify the interactors for non-coding RNA was used to identify interactors of GAS5. Total 10 interactors were found and the network was constructed between the genes with confidence level of 0.7 using STRING database (Figure 4)

Figure 4

Interactors of GAS5 derived from RAIN database. Interaction with confidence level > 0.7 is considered as significant. Among the interactors GAS2 and PMP22 are the proteins involved in interaction with GAS5.

The annotation of the gene with GO and Reactome pathways using DAVID resulted in a single cluster with genes enriched in GO- cellular components of proteinaceous extracellular matrix and extracellular region. GO- Biological process contained the collagen catabolic process and Reactome pathways included Generic transcriptional regulation. Only 3 genes were involved in the enrichment when the condition for enrichment was made high, namely- COL1A2, MMP12, MMP13. The enriched terms are shown in Table .

Discussion

Despite of various studies carried out on YBX1, role of its interactors in lung adenocarcinoma have been less explored. YBX1 have been involved in the regulation of tumor progression in lung adenocarcinoma 83, 84, so understanding its interactions with other gene become a necessary to know their mechanism in disease progression and drug discovery. Our present study finds that YBX1 is up regulated with a low fold change (0.224). GATA 1 and GATA 2 which upregulate YBX1 level in erythroid cells 20 were found to be down regulated that indicates that they don't play major role in regulation of YBX1 in lung adenocarcinoma.

Among the genes which upregulate YBX1, PABPC1 which is involved in the translation of YBX1 mRNA 85 was found to be having high correlation (> 0.4) with YBX1 and it was significantly expressed with a fold change of 1.259. C1QBP which is known to be involved in prostate cancer progression 86 and a highly negative regulator of YBX1 in renal cell carcinoma 87 is highly correlated with YBX1 along with GAS5, a long non-coding RNA which involves in regulating cell death in prostate cancer 88. Among the GAS5 interactors identified from RAIN database there are two proteins namely Growth Arrest Specific 2 (GAS2) and Peripheral myelin protein 22 (PMP22) which are known to be over expressed in colorectal cancer cell and breast cancer patients 89, 90 respectively. However the evidence regarding their role in LUAD disease progression needs to be explored for being a potential biomarker in LUAD identification and for therapeutic application.

Comparison of the differentially expressed genes and the hub genes in the network formation revealed that MMP13, ABCC2, MMP12, TSHR, COL1A2, PABPC1 and YBX1 are common among the both groups. Except for PABPC1, rest of the genes are usually repressed by YBX1, but in our study they are among the highly up regulated genes making a need for further study to understand their relationship between YBX1 in LUAD. Among the three genes involved in the enrichment, MMP12 and MMP13 are reported to be involved in the lung cancer in earlier studies 90- 94. Significantly COL1A2 is involved in progress of gastric cancer 95 and in head and neck cancer 96. It can serve as a potential biomarker in LUAD and its involvement in the disease progress need to be explored.

Conclusion

The study aimed to investigate the role of YBX1 and its interactors in lung adenocarcinoma by looking at their differential expression and network topology study. The study showed PABPC1 can be potential target in regulating YBX1. The lncRNA GAS5, whose role in the LUAD need to be explored and it can be a potential biomarker along with COL1A2.

Supplementary figures and tables.

Click here for additional data file.

Table 5

Gene Enrichment Analysis- Gene Ontology (GO) and Reactome Pathways (RP) enriched by genes common to differential expression and hub gene in the gene regulatory network.

Group	Term	Gene Count	p-value
GO-BP	GO:0030574~collagen catabolic process	3	1.419E-4
GO-CC	GO:0005578~proteinaceous extracellular matrix	3	0.002
GO-CC	GO:0005576~extracellular region	3	0.065
RP	R-HSA-1442490:R-HSA-1442490	3	4.829E-4