The KM12 model cell lines, established by intrasplenic or subcutaneously injecting poorly metastatic KM12C cells derived from a Dukes' B stage colon carcinoma into nude mice to isolate variant lines KM12L4A and KM12SM, both with high liver metastatic potential (15), were obtained from Dr. Isaiah J. Fidler's laboratory at University of Texas MD Anderson Cancer Center, and cultured according to their published protocols (15). The highly metastatic variants are genetically related to the poorly metastatic parental cell line KM12C and gain of 20q has been identified in all the three lines (16). The SW480, SW620 and the HCT116 cells were obtained from the American Type Culture Collection (ATCC) and cultured according to the recommended protocols. The SW620 cells were derived from a metastatic lesion in the same patient whose primary tumor was the source of the SW480 cells unlike the KM cell line model in which the metastatic variants were isolated from mouse xenografts of primary tumor derived cells. The cell line SW480 was cultured from a Duke's Stage B colon carcinoma, while the cell line SW620 was from an abdominal lymph node in the same patient (17). Additionally, the HCT 116 cell line was derived from a poorly differentiated human colon carcinoma (18). Genomic DNA was isolated from these cells using the Qiagen DNeasy Blood and Tissue Kit (cat. no. 69504; Qiagen) according to the manufacturer's protocols and stored at −20°C. Total cellular RNA was isolated from the cells using TRIZOL reagent (cat. no. 15596026; Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol and further purified using the Qiagen RNeasy RNA Extraction Kit with an on-column DNase digestion (cat. no. 74104; Qiagen).

Gene copy number analysis and reverse transcription-quantitative (RT-qPCR) were performed using DNA and RNA of 23 matched tumor-normal tissue samples that were snap-frozen or stored in RNAlater (cat. no. AM7020; Invitrogen; Thermo Fisher Scientific, Inc.), obtained from the laboratories of Drs. Stanley Hamilton and Marsha Frazier at UT MD Anderson Cancer Center. The five polyp samples with matched WBC and normal tissue were obtained from Dr. Stephen Meltzer's laboratory at the Johns Hopkins University School of Medicine. All samples including normal and tumor tissues were obtained from residual samples remaining after clinical diagnosis under IRB approved protocols with ethics approval and participant consent received at the time of admission approved by the Institutional Review Boards of Johns Hopkins University School of Medicine and UT MD Anderson Cancer Center (approval no. LAB90-018).

DNA from KM12, SW, HCT116 cell lines, five tumors and four polyps that passed the QC were used for the arrayCGH analysis. Commercially available normal female genomic DNA (Promega Corporation; containing DNA from various anonymous donors) was used as the source of control DNA for the experiments. DNA from the samples was diluted to a concentration of 25 ng/ml and 1 ml was used in a Phi29 amplification step using the Qiagen Repli-g amplification kit (cat. no. 150023; Qiagen) according to Agilent-modified manufacturer's protocol, omitting the initial denaturing step. DNA labeling, hybridization, and washing were carried out according to manufacturer's protocol. Genomic arrays were scanned with the Molecular Devices GenePix 4000B Microarray scanner and fluorescence intensities, normalization and ratios obtained using the Agilent Feature Extraction Software 7.5.1 at the recommended settings. Dye-swap hybridizations were performed for each cell line, tumor and polyp samples, and the resulting data combined and further analyzed by statistical methods as mentioned below. Heat maps were generated using Agilent's CGH Analytics software Version 3.4 and Nexus Copy Number 7 software.

Total RNA from cell lines and commercially available normal (non-diseased) colonic epithelial RNA (Stratagene; Agilent) was labeled, hybridized and washed using manufacturer's protocol and reagents [Agilent Low RNA Input Labeling Kit (cat. no. 5184-3523), Agilent in situ Hybridization Kit (cat. no. 5190-0404), Agilent's Stabilization and Drying Solution cat. no. 5185-5979)]. Arrays were scanned on Axon 4000B scanner and Agilent DNA Microarray scanner. Dye-swap hybridizations were performed for each cell line and the resulting data combined. The initial fluorescence data was generated with GenePix Pro 5.1 software (Axon; Molecular Devices, LLC) for all images produced with the Axon scanner, and normalization, filtering and ratios calculated as described below.

Total RNA was used to produce cDNA with the Invitrogen First Strand cDNA synthesis kit according to manufacturer's protocol (cat no. 18091050). RT-qPCR was carried out on isolated total RNA in MD Anderson Cancer Center's DNA Core facility. RT-qPCR analysis was performed using ABI TaqMan Assays for the candidate genes, showing gain of expression, as well as for the internal control gene ACTB, with an ABI 7500HT Real-Time PCR System, according to manufacturer's recommendations. All assays were run in duplicate and the relative expression of the selective genes in terms of fold change in the cancer cell lines, tumor and polyp samples, compared with normal tissue samples was calculated using the comparative CT (2^−ΔΔCq) method (19), using the endogenous control to normalize the data. Primers used are summarized in Table SI.

MSI status was evaluated by fluorescence labeled microsatellite marker PCR followed by capillary electrophoresis fragment size analysis using an ABI 3130 sequencer and Genescan software (Applied Biosystems; Thermo Fisher Scientific, Inc.). Five markers (BAT25, BAT26, D2S123, D5S346, D17S250; of the National Cancer Institute panel) were analyzed. The samples were classified as MSI-High (MSI-H) if two or more markers showed altered allelic size, MSI-Low (MSI-L) if one marker showed allelic shift and MSS if none showed allelic shift (20).

For the independent MDACC test dataset, the paired neoplastic and non-neoplastic tissue controls were analyzed for MSI according to the standard-of-care immunohistochemistry test method (21). Immunoperoxidase stains were performed on formalin-fixed paraffin embedded sections with antibodies for the DNA mismatch repair enzymes MLH1, MSH2, MSH6 and PMS2. The samples were classified as MSI-H if two or more markers were altered, and MSS if none of the markers were altered. The samples unevaluated for MSI were filtered out for the analysis.

DNA was extracted from frozen tissue blocks by QIAamp Mini columns according to manufacturer's directions (cat. no. 56304; Qiagen) and was analyzed by the Sequenom MALDI TOF MassArray system and software (Sequenom Inc, San) for the following mutations [PIK3CA (E542A, E542K, E545A, E545K, H1047R/L), AKT1 (E17K), BRAF (V600E), KRAS (284 G12X, 285 G12C, G13D)].

The Agilent Feature Extraction software was used to identify chromosomal gains and losses, with significance being considered at P-value less than 0.01 (22). The ploidy information on the cell lines obtained from spectral karyotyping (SKY) was used to correct for the software's assumption of overall diploidy in both test and reference channels.

For expression array analysis, the Axon generated fluorescence intensity for test and reference channels were treated as follows: median background was subtracted from median signal for each probe and the 75th percentile over probes in each channel was standardized to 1,000 with a multiplicative constant. The resulting expressions were then truncated at a lower minimum value of 25. These values were then base 2 log transformed. A linear model including factors for experimental and biological variation was fit to each probe to account for differential expression by cell line model (Patient), primary vs. metastatic strains nested within patient and date blocks as a random factor: y _g,i,m(i),j,k=µ_g + Patient_g,i + Metastasis_g,m(i) + Date_g,j + ε_g,i,m(i),j,k (1); Date_g,j ~ N(0, σ²_D,g) (2); ε_{g,i,m(i),j,k(i,j)} ~ N(0, σ²_e,g) (3); g=1 to 41,675 i=1, 2 m (1)=1 to 3 and m (2)=1, 2 (4); j=1,2 replicates k=1 to n k(i,j) (5).

Expression analysis was performed based on the above models generated with variable factors, such as, biological variables of the samples, date variables when experiments were performed, g refers to signal differentials, i refers to dye swaps, j refers to replicates.

For the analysis using TCGA sample cohort, a total of 195 patients with complete tumor profiles and clinical data were used as the data set (23). Copy Number Variation and Expression array data (RNAseq) for these patients were downloaded from the website cBioPortal.org (14). All patients who showed an expression z-score >=1.96 for the candidate genes were considered to have a gain of function altered gene status.

The second dataset from MDACC has been profiled using Agilent Microarray, which were all hybridized against a common normal colon reference pool, according to manufacturer's protocol. Lowess based normalization was performed using the Agilent Feature extraction software with background subtraction. Samples that had both expression and mutation profiling data (n=182) were only considered for the subsequent analysis. Expression data corresponding to the four candidate genes of interest was extracted for further analysis. The gene BMP7 has three probes in the dataset. The probe A_24_P91566 was filtered due to poor correlation with the other two probes-A_23_P154643 and A_23_P68487 (Pearson correlation ~0.09). The Mean expression of the remaining two probes is taken as the representative expression of BMP7 for the downstream analysis.

Chi-square or Fisher's exact test were computed to test the differences in the distribution of categorical variables between patients with genetic alteration in any of the four genes in the signature and those without any alteration as the reference group. The categorical variables included Expression Subtypes, aggressiveness (defined as having lymph node spread and/or distant metastasis), MSI status and the BRAF mutation status. To assess the association between clinical outcomes and genetic alteration status, the adjusted odds ratios (ORs) with the 95% confidence intervals (95% CI) were estimated using unconditional multivariable logistic regression analyses, with adjustment for the covariates. All statistical tests were two-sided with a statistical significance level of P<0.05. Data were analyzed using Statistical Analysis System/Genetics software program (SAS/STAT version 9.1.3; SAS Institute, Inc.).

Mitotic shake-off preparations for each of the cell lines were fixed at room temperature in Carnoy's fixative (3:1 methanol: acetic acid) for 30 min and slides were prepared by dropping cell suspensions onto wet slides. SKY was carried out on aged slides using the Human SkyPaint Probes kit according to manufacturer's protocol (Applied Spectral Imaging Inc.). Visualization and analysis were carried out with an Olympus AX70 microscope using Image Analysis (Applied Spectral Imaging Inc.) and karyotype images were generated using SKYVIEW software (Applied Spectral Imaging Inc.). Modality was determined based on an average of 5 metaphases captured for each cell line.

Genomic bacterial artificial chromosome (BAC) and P1-derived artificial chromosome (PAC) clones with inserts ranging in size of 150–200 kb representing the genomic intervals of minimum common region of amplification (MCR) on chromosome 20q and residing near or containing the candidate genes were obtained from the Sanger Institute or generously provided by Dr. Craig Chinault at Baylor College of Medicine. These BAC and PAC clones were cultured and their DNA isolated using the Eppendorf miniprep kit. The DNA was then labeled in a polymerase reaction with digoxigenin-dUTP (cat. no. 11573152910; Millipore Sigma). Human Cot-1 DNA (cat. no. 15-279-011; Invitrogen; Thermo Fisher Scientific, Inc.) was added, and the DNA was precipitated in ethanol, followed by 70% ethanol wash. The DNA pellet was briefly dried in a speed-vac then re-suspended in Hybrisol (cat. no. S4040; EMD Millipore). Probe was added to the slides prepared as stated above and hybridized under coverslip in a humidified chamber for 24 h. The coverslip was then removed, the slides were washed and the probe was detected with anti-digoxigenin-rhodamine FAB (cat. no. 11207750910; Millipore Sigma) and counterstained with DAPI (1 µg/ml). Slides were examined and signals counted in 5–10 interphase cells.

To validate accurate copy number estimates from the intensity normalized aCGH data, the modal number for each cell line was determined by SKY (Fig. S1). In addition, the status of chromosome 20 in the karyotypes of the cells utilized was compared with previously published karyotypes. Notable genomic changes in all the cell lines are listed in Table SII. Most rearranged chromosomes were observed previously (16). Gain of chromosome 13 was identified in KM12C, which had not been previously reported. KM12C (near diploid), KM12L4A (pseudo-triploid) as well as KM12SM (near tetraploid) karyotypes contained the translocation of chromosome 20 to 22. In the primary tumor-derived SW480 cells (pseudo-diploid) as well as the metastasis-derived SW620 (near-diploid), translocation between chromosomes 5 and 20 were consistent with the previous observations (24). Ploidy correction was employed to aid in the normalization of the fluorescence intensity ratios obtained from the array GCH experiments analyzed by the Feature Extraction software from Agilent Technologies (Version 7.5.1).

CGH to the Agilent Human Genome CGH 44A Microarray (resolution of ~75 kb) was carried out on the five cell lines, with dye swaps. Human female genomic DNA was used as the reference. The genome-wide profile was generated after circular binary segmentation analysis for each cell line (data not shown). Overall, KM12C, KM12L4A, and SW480 showed a lower level of noise in the data than KM12SM and SW620. Several genomic regions showed amplifications and deletions in more than one cell line, but since our primary aim was to identify candidate genes on the 20q amplicon, our analysis was focused on chromosome arm 20q. The high-resolution aCGH profile of chromosome 20 for each of the five cell lines is represented in a heat map (Fig. 1A), demonstrating copy gains along the q arm in all the five cell lines. After considering the copy gain profile of all the cell lines, three MCRs of amplification across the 20q arm were identified spanning 29–34 megabases (Mb), 40–45 and 54–57 Mb. Copy number verification for these intervals was carried out through interphase FISH using BAC and PAC probes located within these intervals (Fig. 1B-G, Table I). Based on FISH analysis of each probe, the genomic copy number data estimated from the aCGH analysis was confirmed.

The present study was extended to primary tumor samples by performing aCGH analyses on a set of five colorectal tumor and four polyp samples. Differential amplification of varying genomic intervals within the three MCRs was detected in the tumor and the polyp samples (Fig. 2). The overall significant gains observed in polyps and tumors are included in Table SIII. The HCT116 CRC cell line, with no known chromosome 20q gain (25,26), was also investigated as a reference cell line for amplification profile and as expected, did not show any significant gains or losses on chromosome 20q.

To identify the genes with gain of expression localized in the MCRs identified above, genome-wide microarray expression analysis of the five cell lines compared with commercially available normal colon total RNA was carried out, and the expression patterns of the 1061 oligo probes for chromosome 20 were examined. It was found that there was a variety of expression patterns within these intervals, including significantly downregulated genes as well as upregulated and unchanged gene expression profiles. To focus on the genes residing within the MCRs, the data was filtered based on the chromosomal position of the genes represented in the list of 86 probes that displayed overexpression on chromosome 20q in the five cell lines. This reduced the list to 18 probes representing 14 known genes and 2 open reading frames. Further filtering of the list to include the probes showing elevated expression in the metastasis-derived cell lines led to the identification of four genes (BMP7, DNMT3B, UBE2C and YWHAB) as the candidate gene signature for further investigation. The RPL13A transcript was selected as the endogenous control for its consistent expression detected across the cell lines and reference normal colon total RNA.

To validate the expression array data, RT-qPCR for the gene signature consisting of BMP7, DNMT3B, UBE2C and YWHAB was carried out in the cell lines and the training set of twenty-three tumors and five polyps. The expression patterns detected with both methods (expression array and RT-qPCR) were consistently altered in the same direction, with varying amplitude in the cell lines, tumors and polyps. The tumor samples showed elevated expression of the candidate genes compared with the matched normal, though the degree of increase was relatively less compared with the cell lines. A reference cell line HCT116 also revealed increased expression of the candidate genes following RT-qPCR analyses (Fig. 3) suggesting that alteration in the expression levels of these genes do not always result from gain of gene dosage. Furthermore, the tumor samples were analyzed for their MSI status. Of the 20 tumors analyzed, 16 tumors had altered status of the four gene panel and of these, eleven revealed MSS or MSI-L phenotypes.

To validate the statistical significance of the data obtained with the training set in a larger cohort of test set patient samples, copy number and mRNA expression data were analyzed in patients with CRC included in TCGA database (23), accessible through the cBioPortal for Cancer Genomics (14), as well as an independent set of colorectal adenocarcinoma samples that were annotated in-house at MD Anderson Cancer Center. Our primary aim was to investigate whether increased expression of the gene signature was associated with aggressive tumor behavior such as lymph node spread and/or presence of distant metastasis, as well as their degree of association with MSS or MSI. Table II shows the total number of patients from the TCGA dataset with an increased expression of the four candidate genes. Of the 182 patients that had complete tumor profiles, ~56% of the patients showed alteration in one or more of the four genes. A total of 74 patients showed lymph node spread and/or distant metastasis, with a significant number of these patients overexpressing the gene signature (P<0.05, TCGA data as of August 2013) (Table III). Specifically, a 2-fold increased risk for lymph-node spread/distant metastasis was observed in patients with at least one gene showing alteration (defined as ‘altered’ group; adjusted OR=1.97, 95% CI=1.05–3.67) after adjustment for sex, cancer site, family history of cancer and prior diagnosis of cancer. Furthermore, a significant association between the gene signature and absence of MSI in these patients was also observed (Table III), with majority of these patients showing the MSS phenotype, as observed in the 20 tumor samples of the training set analyzed in this study. Overall, the results of the TCGA data analyses corroborated the results obtained from the training set samples. A similar trend of overexpression associating with MSS status was also observed in the MDACC dataset, although statistical significance for this association was evident at lower stringency of z>1.28, 80% CI (P=0.03) (Fig. S2).

The four gene signature showed statistically significant association with the outcome ‘Expression subtypes’ in the TCGA data defined by the gene expression profile correlated with genomic instability. Furthermore, the gene signature showed association with additional CRC characteristics, such as, methylation status (P<0.0001) and stage of cancer (P=0.048). These findings were not surprising since expression subtype, MSI status, lymph node spread/distant metastasis and stage of cancer outcomes were themselves associated with each other (P<0.01, data not shown). Additionally, the percentage of patients with overexpression of gene signature was examined in relation to the status of known CRC ‘driver gene’ mutations. The gene signature expression tended to correlate with mutations in the APC gene (P=0.065; Table II). It was noteworthy that the frequency of patients with overexpression in the gene signature was markedly higher compared with the frequency of patients with alterations in other well-characterized genes associated with CRC such as BRAF, PIK3CA, SMAD4, FBXW7 (23) (data not shown). Interestingly, association of gene signature with absence of BRAF mutation was also observed in both the TCGA (P=0.06, z>1.96, 95% CI) and MDACC (P=0.07, z>1.28, 80% CI) datasets (Fig. S3).

Taken together, these observations suggest that expression of the chromosome 20q gene signature, identified in the present study, represents functionally significant genetic alterations cooperating with the driver gene alterations in the progression of CRC.