- Research article
- Open Access
Analysis of RNA expression of normal and cancer tissues reveals high correlation of COP9 gene expression with respiratory chain complex components
BMC Genomicsvolume 17, Article number: 983 (2016)
The COP9 signalosome, composed of eight subunits, is implicated in cancer genetics with its deneddylase activity to modulate cellular concentration of oncogenic proteins such as IkB and TGFβ. However, its function in the normal cell physiology remains elusive. Primarily focusing on gene expression data of the normal tissues of the head and neck, the cancer genome atlas (TCGA) database was used to identify groups of genes that were expressed synergistically with the COP9 genes, particularly with the COPS5 (CSN5), which possesses the catalytic activity of COP9.
Expressions of seven of the COP9 genes (COPS2, COPS3, COPS4, COPS5, COPS6, COPS7A, and COPS8) were found to be highly synergistic in the normal tissues. In contrast, the tumor tissues decreased the coordinated expression pattern of COP9 genes. Pathway analysis revealed a high coordination of the expression of the COPS5 and the other COP9 genes with mitochondria-related functional pathways, including genes encoding the respiratory chain complex.
The results indicate that mRNA expression data for the matched normal tissues available in TCGA are statistically reliable, and are highly useful to assess novel associations of genes with functional pathways in normal physiology as well as in the cancer tissues. This study revealed the significant correlation between the expressions of the COP9 genes and those related to the mitochondrial activity.
The constitutive Photomorphogenesis (CSN) 9, or COP9 signalosome, is a critical multi-functional protein complex in cells [1, 2]. COP9, originally discovered in Arabidopsis thalania as a negative regulator of photomorphogenesis, was later found to be highly conserved across many species [3, 4]. COP9 is a 450–550 kDA holocomplex composed of eight subunits whose official gene symbols are listed in Table 1 [1, 5]. COP9 has catalytic activity to remove Nedd8, a small-ubiquitin like modifier, to regulate cullin-RING ubiquitin E3 ligases, and thus protein degradation pathways mediated by cullin complexes . COP9 regulates phosphorylation of proteins such as tumor suppressor p53, and neddylation essential to ubiquitin-mediated proteolysis of key proteins including tumor suppressor HIF-1α [6–10]. Thus, COP9 has been an important focus in cancer genetics [5, 10–13].
Despite intense studies in the past, the exact essential role of COP9 still remains unknown. This is largely due to COP9’s promiscuous property to interact with a variety of cullin complexes and others, and to affect stabilities of a number of cellular factors. Additionally, subforms of COP9 complexes were reported [2, 14, 15]. Dubiel et al. proposed that such flexible conformation enables the broad range of substrates.
Among the eight subunits, COPS5 is the most extensively studied subunit that has numerous critical functions including deneddylation reactions [1, 10]. It is commonly over-expressed in cancers leading to increased deneddylation within cells . Other subunits have been shown to possess unique functions. For example, CSN1 (GPS1) is involved in stabilization of p53 [10, 13]. Over-expression of COPS2 is linked to chromosome instability [10, 13]. Both COPS2 and COPS3 knockout are embryonic lethal at E3.5 and E8.5 respectively [10, 13]. COPS6 is frequently over-expressed in cancer . COPS8 knockout is also embryonic lethal (E7.5) with impaired growth and differentiation . Although these reports convincingly indicate that all the COP9 subunits have essential functions in mammalian cells, how these reported functions contribute to COP9 activity remains largely unknown.
Identification of functional groups with which COP9 signalosome may interact in normal physiology may provide an important clue to understand how dysregulation of COP9 may facilitate tumorigenesis.
Deficiency in even one subunit leads to COP9 destabilization. Down-regulation of COPS5 in mouse embryonic fibroblasts drastically decreased the stability of COPS1, COPS3, and COPS8 [17–19]. Therefore, it is plausible that coordinated expression of COP9 genes is needed to maintain necessary amounts of COP9 and imbalance may be associated with disease . To elucidate this possibility, a statistically adequate number of normal tissues must be analyzed with precision for their gene expressions. Using cultured cell lines is not appropriate for this investigation, because the cell lines established in vitro are all transformed and thus expressions of the COP9 genes are likely altered from cells in normal in vivo conditions.
The Cancer Genome Atlas (TCGA) has procured specimens from more than 500 human tumor tissues from various cancers for comprehensive analyses of tumor tissues and advanced cancer genomics. Importantly, matched normal tissues have been analyzed for about 10% of each cancer type. Using the RNAseq expression data in TCGA, we elucidated the possibility that expression of the COP9 genes are synchronized in the normal cells. Indeed the COP9 gene expressions were in synergy in the normal tissues, but this fine coordination was diminished in the tumor tissues. Furthermore, the analysis led to the unexpected observation that the expressions of the COP9 genes were highly synchronized with those involved in mitochondrial biogenesis, particularly those of respiratory chain complex.
Nomenclature of COP9 genes
Although the gene names from CSN1 to CSN8 are intuitive and have been accepted as the names for the COP9 subunits , the TCGA uses the official gene symbols. For the simplicity, this study used the official gene symbols as the names of COP9 subunits (Table 1).
Data obtained from TCGA
The mRNA expression data of normal and tumor tissues were obtained through TCGA’s online data portal site (http://cancergenome.nih.gov). The cancer tissues chosen for this study are 522 cases of head and neck squamous cell carcinoma with 44 normal tissues, as well as lung adenocarcinoma (50 normal and 502 tumor). Normalized gene expression results (mRNASeq2 level3) were collected. Each set of mRNA expression data contained 20,531 gene entries, which were linked to annotations through DAVID functional gene annotation database (Database for Annotation, Visualization and Integrated Discovery; https://david.ncifcrf.gov), yielding 20,154 analyzable gene lists. The RNAseq data were linked to information about the chromosome cytoband and gene ontology (GOTERM_BP_FAT, GOTERM_CC_FAT, GOTERM_MF_FAT, and KEGG_PATHWAY).
Similarly, the de-identified clinical information including anatomical sites and smoking history were retrieved from TCGA. Data processing and analysis were carried out using R (ver 3.1.2) and series of command line programs developed with the Swift language using Xcode ver. 6.
Generation of correlation coefficient table specific for individual genes
Gene expression data of normal and tumor tissues from donors were assembled into a single table file consisting of columns of tissues and rows of genes. Genes that show zero mRNA values in more than 5% of the tissues analyzed were excluded unless noted otherwise. Pearson’s correlation coefficient (r) values between expressions of a target gene (e.g., COPS5) and each of all the other genes were determined, and the genes were then sorted based on the r values. Significance of r values were assessed with p values calculated with cor.test function in R. Pair and boxplots were drawn for particular gene sets using R.
Expressional correlation of genes closely localized on the chromosomes
Generation of cytoband plots
Genes ranked within 500 highest |r| values were pooled, and the cytoband information was processed into float values by the simple formula as following. For a gene with cytoband of ApB.N or AqB.N, a float value was generated by A - (B.N)/100 or A + (B.N)/100, respectively. For example, 8p11.2 was converted to 8–11.2/100 = 7.888, and 8q11.2 was to 8.112. Values of (x, y) = (cytoband, r) was plotted using R. X and Y chromosomes were shown in one column as all analyzed genes mapped in X chromosome were also found in Y chromosome.
Consideration of exact transcription start sites
Transcriptional start sites of the COP9 genes were retrieved from the Ensemble database. For a particular geneA in the RNAseq data, a group of genes (Gi) localized within 1 Mbp of the geneA were pooled (the number of genes in Gi = Nh). Then, an enrichment score (ES) for the Gi was calculated based on the correlation coefficient table for geneA . Prior to this calculation, geneA itself (r = 1) was removed from the table. Separately, genes localized on the same chromosome but outside of the 1 Mbp of geneA were pooled (Go) and used to calculate control mean ES and standard deviation (sd) for the Nh number of genes that were randomly selected from the Go gene pool. This calculation was repeated for 1000 times to obtain mean and sd values which were used to generate a p value based on pnorm function with R. To compare these ES values with other genes, the above calculation was performed with all available genes in the RNAseq data, and the p-values were used to calculate false discovery rates (FDR; Benjamini–Hochberg) using p.adjust function in R.
Interpretation of functional gene annotation
The top 1000 genes with the highest |r| values were pooled with each gene accompanying gene ontology (GO) and KEGG pathway information. Over-represented pathways are sorted by the p values based on binomial distribution: 1 - pbinom(n, 1000, μ) (in R), where n is the actual count occurring in the 1000 genes, and μ is the average of finding the particular pathways, and calculated with the whole gene set of 20,154.
Derivation of enrichment score (ES) for KEGG functional pathways
ES scores were calculated by the method described previously . For example, an ES for COPS5 to oxidative phosphorylation pathway was calculated using the COPS5 correlation coefficient table with gene annotations based on the association with “hsa00190” (KEGG i.d. for oxidative phosphorylation) to generate Phit - Pmiss. ES was provided as the maximum value of Phit - Pmiss as described . A histogram of ES values was calculated with randomly chosen Nh genes for iteration of 1000, and the distribution was used to assess the statistical significance.
ESs and corresponding p values were calculated for all genes and then false discovery rates were determined as described above.
A single gene expression analysis
COPS5 (official gene symbol of CSN5) is responsible for the deneddylation activity of the COP9 signalosome, and also known as Jab1 to activate c-Jun (Jun-activating/binding protein). Thus, most COP9 studies have focused on the activity of COPS5  and its influence on biological activities including tumorigenesis. Because of the pivotal role of COPS5 for the COP9 signalosome in mammalian cells, this study primarily focused on COPS5 to investigate the possibility of coordinated expression of the COP9 subunits in the normal and tumor tissues. Our primary interest was the head and neck squamous cell carcinoma (HNSCC), where effects of smoking and influences of major genetic factors, including p53, p16, AKT1 and K-Ras could be addressed due to its established cancer etiology involving these genes. In addition, HNSCC data in TCGA included a few distinct anatomical sites including oral cavity, tongue and larynx, and thus effects of tissue sites for the gene expression could be elucidated with the clarified clinical record. For the simplicity, the official gene symbols will be used in this study instead of more commonly used protein names (Table 1).
Expression of the COPS5 gene was examined using the data of 44 normal head and neck tissues available in TCGA. COPS5 expression in normal tissues showed mean value of 1261.0 with standard deviation (sd) 328.83 and coefficient of variation (sd/mean) 0.26. As a comparison, mean coefficient of variation of the 15,660 gene set was calculated, and found to be 0.693 (±0.575). Therefore, the COPS5 gene expression appeared to be relatively stable among normal oral tissues.
Coordinated expression of the COP9 genes in normal tissues
Expression patterns of COPS5 and other COP9 genes were elucidated by comparing correlation coefficients (expressed as r thereafter), using the RNA expression data (RNAseqV2 level 3). High positive correlations were found among the COP9 gene expressions. Specifically, relative to COPS5, all the other COP9 genes except COPS7B showed statistically significant positive r values (Table 2 and Fig. 1). GPS1 (CSN1) showed a smaller r value than the others, namely COPS2, COPS3, COPS4, COPS6, COPS7A, and COPS8. The lower correlation of gene expression of GPS1 vs. COPS5 was found to be a general characteristic through our analyses in this study. However, the positive correlation between GPS1 and COPS5 was still statistically significant (p = 0.011). In addition to the low r value relative to the other COP9 genes, the level of COPS7B expression was significantly lower compared to COPS7A (Additional file 1: Figure S1).
These results indicated that the expressions of COP9 genes were synchronized in human normal oral tissues. Leppert et al. recently observed that mRNA levels of COP9 genes were regulated by let-7 miRNA . To probe the importance of the regulatory mechanism by miRNA, the TCGA miRNA dataset was used to determine the gene expression correlation of COP9 genes and the let-7 miRNA species (Additional file 2: Table S1). There were statistically significant negative correlations between let-7 isoforms (7a, 7b, 7c, 7e, 7f-2, 7 g) and COP9 genes (GPS1, COPS5, COPS6, and COPS8), which was consistent with the previous work.
Loss of synergistic expressions of the COP9 genes in cancer tissues
To determine whether the synergistic expressions of the COP9 genes are maintained in the cancer tissues, 42 tumor tissues matched to the normal cases were pooled (2 normal cases could not be matched based on the id tag), and determined their correlation coefficients (Table 3). High positive r values were observed among COP9 genes in normal tissues as expected. In contrast, the correlations of expressions among COP9 genes diminished in the matched tumor tissues (Table 3 and Additional file 2: Table S2). Moreover, the negative effect of let-7 miRNA observed in the normal tissue dataset was diminished in the matched tumor tissues (Additional file 2: Table S1). These results demonstrated that the coordinated expression of COP9 genes was diminished in these tumor tissues.
Validation of the results with subgroups (age, anatomical sites, and smoking history)
HNSCC is an age-related disease. The majority of specimens in TCGA are from elderly patients (median = 61; the first quartile = 53, the third quartile = 69), and less than 4% of the patients are younger than 40 years old. Thus, the loss of the coordination of the COP9 gene expressions may be due to the age-related degeneration of the coordination. However, we analyzed the COP9 expression in the 31 tumor tissues that were 42 years old or younger to the 30 tumor tissues that were 80 years old or older (Additional file 2: Table S3 A and B). With this analysis, we found similar loss of the coordinated expressions in specimens from both the young and elderly patients. In addition, there were no significant differences of the expression levels of each COP9 gene between the young and aged groups (Additional file 2: Table S3C). These results support that the age is not a factor causing the low coordination of the COP9 genes in the tumor specimens.
The oral tissues in TCGA were taken from different anatomical sites (oral cavity, tongue, etc.). The synergistic expression of COP9 subunit genes may vary depending on the exact site of oral tissues. To assess this possibility, the data set was arranged into subgroups based on the sites. The major groups were the oral cavity (14 specimens) and the oral tongue (13 specimens). The data sets of the two subgroups were used to determine correlation coefficients of the COP9 genes. These values were similar to those of the entire dataset (Additional file 2: Table S4 and Additional file 2: Table S5), indicating that the synergistic expression of the COP9 genes are maintained regardless of the sites of the oral tissues. Similarly, there appeared to be no clear effects of smoking, because analyzing the normal tissues from eleven current smokers identified a highly similar synergistic RNA expressions of COP9 genes (Additional file 2: Table S6).
Chromosome mapping of genes of high synergy with the COPS5 expression
To elucidate the tumor-associated disintegration of the COP9 gene expression in detail, we focused on the COPS5 gene. We generated a “correlation table” of the COPS5, in which expressional correlation coefficients relative to the COPS5 gene (r COPS5 ) were calculated for all the analyzable genes (about 14,000 genes). By cross-referring the results to gene ontology provided with DAVID, we noticed that genes with high r COPS5 of the tumor tissues were mapped in the same chromosomal location as the COPS5 gene (8q13.2, Table 1). We thus speculated that the loss of the coordination of the COP9 expressions in the tumor tissues was not due to unpredictably disorganized expressions of the COP9 subunit genes, but rather there was a dominant effect of the physical locations of the genes on their expressions. To illustrate this possibility, chromosomal locations of 500 genes with top | r COPS5 | values were plotted on chromosomes based on the cytoband information (Fig. 2, top). The genes with high |r COPS5| in the normal tissues were distributed through all the chromosomes; among the top 100 genes, only seven genes were mapped at 8q, near the COPS5 chromosomal location (8q13.2). In contrast, genes with high |r COPS5| in the tumor tissues were concentrated near COPS5’s chromosomal region (Fig. 2, bottom); 56 genes out of the top 100 genes were bound to 8q.
We sought to assess the chromosomal cis-effect on the COPS5 gene expression more rigorously. To this end, based on the method by Subramanian et al. we determined the enrichment score (ES) based on the r cops5 for the groups of genes that are localized within 1 mega base pairs (Mbp) from the COPS5, and obtained a p-value based on the null hypothesis that these ES values were not different from that calculated for genes outside 1 Mbp from the COPS5 . We have also determined such ES values for all the genes in the chromosome 8 to determine the false discovery rate (FDR, q). The ES value for the COPS5 with the tumor tissues (ES = 0.65, p = 0.02, q = 0.043) was significantly higher than that with the normal tissues (ES = 0.39, p = 0.27, q = 0.49). The procedure was applied for the other COP9 genes, and revealed significant links to the chromosomal cis-effect in the case of the tumor tissues (Additional file 2: Table S7). Therefore, we conclude that the COP9 expression in the tumors is disintegrated because the chromosomal cis-effect becomes a dominant force to express each subunit gene. Furthermore, these results indicated that analysis of the COP9 expressions in the normal tissues rather than tumor tissues is pivotal for understanding the COP9’s role in vivo, and such an analysis is possible using the RNAseq data in the normal tissues available from TCGA.
Pathway analysis based on KEGG and GO revealed association of COPS5 expression with mitochondrial pathways
Genes that belong to particular cell pathways may be over-represented in the COPS5 correlation coefficient table. Identifying such pathways should provide important clues as to COP9’s roles under the normal physiology. Based on the results above, such functional association should be elucidated with the data derived from normal tissues, and the large pool of high quality gene expression data in the TCGA database prompted us to investigate this possibility. Genes with the top 1000 |r COPS5| values in the normal tissues were pooled, and functional pathways available through GO and KEGG database were associated to each gene (Additional file 3: Table S10). The GO and KEGG pathways were then ranked based on the over-representation, and the 20 most over-represented pathways are shown in Table 4 (the expanded list is shown in Additional file 4: Table S11). Unexpectedly, the pathways related to mitochondrial functions and components were highly over-represented. For example, 61 genes belonging to the oxidative phosphorylation KEGG pathway (hsa00190) were among the genes with high expression correlation with COPS5 out of 106 total genes (p < 1.1e-16). A number of GO pathways involving mitochondrial functions were also found among the over-represented pathways. These included respiratory chain (GO:0070469), mitochondrial electron transport (GO:0006120), NADH dehydrogenase activity (GO:0003954), mitochondrial respiratory chain complex I (GO:0005747). In fact, all of the 20 most over-represented pathways are directly related to major mitochondrial and energy derivation functions. At the time this study was carried out, gene ontology database had not yet linked COPS5 to any of mitochondrial-related pathways. Expression values were examined in pair-plots for COPS5 and mitochondrial genes of the five highest r COPS5, NDUFB6, MRPS36, ATP5F1, TMEM126A, NDUFB3 (Fig. 3). The results illustrate high integrity of the data and the coordination of the expression of the COPS5 with the mitochondria related genes.
Dependency on the anatomical sites and smoking history were tested with the same subgroups as described above (Additional file 5: Table S12). The mRNA expression data from the two major anatomical sites, oral cavity and oral tongue, resulted in almost identical over-representation of the mitochondrial pathways with COPS5. The same conclusion was made with the subgroup containing only the current smokers.
Alteration of COPS5 expression coordination in cancer tissues
The gene expression coordination between the COPS5 gene and the mitochondrial pathways was elucidated with tumor tissues (Table 5 and Additional file 4: Table S11). Although a few GO pathways related to mitochondrial functions were over-represented, major pathways associated with the COPS5 expression were those involved in ribosome processing and RNA binding as well as membrane lumen pathways.
The possibility that tissues from a particular anatomical site with high dysregulation caused the bias was addressed. Oral cavity, tongue, and larynx are the major three sites that together contribute to about 80% of the donated head and neck cancer tissues. COPS5-associated GO and KEGG pathways were determined for the tissues from the particular sites, and confirmed that correlation of COPS5 expression with mitochondria associated genes were not as predominant as in the normal tissues.
To further support the finding, enrichment score (ES) was calculated for the genes in the oxidative phosphorylation KEGG pathway (hsa00190) using the RNAseq data of 44 normal head and neck tissues  (Fig. 4a). The obtained ES value, 0.797, was much higher than the average ES value based on randomly assigned gene pools (Fig. 4b).
In addition, ES values corresponding to all the other KEGG pathways were calculated for the 42 matched normal and tumor tissues (Table 6). The ES value of the oxidative phosphorylation pathway (hsa00190) found to be among the highest. The ubiquinone biosynthesis pathway generated the highest ES value, which is an essential factor in the respiratory chain and requires phenylalanine (generating the second highest ES) for its synthesis. Thus, the enrichment of genes involved in mitochondria-related pathway based on the COPS5 correlation was exceptionally high.
To elucidate the significance of the COP9 coordination with oxidative phosphorylation compared to other genes, we calculated ES values relative to the hsa00190 KEGG pathway for the entire gene list in the RNAseq data of the matched normal and tumor tissues to calculate FDR (Additional file 2: Table S9). Except for the COPS7B gene, all the COP9 subunit genes showed highly significant correlations with the genes in the oxidative phosphorylation pathway. Not surprisingly, most of the genes in the hsa00190 oxidative phosphorylation pathway resulted in the high ES/low FDR values, validating our approach to detect a linkage of a gene to the oxidative phosphorylation. Remarkably, the COP9 genes (except for COPS7B) were ranked comparably with the hsa00190 genes. The p as well as FDR values were increased for all the COP9 genes (but COPS7B) with the tumors, indicating significant decrease in the association of the COP9 expression to the oxidative phosphorylation related genes. Taken all together, we conclude that there is a high synergy between the expression of the COP9 and mitochondria functions, which implies that the normal physiological function of the COP9 may be inseparable from the mitochondrial activity.
Analysis of gene expression of other cancers in the TCGA database
Finally, the question was addressed whether the association of the expressions of COP9 genes with those of mitochondrial genes was unique in head and neck tissues. To probe this, the RNA expression data of TCGA for the lung squamous cell carcinoma tissues were analyzed. Tumors originating from the lung upper airway epithelia share similarity with HNSCC in its etiology and its involvement of smoking as an important risk factor. First, the synergistic expression pattern of COP9 genes in the normal tissues (n = 50) was examined, and found to be highly synergistic with one another (Table 7 and Additional file 1: Figure S2). Exceptions were GPS1 and COPS7B, which showed lower correlation coefficients to other COP9 factors. The high positive correlations were again lost in the tumor tissues. These results showed a remarkable coincidence with those of head and neck normal and tumor tissues.
The functional pathway analysis with the entire gene list resulted in the high correlation of the COPS5 expression with the mitochondrial genes (Additional file 6: Table S13). Unlike the case of head and neck squamous cell carcinoma, however, the coordinated expression of COP9 and the mitochondrial genes was maintained well in the lung squamous cell carcinoma tissues (Additional file 6: Table S13).
Cancer genetics and analyses of the TCGA database mainly focus on the loss of genomic integrity. Mutations and single nucleotide polymorphism in DNA show the high rate of permanent alteration of the genome, epigenetic modification of the genome, RNA and miRNA expression alterations. TCGA data sets are available publicly for these types of analyses. The database is highly reliable and allowed us to examine cancer genomics with more than 500 tumor specimens .
In addition to the tumor specimens, data from about 40 to 50 normal tissues are made available for most types of cancers with a notable exception for glioblastoma. Although the data sizes are smaller than those for the cancer tissues, they are arguably one of the most comprehensive databases for genomics of human specimens of non-disease status with pathological and clinical scrutiny, which should provide valuable information for elucidating functions of a gene under the normal physiology. We hypothesized that a synergistic expression of genes may be observed in the normal tissues, when expressions of a group of genes need to be coordinated for their proper cellular activities. COP9 signalosome is an ideal subject to test this possibility, because all eight subunits are essential for the proper COP9 function; deficiency in one of the COP9 genes decreased stability of the other subunits and caused lethality [23, 25]. In this study the COP9 signalosome genes showed highly synergistic expressions in the normal head and neck tissues, except for GPS1 and COPS7B. Remarkably, these characteristics, including the lower r values for GPS1 and COPS7b, were maintained well even in the lung and breast normal tissues (data not shown). The coordinated expression of the COP9 genes decreased in the tumor tissues. Therefore, this coordination could not have been detected by analyzing tumor tissues. Likewise, because all cells cultured in vitro are under non-physiological conditions and thus transformed at some extent, it would have been difficult to obtain convincing evidence for the expression correlation of the COP9 subunit genes based on in vitro cell biological studies alone.
It is interesting to observe the lowest coordination of GPS1 expression compared to the other COP9 components (excluding COPS7B). Besides being a subunit of COP9, GPS1 is known for its function in suppression of G-protein and mitogen-activated signal transduction in mammalian cells . This COP9-independent function quite likely requires a regulation of the GPS1 gene separate from the other COP9 genes.
COPS7B not only had the lowest r value, but also significantly lower expression compared to COPS7A (Additional file 1: Figure S1). While COPS7A and COPS7B genes are encoded in different human chromosomes, chromosome 12 and 2 respectively, there is a high homology between COPS7A and COPS7B in the amino acid sequence (78% identical). Thus, COPS7B may have only a minor functional role for COP9, which may be compensated by COPS7A.
To our knowledge, using the exact physical map of the transcription start sites available from the Ensemble genome database has not been applied previously to elucidate effects of the relative locations of genes on the gene expression correlations. In this study, we used this approach to find the chromosomal locations of the genes as a dominant factor causing the loss of the coordinated expressions of the COP9 genes. Much fewer number of genes were under the control of this cis-effect in case of normal tissues, indicating that most genes in normal tissues are regulated independently of their chromosomal locations and perhaps more based on their functional groups, as supported by this study. Therefore, important implications derived from this analysis are that 1) analyzing tumor tissue data may not identify a normal function of a protein family such as COP9, and 2) the qualities of the normal tissue data available in TCGA are high enough to help identifying unknown expressional linkages with functional pathways. Although the primary focus of this study was COP9 expression profiling, we expanded the analysis of the influence of the chromosomal location for all the genes that were analyzable in the RNAseq data, and the data for cancer genes (hsa05200) are shown in Additional file 7: Table S14. It may be interesting to elucidate whether these proto-oncogenes have tight association with the chromosomal locations in other types of tumor tissues available in TCGA.
The observation that the expression of the COPS5 gene (and therefore the other COP9 genes) was synchronized with those involved in mitochondrial function was unexpected. The correlation coefficients of the COPS5 and top ranked genes such as NDUFB6 and MRPS36 (Fig. 3) were remarkably high with statistical significance. Moreover, not only a handful of genes involved in mitochondria but also the total pathway analysis revealed that mitochondria-related cellular functions and pathways were predominantly over-represented in the COPS5 correlation coefficient table. Although there have not been reports that show direct function of COP9 for mitochondria, correlation analyses with the TP53 of which functions have been extensively studied linked it with known cellular functions including cell cycle and DNA damage response pathways (Additional file 8: Table S15). The result validates the approach in this study, and support the prediction that the COP9 signalosome has a novel functional link to mitochondrial activities.
Few studies have linked COP9 to mitochondria-related cellular activities. Jab1/COPS5 was previously identified as an apoptosis inducing factor through modulation of mitochondrial BH3-domain containing protein BclGs . A study found that in Neurospora crassa, three subunits of the COP9 were required to recover from deficiency of the alternative oxidase function, which is important for mitochondrial functions . However, no follow up studies have been carried out for mammalian cells. On the other hand, pathways involved in ribosomal functions were associated with the COPS5 expression in the tumor tissues. The significance of this link is currently unclear. However, the COP9 structure is homologous to translation initiation factor 3 (eIF3), another high-ordered complex in the cells , and our speculation is that COPS5 may have an overlapping function in the protein synthesis in transformed cells.
The quantity and quality of the genomics database will only be improved. The present study underscores the importance and usefulness of examining the genomics of normal tissues. The results would provide novel functional links of a gene, which was not possible by analyzing the data set of pathophysiological specimens. This approach can reasonably be applied to other genes to identify pathways that have not been known due to scarce information for genomics information of normal tissues.
The present study investigated COPS5 gene expression using the TCGA database. The analysis revealed the highly synergistic mRNA expression pattern among COP9 genes in normal oral tissues. The coordination was abrogated in the tumor tissues, and the regulatory mechanism was taken over by a cis-acting gene regulation. Further analysis revealed an unexpected expression correlation between the COPS5 gene and a number of mitochondria-related genes, postulating the possible functional role of the COP9 signalosome for mitochondrial activities. The quality of the RNA expression data available in TCGA is high enough to allow us gene expression analysis in normal tissues to identify gene functions under normal physiology, which may not be recognized by the data obtained from tumor tissues.
Dubiel D, Rockel B, Naumann M, Dubiel W. Diversity of COP9 signalosome structures and functional consequences. FEBS Lett. 2015;589(19 Pt A):2507–13.
Lingaraju GM, Bunker RD, Cavadini S, Hess D, Hassiepen U, Renatus M, Fischer ES, Thoma NH. Crystal structure of the human COP9 signalosome. Nature. 2014;512(7513):161–5.
Chamovitz DA, Wei N, Osterlund MT, von Arnim AG, Staub JM, Matsui M, Deng XW. The COP9 complex, a novel multisubunit nuclear regulator involved in light control of a plant developmental switch. Cell. 1996;86(1):115–21.
Emberley ED, Mosadeghi R, Deshaies RJ. Deconjugation of Nedd8 from Cul1 is directly regulated by Skp1-F-box and substrate, and the COP9 signalosome inhibits deneddylated SCF by a noncatalytic mechanism. J Biol Chem. 2012;287(35):29679–89.
Shackleford TJ, Claret FX. JAB1/CSN5: a new player in cell cycle control and cancer. Cell Div. 2010;5:26.
Schweitzer K, Bozko PM, Dubiel W, Naumann M. CSN controls NF-kappaB by deubiquitinylation of IkappaBalpha. EMBO J. 2007;26(6):1532–41.
Huang X, Langelotz C, Hetfeld-Pechoc BK, Schwenk W, Dubiel W. The COP9 signalosome mediates beta-catenin degradation by deneddylation and blocks adenomatous polyposis coli destruction via USP15. J Mol Biol. 2009;391(4):691–702.
Bae MK, Ahn MY, Jeong JW, Bae MH, Lee YM, Bae SK, Park JW, Kim KR, Kim KW. Jab1 interacts directly with HIF-1alpha and regulates its stability. J Biol Chem. 2002;277(1):9–12.
Tomoda K, Kato JY, Tatsumi E, Takahashi T, Matsuo Y, Yoneda-Kato N. The Jab1/COP9 signalosome subcomplex is a downstream mediator of Bcr-Abl kinase activity and facilitates cell-cycle progression. Blood. 2005;105(2):775–83.
Lee MH, Zhao R, Phan L, Yeung SC. Roles of COP9 signalosome in cancer. Cell Cycle. 2011;10(18):3057–66.
Bech-Otschir D, Kraft R, Huang X, Henklein P, Kapelari B, Pollmann C, Dubiel W. COP9 signalosome-specific phosphorylation targets p53 to degradation by the ubiquitin system. EMBO J. 2001;20(7):1630–9.
Wei N, Deng XW. The COP9 signalosome. Annu Rev Cell Dev Biol. 2003;19:261–86.
Kato JY, Yoneda-Kato N. Mammalian COP9 signalosome. Genes Cells. 2009;14(11):1209–25.
Kotiguda GG, Weinberg D, Dessau M, Salvi C, Serino G, Chamovitz DA, Hirsch JA. The organization of a CSN5-containing subcomplex of the COP9 signalosome. J Biol Chem. 2012;287(50):42031–41.
Sharon M, Mao H, Boeri Erba E, Stephens E, Zheng N, Robinson CV. Symmetrical modularity of the COP9 signalosome complex suggests its multifunctionality. Structure. 2009;17(1):31–40.
Christmann M, Schmaler T, Gordon C, Huang X, Bayram O, Schinke J, Stumpf S, Dubiel W, Braus GH. Control of multicellular development by the physically interacting deneddylases DEN1/DenA and COP9 signalosome. PLoS Genet. 2013;9(2):e1003275.
Yoshida A, Yoneda-Kato N, Panattoni M, Pardi R, Kato JY. CSN5/Jab1 controls multiple events in the mammalian cell cycle. FEBS Lett. 2010;584(22):4545–52.
Schutz AK, Hennes T, Jumpertz S, Fuchs S, Bernhagen J. Role of CSN5/JAB1 in Wnt/beta-catenin activation in colorectal cancer cells. FEBS Lett. 2012;586(11):1645–51.
Meir M, Galanty Y, Kashani L, Blank M, Khosravi R, Fernandez-Avila MJ, Cruz-Garcia A, Star A, Shochot L, Thomas Y, et al. The COP9 signalosome is vital for timely repair of DNA double-strand breaks. Nucleic Acids Res. 2015;43(9):4517–30.
Leppert U, Henke W, Huang X, Muller JM, Dubiel W. Post-transcriptional fine-tuning of COP9 signalosome subunit biosynthesis is regulated by the c-Myc/Lin28B/let-7 pathway. J Mol Biol. 2011;409(5):710–21.
Deng XW, Dubiel W, Wei N, Hofmann K, Mundt K, Colicelli J, Kato J, Naumann M, Segal D, Seeger M, et al. Unified nomenclature for the COP9 signalosome and its subunits: an essential regulator of development. Trends Genet. 2000;16(5):202–3.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.
Tomoda K, Yoneda-Kato N, Fukumoto A, Yamanaka S, Kato JY. Multiple functions of Jab1 are required for early embryonic development and growth potential in mice. J Biol Chem. 2004;279(41):43013–8.
Cancer Genome Atlas N. Comprehensive genomic characterization of head and neck squamous cell carcinomas. Nature. 2015;517(7536):576–82.
Mundt KE, Liu C, Carr AM. Deletion mutants in COP9/signalosome subunits in fission yeast Schizosaccharomyces pombe display distinct phenotypes. Mol Biol Cell. 2002;13(2):493–502.
Spain BH, Bowdish KS, Pacal AR, Staub SF, Koo D, Chang CY, Xie W, Colicelli J. Two human cDNAs, including a homolog of Arabidopsis FUS6 (COP11), suppress G-protein- and mitogen-activated protein kinase-mediated signal transduction in yeast and mammalian cells. Mol Cell Biol. 1996;16(12):6698–706.
Liu X, Pan Z, Zhang L, Sun Q, Wan J, Tian C, Xing G, Yang J, Liu X, Jiang J, et al. JAB1 accelerates mitochondrial apoptosis by interaction with proapoptotic BclGs. Cell Signal. 2008;20(1):230–40.
Nargang FE, Adames K, Rub C, Cheung S, Easton N, Nargang CE, Chae MS. Identification of genes required for alternative oxidase production in the Neurospora crassa gene knockout library. G3 (Bethesda). 2012;2(11):1345–56.
des Georges A, Dhote V, Kuhn L, Hellen CU, Pestova TV, Frank J, Hashem Y. Structure of mammalian eIF3 in the context of the 43S preinitiation complex. Nature. 2015;525(7570):491–5.
This work was supported by the NIH (NIEHS) Training grant T32ES007266 and NCI grants CA098664, Cancer Center Supporting Grant P30CA17758, and CTSA UL1TR000117.
Availability of data and materials
All the data discussed in this study are available as in main figures or in the supplemental data.
TI conceived the project idea and wrote necessary programs. TI and CAW performed the data analysis. TI and CAW wrote and edited the manuscript, and approved the final manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
(Additional files Figures and Tables) for Analysis of RNA expression of normal and cancer tissues reveal high correlation of COP9 gene expression with respiratory chain complex components by Christina A. Wicker and Tadahide Izumi.
Expressions of COP9 genes in normal oral tissues. The RNA expression levels (y-axis) of COP9 genes from 44 normal oral tissues are shown. Horizontal lines in the boxes: median values. Top and bottom whiskers denote the first and the third quartile (Q1 and Q3). Outliers (> Q3 + 1.5 * (Q3 - Q1) and < Q1 - 1.5 * (Q3 - Q1)) are shown as individual dots. Figure S2. Pair-wise plot for expressions of COP9 genes in the matched lung tissues of normal and squamous cell carcinoma. The plots were generated as described in Fig. 1, except for the data set of the normal (A) and their matched lung squamous cell carcinoma (B) from TCGA was used. The number of samples = 50. (PDF 1508 kb)
The let-7 expression profilesof the matched head and neck normal and tumor tissues. Correlation coefficients for individual let-7 isoforms relative to COP9 genes are shown. Asterisks: p < 0.05. Table S2. To support the result of Table 3, 2 additional normal head and neck tissues were included to calculate correlation of expressions among the COP9 genes. Correlation coefficients of COP9 genes versus COPS5 (rCOPS5) were determined from randomly selected 44 tumor specimens from the entire pool of 511. The calculations were repeated for 10,000 times to obtain their mean and standard deviation (sd) values, which were used to compare the difference from rCOPS5 of normal tissues. Table S3. RNAseq data from young (30 specimens) and elderly patients (31 specimens) were used to calculate the correlation coefficients among the COP9 genes. Table S4. 14 normal oral cavity tissues were pooled based on the TCGA clinical data, and the correlation coefficients and their p values for COPS5 were calculated as described in the main article. Table S5. Data based on the 13 normal tongue tissues were pooled based on the clinical data in the TCGA database, and the correlation coefficients and their p values for COPS5 were calculated. Table S6. 11 normal tissues from smokers were pooled based on the clinical data in the TCGA database, and the correlation coefficients and their p values for COPS5 were calculated. Table S7. Nh: number of genes found within 1 Mbp of each COP9 gene; mean, sd, p and q values were determined as described in Materials and Methods. Table S8. ES for COP9 genes for all KEGG pathways. Table S9. ES for COP9 genes for the oxidative phosphorylation pathway. (XLS 80 kb)
The COPS5 correlation coefficient table list includes 1000 highest |r| values relative to the expression of COPS5 in the normal and tumor tissues of head and neck. (XLSX 621 kb)
GO and KEGG pathways that are over-represented in the normal and tumor oral tissues based on the COPS5 correlation coefficient table. (XLSX 147 kb)
GO and KEGG pathways over-represented by COPS5 expression correlation in the two major anatomical sites in the head and neck (oral cavity and tongue). (XLSX 57 kb)
GO and KEGG pathways over-represented in the normal and tumor lung tissues based on the COPS5 correlation coefficient table. (XLSX 44 kb)
Genes under the chromosomal cis-effect in the HNSCC tissues (XLS 2117 kb)
KEGG pathway analysis for the TP53 gene (XLS 70 kb)