Diverse gene expression patterns in response to anticancer drugs between human and mouse cell lines revealed by a comparative transcriptomic analysis.

The aim of the present study was to perform comparative genomics using gene expression profile datasets of mice and humans who had been treated with anticancer drugs, to determine the similarities and differences in the antitumor mechanisms in the two mammals. This involved data mining of antitumor gene expression regulation, and screening of genetic loci from experimental mouse models of antitumor targets, to provide a theoretical basis of drug design. Subsequently, 9 overlapping genes with opposite expression patterns were identified across mouse and human cell lines that were treated with a specific cyclin‑dependent kinase 4/6 inhibitor, PD0332991. These genes included LIM homeobox 2, adenomedullin, bone marrow stromal cell antigen 1, caveolin 1, histone cluster 1 (HIST1) H2B family member C, HIST1 H3 family member F, low density lipoprotein‑receptor related protein 11, prolyl 4‑hydroxylase subunit α1 and torsin family 3 member A. In addition, the janus kinase‑signal transducer and activator of transcription signaling pathway, Toll‑like receptor signaling pathway, T cell receptor signaling pathway and the nucleotide‑binding oligomerization domain‑like receptor signaling pathway were identified as candidate pathways for explaining antitumor mechanisms.

Introduction

Disorders of the cell cycle are closely associated with the occurrence and development of tumors and are regulated by complex regulatory networks, particularly those involving cyclin-dependent kinases (CDKs), including CDK1, CDK2, CDK3, CDK4, CDK5, CDK6 and CDK7 (1,2). Of these, CDK4/6 has been previously reported to be expressed in various types of cancer, including lung, oral and breast cancer, and hepatocellular carcinoma (3–6). The role of CDK4/6 in tumor cell proliferation is more important than other CDKs and has been established to be one of the essential signaling molecules involved in cell differentiation and metastasis (7). Aberrant CDK activity may lead to uncontrolled proliferation in many tumors, which suggests that inhibiting the activity of CDKs may have a therapeutic benefit. The CDK4/6 inhibitor, PD0332991, has been demonstrated to inhibit tumor cell cycle replication and proliferation in mice and rat animal models in addition to human cancer cell lines including MCF7 (8–11). It has become the prospect for the development of new antitumor drugs. Although PD0332991 is particularly effective in specific antitumor models that induced Rb-dependent cytostasis, identifying the underlying mechanistic pathways of CDK inhibitors is required to fully understand therapeutic responses and this has been investigated in a panel of breast cancer cell lines (12).

Microarray technology is widely used to investigate the biological mechanisms in response to treatment with PD0332991 in humans and mice (8–11). In the present study, a comparative genomics approach was used with relevant gene expression profile datasets of mice and humans that had undergone anticancer drug treatments. The similarities and differences in the antitumor mechanisms in the two mammals at the transcriptomic level were determined. This was achieved by data mining the process of antitumor gene expression regulation, and by screening of genetic loci in a mouse model of an antitumor target, to provide the theoretical basis of drug design.

Materials and methods

Microarray data collection and preprocessing

The microarray dataset of the anticancer efficiency of PD0332991 was obtained from the National Centre for Biotechnology Information Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo). Datasets were reanalyzed if they met the following conditions: i) The data was genome-wide; ii) the dataset contained human and mouse samples; iii) the number of cases and controls in each dataset was ≥3; and iv) complete raw microarray data were available. Based on these criteria, the GSE40514 dataset [contributed by Choi et al (13)] was selected for reanalysis. This superseries is composed of two subseries including GSE40512 and GSE40513. In the dataset of GSE40512, human T-ALL cell line KOPTK1 cells were cultured in the presence of the CDK4/6 inhibitor PD 0332991 with 1 mM or vehicle (VO) for 48 h and the experiment was performed in triplicate. A total of 6 RNA samples (3 vehicle treated and 3 PD0332991 treated samples) were tested by using the platform of Affymetrix Human Genome U133 Plus 2.0 Array (HG-U133_Plus_2, GPL570). In the dataset of GSE40513, mouse breast cancer V720 cells were cultured in the presence of the CDK4/6 inhibitor PD 0332991 with 1 mM or vehicle (VO) for 24 h and the experiment was performed in triplicate. A total of 6 RNA samples (3 vehicle treated and 3 PD 0332991 treated samples) were tested by using the platform of Affymetrix Mouse Genome 430 2.0 Array (Mouse430_2, GPL1261).

R version 3.2.0 (www.r-project.org) and Bioconductor version 3.5 (www.bioconductor.org) were used for data preprocessing (14). The Robust Multichip Average (RMA) algorithm was used in the oligo package to normalize the raw expression data and to generate normalized gene expression intensity in human and mouse cell lines (15). Gene annotation of the probes of human and mouse cell lines was performed using custom written Python code version 2.7 (www.python.org). Probes with no gene annotation or ones that matched multiple gene symbols were removed. Then, the average expression value of multiple probe identities that matched to an official gene symbol was calculated, and this value represented the expression intensity of the corresponding gene symbol. Finally, 20,307 human gene symbols and 20,968 mouse gene symbols were identified. There were 13,976 homologous genes in both human and mouse cell lines. These 13,976 genes were selected for further analysis.

Differential expression gene analysis

Differential expression gene analysis was performed using R version 3.2.2 and the Bioconductor library. The empirical Bayes algorithm (eBayes) in the limma package was used to detect differentially expressed genes between PD0332991-treated and vehicle-treated samples in human and mouse cell lines (16). Fold-change (FC) was calculated as the mean gene expression value in PD0332991-treated samples divided by the mean gene expression value in the vehicle-treated samples. Upregulated genes were considered as a logarithmic transformed fold-changes (log2(FC))≥log2(1.5) and a false discovery rate (FDR) was the adjusted P-value of ≤0.05. Downregulated genes were considered as log2(FC)≤-log2(1.5) and a FDR P-value of ≤0.05.

Gene set enrichment analysis

Java Gene set Enrichment Analysis (GSEA) version 2.2.2 (http://software.broadinstitute.org/gsea/index.jsp) was used for human and mouse samples. The curated Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway software version 5.1 (http://www.genome.jp/kegg) was chosen as the gene sets database to perform the enrichment analysis. Gene sets with <15 or >500 genes were excluded. Phenotype label was set as breast cancer vs. control. The t-statistic mean of the genes was calculated for each KEGG pathway using a permutation test 1,000 times. The upregulated pathways were considered as a normalized enrichment score (NES) >0 and the downregulated pathways were considered as a NES <0. Pathways with a FDR P-value of ≤0.05 were considered to be significantly enriched. A Venn diagram was produced using InteractiVenn (www.interactivenn.net) to demonstrate the enriched KEGG pathways in the human and mouse groups (17).

Results and Discussion

Differentially expressed genes in human T-ALL cells and mouse breast cancer cells

In total, there were 13,976 homologous genes expressed in human and mouse gene expression profiles. Volcano plots are presented in Fig. 1. A greater number of dysregulated genes were identified in the human group compared with the mouse group. A total of 585 upregulated genes and 715 downregulated genes of PD0332991-treated samples were identified in human KOPTK1 cells. By comparison, 125 upregulated genes and 284 downregulated genes were identified in mouse V720 cells. There were 233 common downregulated genes in both human and mouse samples, which was presented in the Venn diagram in Fig. 2. However, there were only 18 upregulated genes that overlapped in these two groups (Fig. 2). The overlapping genes which were significantly up or downregulated across human and mouse cells are presented in Tables I and II.

Figure 1.

Volcano plots of the gene expression profiles of (A) human T-lymphoblast acute lymphoblastic leukemia cells and (B) mouse breast cancer cells. The red points represent upregulated genes and the green points represent downregulated genes. The vertical dotted grey lines represents the log(FC) cutoff and the horizontal grey line represents the logarithmic transformed FDR P-value cutoff. FC, fold-change; FDR, false discovery rate.

Figure 2.

Differentially expressed genes in human T-lymphoblast acute lymphoblastic leukemia cells and mouse breast cancer cells. (A) The log (fold-change) of the dysregulated genes in human and mouse cells. Red lines represent the upregulated genes and the green lines represent the downregulated genes. All genes were sorted by name. Venn diagram of (B) upregulated genes in human and mouse cells, (C) upregulated genes in human cells and downregulated genes in mouse cells, (D) downregulated genes in human cells and the upregulated genes in mouse cells and (E) downregulated genes in human and mouse cells.

Table I.

Enriched Kyoto Encyclopedia of Genes and Genomes pathway analysis in human T-lymphoblast acute lymphoblastic leukemia cells.

A, Upregulated
Pathway	NES	FDR
Allograft rejection	2.22	<0.001
Cytokine-cytokine receptor interaction	2.18	<0.001
JAK-STAT signaling pathway	2.10	<0.001
Toll-like receptor signaling pathway	2.06	<0.001
Type I diabetes mellitus	2.04	0.001
Lysosome	2.00	0.001
Autoimmune thyroid disease	2.00	0.001
Intestinal immune network for IgA production	1.96	0.001
Hematopoietic cell lineage	1.92	0.003
T cell receptor signaling pathway	1.85	0.007
Asthma	1.82	0.009
Natural killer cell mediated cytotoxicity	1.78	0.014
O-glycan biosynthesis	1.69	0.032
Leishmania infection	1.70	0.033
SNARE interactions in vesicular transport	1.65	0.050
B, Downregulated
Pathway	NES	FDR
Cell cycle	−2.49	<0.001
DNA replication	−2.46	<0.001
Oocyte meiosis	−2.26	<0.001
Mismatch repair	−2.14	<0.001
Homologous recombination	−2.12	<0.001
Pyrimidine metabolism	−1.98	0.001
Base excision repair	−1.98	0.001
Nucleotide excision repair	−1.94	0.001
Progesterone-mediated oocyte maturation	−1.88	0.003
One carbon pool by folate	−1.84	0.005
Glycolysis/gluconeogenesis	−1.68	0.034

JAK-STAT, janus kinase-signal transducer and activator of transcription; SNARE, soluble N-ethylmaleimide-sensitive factor activating protein receptor; FDR, false discovery rate; NES, normalized enrichment score.

Table II.

Enriched Kyoto Encyclopedia of Genes and Genomes pathway analysis in mouse breast cancer cells.

A, Upregulated
Pathway	NES	FDR
Cytosolic DNA-sensing pathway	2.27	<0.001
NOD-like receptor signaling pathway	2.17	<0.001
Hematopoietic cell lineage	1.98	0.003
Focal adhesion	1.99	0.004
Toll-like receptor signaling pathway	1.83	0.018
Leishmania infection	1.85	0.019
Prion diseases	1.84	0.020
Cytokine-cytokine receptor interaction	1.79	0.024
JAK-STAT signaling pathway	1.77	0.026
ECM-receptor interaction	1.77	0.028
Glycosaminoglycan degradation	1.73	0.037
Type I diabetes mellitus	1.72	0.040
Complement and coagulation cascades	1.70	0.045
Lysosome	1.69	0.046
B, Downregulated
Pathway	NES	FDR
Cell cycle	−2.71	<0.001
DNA replication	−2.67	<0.001
Homologous recombination	−2.47	<0.001
Oocyte meiosis	−2.41	<0.001
Spliceosome	−2.37	<0.001
Nucleotide excision repair	−2.32	<0.001
Mismatch repair	−2.27	<0.001
Base excision repair	−2.20	<0.001
Progesterone-mediated oocyte maturation	−2.11	<0.001
RNA degradation	−2.10	<0.001
Pyrimidine metabolism	−1.99	0.001
RNA polymerase	−1.71	0.019
Cysteine and methionine metabolism	−1.62	0.043

JAK-STAT, janus kinase-signal transducer and activator of transcription; NOD, nucleotide-binding oligomerization domain; ECM, extracellular matrix; FDR, false discovery rate; NES, normalized enrichment score.

To investigate the anticancer efficiency of PD0332991, 9 genes with diverse regulatory patterns between mouse and human were identified (Fig. 3). In the PD0332991 group, LIM homeobox 2 (LHX2) was overexpressed in humans compared with the vehicle group, whereas in mice, LHX2 expression was reduced (Fig. 3). Dysregulated expression of LHX2 in human cancers has been previously reported (18). In addition, cancer-associated epigenetic alterations were significantly associated with LHX2 regulation, including DNA methylation and microRNA regulation (19–22). In the present study, adrenomedullin (ADM), bone marrow stromal cell antigen 1 (BST1), caveolin 1 (CAV1), histone cluster 1 H2B family member C (HIST1H2BC), histone cluster 1 H3 family member F (HIST1H3F), low density lipoprotein-receptor related protein 11 (LRP11), prolyl 4-hydroxylase subunit α1 (P4HA1) and torsin family 3 member A (TOR3A) were downregulated in human KOPTK1 cells, and were upregulated in mouse V720 cells. HIST1H3F demonstrated the largest difference between PD0332991- and vehicle-treated samples in human and mouse cell lines (log2(FC)=−1.77, FDR=5.99E-6 in human; log2(FC)=1.21, FDR=8.20E-4 in mice) from the 8 genes investigated. The histone HIST1H3F has been previously identified as a biomarker for predicting cancer risk (23). In addition, the mean log2 gene expression value of HIST1H2BC in humans was ~4. However, the mean value in mice was ~12. This suggested that there was a large difference between human and mouse cell lines. HIST1H2BC is one of the most significant susceptible locus in the transcriptome study of schizophrenia, which has been also reported to be a suitable reference gene in non-small cell lung cancer (24,25). The association between ADM and epithelial-mesenchymal transition, and the underlying signaling pathways in intrahepatic cholangiocellular carcinoma has been previously elucidated (26). The overexpression of ADM in tumor cells may accelerate breast cancer bone metastasis, which suggests that it may be a potential target for therapeutic intervention against bone metastases (27). The association between movement disorders and genetic polymorphisms of BST1 has been confirmed in different populations, including Parkinson's and Alzheimer's disease (28–30). BST1 has been also investigated in a previous study on the drug resistance of tumors, including anticancer drugs for artesunate resistance (31). CAV1 may have a tumor suppression role in breast cancer and its expression is regulated by CpG island shore methylation (32). LRP11 is a member of the low-density lipoprotein receptor family, and has a potential role in mediating cellular drug uptake (33). One of the isoforms of P4HA, P4HA1, has an important role in prostate cancer progression and is associated with microRNA-124 (34). TOR3A functions with the gene encoding the Ras-related protein RAB1B, which regulates molecular transport, protein trafficking and developmental disorders (35).

Figure 3.

Expression profiles of nine target genes in human T-lymphoblast acute lymphoblastic leukemia cells and mouse breast cancer cells. The dark color bars represent the gene expression value in the PD0332991-treatment group and the light color bars represent the gene expression values in vehicle-treatment group. Error bars represent standard deviation of the mean expression values of each gene of three replicates. LHX2, LIM homeobox 2; ADM, adenomedullin; BST1, bone marrow stromal cell antigen 1; CAV1, caveolin 1; HIST1H2BC, histone cluster 1 H2B family member C; HIST1H3F, histone cluster 1 H3 family member F; LRP11, low density lipoprotein-receptor related protein 11; P4HA1, prolyl 4-hydroxylase subunit α1; TOR3A, torsin family 3 member A.

Anticancer efficiency of PD0332991 on multiple pathways levels

GSEA analysis was performed to identify the signaling pathways that were significantly associated with anticancer effect of PD0332991. The findings of the present study revealed that there were 15 upregulated pathways and 11 downregulated pathways in PD0332991-treated human KOPTK1 cells. In addition, there were 14 upregulated pathways and 13 downregulated pathways in PD0332991-treated mouse V720 cells (Fig. 4A and B). There were 7 common upregulated pathways and 9 common downregulated pathways in the human and mouse samples. However, oppositely dysregulated pathways in the two groups were not identified. The details of these pathways are presented in Tables I and II. The downregulated pathways in human and mouse cell lines were primarily associated with cell cycle, DNA replication and nucleotide metabolism. The upregulated pathways in the two groups treated with PD0332991 were involved in various intracellular signaling pathways, including the JAK-STAT, Toll-like receptor, T cell receptor and nucleotide-binding oligomerization domain (NOD)-like receptor signaling pathways.

Figure 4.

The comparison of enriched Kyoto Encyclopedia of Genes and Genomes pathways analysis in response to PD0332991 treatment between human and mouse. (A) Upregulated pathways in human T-lymphoblast acute lymphoblastic leukemia cells (red) and mouse breast cancer cells (orange). (B) Downregulated pathways in human (blue) and mouse (cyan) cells.

It has been previously reported that the JAK/STAT pathway may be aberrantly activated in solid tumors including breast cancer, whose activation in malignant and nonmalignant cells contributes to the cancer pathogenesis and therapeutic response (36,37). Toll-like receptor (TLR)3 signaling is an integral component of solid tumors, and inhibition of this pathway may increase the effectiveness of current treatments of breast cancer in combination with anticancer drugs (38). The TLR8 signaling pathway has been associated with the functional regulation of tumor-specific ∆γ T cells, which are important contributors to innate immunity against cancer (39). Gene polymorphisms of NOD1/caspase recruitment domain (CARD)4 and NOD2/CARD15 may be associated with altered risk of diverse malignancies, including breast cancer (40).

In conclusion, the present study identified 9 overlapping genes with opposite expression patterns across PD0332991-treated human T-ALL cells and mouse breast cancer cells. In addition, various cellular signaling pathways were identified to be closely associated with the efficacy of anticancer drugs. Therefore, this data may aid future breast cancer drug screening and design.