Genomic Characterization of Esophageal Squamous Cell Carcinoma from a High-Risk Population in China

Genomic instability plays an important role in most human cancers (1–3). Several questions arise regarding genomic instability in esophageal squamous cell carcinomas (ESCC): (a) How prevalent is genomic instability in these tumors? (b) Is there a relation between different types of genome-wide instability in ESCC, such as loss of heterozygosity (LOH) and copy number (CN) loss/gain? (c) What is the association between genomic instability and risk factors and clinical phenotypes? The ability to answer some of these questions may lead to a better understanding of tumorigenesis and the development of new strategies for prevention, early detection, and therapy.

A combined analysis of changes in both DNA and RNA from tumors is a useful approach to identifying DNA alterations that are important for tumor development, and genome-wide genomic instability and gene expression have been simultaneously evaluated for several different cancer types (4, 5). Such integrated analysis for ESCC has been examined only once previously, and then only in a study of cell lines (6). Identification of DNA alterations in early stage tumors can be used for cancer diagnosis as well as etiology and prevention. High-throughput identification of genetic alterations that affect gene expression remains a challenging task. Recently, several reports focused on the relation between DNA variants and gene expression in the human genome and tumors using several methods, most commonly comparative genomic hybridization (7–16). Thus far, analyses based on the identification of genetic loci have not been resolved at the level of individual genes and polymorphic alleles that affect gene expression (17).

ESCC is a common malignancy worldwide and one of the most common cancers in the Chinese population; Shanxi Province in north central China has some of the highest esophageal cancer rates in the world (18, 19). Previously, we identified several regions of LOH and CN alteration in ESCC using microsatellite markers and low-density single nucleotide polymorphism (SNP) arrays (20–24). Here, we analyzed DNA from 30 microdissected ESCC tumors and compared them to germ-line DNA from the same case using the Affymetrix 500K SNP array. First, we examined relations among three types of genome-wide instability—LOH, CN loss, and CN gain—and then examined associations of these genomic instability measures to ESCC risk factors, clinical characteristics, and prognosis. Second, we identified regions with particularly high LOH and DNA CN change that may contain specific tumor suppressor genes or oncogenes related to ESCC tumorigenesis or disease progression. Third, we compared findings from the current study with our two previously published studies of genome-wide ESCC LOH in the same population. Finally, we compared individual gene CN alteration and mRNA expression.

Case selection. This study was approved by the Institutional Review Boards of the Shanxi Cancer Hospital and the US National Cancer Institute (NCI). Cases diagnosed with ESCC between 1998 and 2001 in the Shanxi Cancer Hospital in Taiyuan, Shanxi Province, PR China, and considered candidates for curative surgical resection were identified and recruited to participate in this study. None of the cases had prior therapy and Shanxi was the ancestral home for all. After obtaining informed consent, cases were interviewed to obtain information on demographics, cancer risk factors (smoking, alcohol drinking, and detailed family history of cancer), and clinical information. Cases were followed for survival status to the end of 2003. The cases evaluated here were part of a larger case-control study of upper gastrointestinal cancers conducted in Shanxi Province (25).

Biological specimen collection and processing. Venous blood (10 mL) was taken from each case before surgery and germ-line DNA from whole blood was extracted and purified using the standard phenol/chloroform method.

Tumor and adjacent normal tissues were dissected at the time of surgery and stored in liquid nitrogen until use. One 5-μm section was H&E stained and reviewed by a pathologist from NCI to guide the microdissection. Five to 10 consecutive 8-μm sections were cut from fresh frozen tumor tissues. Tumor cells were manually microdissected under light microscopy. DNA was extracted from microdissected tumor as previously described (26) using the protocol from Puregene DNA Purification Tissue kit (Gentra Systems, Inc.). RNA from tumor and matched normal tissue was extracted using the protocol from PureLink Micro-to-midi total RNA purification system (Invitrogen). RNA quality and quantity were determined using the RNA 6000 Labchip/Aligent 2100 Bioanalyzer (Agilent Technologies).

Target preparation for GeneChip Human Mapping 500K array set. The Affymetrix GeneChip Human Mapping 500K array set contains ∼262,000 (Nsp I array) and ∼238,000 (Sty I array) SNPs (mean probe spacing, 5.8 Kb; mean heterozygosity, 27%) according to the manufacturers protocol.

Experiments were conducted according to the protocol (GeneChip Mapping Assay manual) supplied by Affymetrix, Inc. Briefly, DNA samples were diluted to ∼50 ng/μL in reduced EDTA TE Buffer (0.1 mmol/L EDTA) and assayed according to the GeneChip Mapping Assay manual. A total of 250 ng of DNA was digested with Nsp I or Sty I for 120 min at 37°C, and the reaction was inactivated at 65°C for 20 min. The digested DNA was then ligated to Nsp I or Sty I adaptors for 180 min at 16°C, followed by 20 min at 70°C before subsequent PCR amplification. All aforementioned steps were carried out in the pre-PCR clean room. The PCR protocol consisted of the following: 94°C for 3 min, followed by 30 cycles of 94°C for 30 s, 60°C for 30 s, and 68°C for 15 s, with a final extension at 68°C for 7 min. PCR was done with the DNA Engine Tetrad PTC-225 (MJ Research). After PCR, a mixture of 3 μL of PCR product and 3 μL of the 2× Gel Loading Dye was electrophoresed on a 2% Tris-borate EDTA gel at 120 V for 30 min to assess successful amplification. If the expected product sizes (200–2,000 bp) were observed, purification and elution of the PCR products were performed using Qiagen MiniElute 96 (Qiagen), followed by DNA quantification using spectrophotomeric analysis. Samples were diluted to a final concentration of 90 μg in 45 μL volume for fragmentation at 37°C for 35 min, followed by 95°C for 15 min. Fragmentation was verified by performing electrophoresis of 4 μL of fragmented DNA from each sample in a 4% Tris-borate EDTA gel at 120 V for 30 min. Successful fragmentation was confirmed by the presence of a smear between 50 to 200 bp. The samples were end-labeled with biotin and hybridized onto the array. The chip was incubated at 49°C for 18 h in the Affymetrix Genechip system hybridization oven, then washed and stained in the Genechip Fluidics Station 450 (Affymetrix) following the manufacturer's instructions. The chip was scanned with the Affymetrix GeneChip Scanner 3000 using GeneChipOperating System 1.4, and the data files were automatically generated. Genotype calls were generated by GTYPE v 4.0 software (Affymetrix). Germ-line and tumor DNA from each case were run together in parallel in the same experiment (i.e., same batch, same day). The Gene Expression Omnibus accession number for these array data is GSE15526.

Probe preparation and hybridization for Human Genome U133A 2.0 array. The Human U133A 2.0 array is a single array representing 14,500 well-characterized human genes (Affymetrix). The array experiment was performed using 1 to 5 μg total RNA; reverse transcription, labeling, and hybridization followed the protocol provided by the manufacturer (Affymetrix).

Genechip 500K array data analysis. Probe intensity data from the Affymetrix 500K SNP array was used to identify autosomal alterations in the present study.

LOH was defined in a traditional manner as a change in genotyping call from heterozygous (AB) in the germ-line DNA to homozygous (AA or BB) in the matched microdissected tumor DNA (all genotype calls generated by using GTYPE, Affymetrix).

CN loss or gain was based on a comparison of tumor with germ-line DNA. Microarray data were first normalized using the gtype-probeset-genotype package included in Affymetrix Power Tools version 1.85. Each tumor sample was individually normalized via the BRLMM algorithm along with 99 blood samples. These blood samples were obtained from the 30 ESCC cases evaluated in the present study plus 69 healthy controls (age-, sex-, and region-matched to cases) who were all part of a larger case-control study of upper gastrointestinal cancers conducted in Shanxi Province. Paired CN analysis was then performed on each tumor sample using the Affymetrix Copy Number Analysis Tool. DNA obtained from the case's blood served as the normal control; a window of 100 kb was chosen to optimize the identification of extended regions of CN alteration. The output of the Copy Number Analysis Tool program is CN state rather than an absolute CN prediction: normal CN corresponds to a state of 2; zero and 1 correspond to CN loss; and states 3 and 4 correspond to CN gain. Therefore, we treated CN loss or gain as a qualitative trait.

Case genomic instability. Each case was categorized as high or low status for each of the three genomic instability measures (LOH, CN loss, CN gain) based on his/her frequency across the entire genome using a cutoff of 10% (≥10%, high; <10%, low), which was approximately a median split for each measure (Table 1A). An overall genomic instability score was calculated for each case as follows: a case was coded “0” for no high (i.e., ≥10%) genomic alterations (LOH, CN loss, or CN gain), “1” for one high measure, “2” for two high measures, and “3” for all three measures high. Cases with 0 genomic instability score were called the low-frequency genomic instability group; cases with scores 1, 2, or 3 were the high-frequency genomic instability group.

Chromosomal arm high genomic instability. A chromosomal arm was considered to show high genomic instability (i.e., high LOH or high CN loss or high CN gain) when the average instability frequency among all cases combined for all SNPs (informative SNPs only for LOH, all SNPs for CN loss/gain) across the entire chromosomal arm was ≥10% for the specific genomic instability measure of interest. To maximize the likelihood that we would identify only true, “nonrandom” areas of change, for a chromosomal arm to be called high genomic instability, we also required that the instability measure be present in ≥50% of SNPs in at least 20% (n = 6) of cases.

Subchromosomal region LOH. To identify focal regions of LOH, we also evaluated SNPs in cytobands. For consideration, a region had to have a minimum size of 15 informative SNPs. We declared an LOH region only when all informative markers showed LOH in at least 10% (n = 3) of cases, with LOH in a minimum one-third of informative SNPs in at least 15% (n = 4) of cases.

Subchromosomal region CN loss/gain. We also identified focal regions of CN loss/gain in cytobands. We declared a CN loss (gain) region when, within a sliding windows of 30 consecutive SNPs (regardless of genotype call), all 30 markers showed CN loss (gain) in at least 30% (n = 9) of cases, and at least one-sixth of the markers showed CN loss (gain) in at least 50% (n = 15) of cases.

Human Genome U133A 2.0 array data analysis. The Affymetrix GeneChip Human U133A 2.0 array is a single array with 22,000 probe sets representing 14,500 well-characterized human genes. Robust Multiarray Average algorithm (27, 28) implemented in Bioconductor in R was used for background correction and normalization across all samples. We applied paired t tests to each of the 22,000 probe sets to identify genes differentially expressed between tumor and matched normal samples. To account for multiple comparisons, we selected genes that showed significant differences with P values of <0.05 after Bonferroni adjustment.

Association between genomic instability and risk factor/clinical data and survival. All statistical analyses were performed using Statistical Analysis Systems (SAS Corp) for assessment of the relations between genomic instability and risk factors (smoking, family history), clinical characteristics (stage, grade, metastasis), and survival. LOH and CN alterations were modeled in several different ways. For summary, genomic instability measures over the entire genome (22 autosomes) outcomes were evaluated as high versus low (defined as median split) for LOH, CN loss, and CN gain, as well as a combined average score, all in relation to individual risk and clinical factors. LOH, CN loss, CN gain, and the combined average score were also examined as continuous variables in relation to risk/clinical factors using general linear models. Additionally, we analyzed relations of LOH, CN loss, and CN gain to risk/clinical factors separately for each chromosomal arm. For the chromosomal arm-specific analyses, we divided LOH and CN alterations into high or “nonrandom” frequency (≥50% of SNPs affected) versus low or “random” frequency (<50% of SNPs affected). Overall survival was examined by LOH and CN alteration status (high versus low frequency categorized as ≥50% versus <50% of SNPs affected) per chromosomal arm with Kaplan-Meier curves; differences were tested using the log-rank test. Pearson correlation analyses were used for comparisons between CN alterations and gene expression levels in tumors. All P values were two-sided and considered statistically significant if P value was <0.05.

Risk factors and clinical characteristics for cases are shown in Table 1A. The average age of cases was 54 years (range, 39–67 years), females predominated (19 of 30), approximately one-third smoked (nearly all males), and about one-third had a positive family history of upper gastrointestinal cancer. Tumors were most commonly stage 3 (24 of 30) and grade 2 (24 of 30), one-half the cases had metastatic disease at the time of surgery, and the median survival from the time of surgery was 666 days.

The overall average genotype call rate was 96% (89–99%) based on a total of 126 SNP array chips, including three cases whose blood DNAs were repeated on both the Nsp I and Sty I SNP arrays for quality control purposes. The average call rate for the 250K Nsp I array was 96% (90–98%) and for the 250K Sty I array was 96% (89–99%). The genotype call rates on microdissected tumor DNA (95% for Sty I and 96% for Nsp I) and germ-line DNA (96% for Sty I and 99% for Nsp I) were similar for both chips. The average present call rate on the Human Genome U133A array was 53% (range, 51–61%) for the 34 chips from the 17 sample pairs with sufficient tissue for RNA isolation and testing.

Case genomic instability. Genome-wide LOH, CN loss, and CN gain in the 30 ESCC cases studied here are shown in Table 1A. Using a frequency of ≥10% as a cutoff for high-frequency instability for each of these three measures, one-half of the cases showed high LOH, 11 cases had high CN loss, and 19 cases had high CN gain (Table 1A and B). Based on our criteria for categorizing cases into high or low genomic instability groups (see Materials and Methods), 11 cases had no high genomic instability measure, 3 cases had one high measure, 7 cases had two high measures, and 9 cases had all three high measures. Altogether, 11 cases had low-frequency and 19 cases high-frequency genomic instability.

LOH and CN loss/gain were analyzed in relation to the case risk factors and clinical characteristics. None of the risk factors or clinical characteristics we examined showed a significant association with any of the three genomic instability measures, with one exception: CN gain on chromosome 3q was positively associated with metastasis (nominal P = 0.025).

Overall and chromosomal arm genome-wide LOH. The overall LOH frequency for all 30 ESCC cases across all chromosomal arms combined was 10.5% (median; range of <1–39%; Table 1A, cases ordered by increasing LOH frequency), whereas the LOH frequency for the 39 individual (autosomal) chromosomal arms (Chr 1-22) ranged from <1% to 33% (Table 1C). Nine chromosomal arms showed high LOH frequency, including 3p, 4p/q, 9p/q, 13q, 17p/q, and 21q, based on the criteria described (Table 1C). LOH frequency overall was highest on chromosomal arm 9p (33%), where 40% of cases had LOH in at least 50% of informative SNPs (Supplementary Table S1A, cases ordered by ID).

Overall and chromosomal arm genome-wide CN loss/gain. Overall genome-wide CN loss and gain in all 30 ESCC cases are shown in Table 1A. CN loss (median, 5.5%; range, <1–24%) was less frequent than CN gain (median, 11.5%; range, <1–20%). Six chromosomal arms had high CN loss, including 3p, 4p/q, 8p, 9p, and 11q (Table 1C). Supplementary Table S1B shows the frequency of CN loss on these six chromosomal arms for each case; the highest CN loss occurred on 3p where 33% of cases had CN loss in at least 50% of SNPs.

Six chromosomal arms were identified with high CN gain: 3q, 7p, 8q, 14q, and 20p/q (Table 1C). The number of SNPs that showed CN gain varied widely. Many cases exhibited CN gain so extensive as to essentially encompass the entire chromosomal arm. CN gain on at least 50% of SNPs was seen in 43% of cases for 3q and 40% of cases for 8q, whereas almost half the cases showed no CN gain at all on the remaining four chromosomal arms (i.e., 7p, 14q, 20p/q) with high CN gain overall (Supplementary Table S1C).

Chromosomal arm comparisons between LOH and CN loss/gain. We noted both high LOH and high CN loss on four chromosomal arms (i.e., 3p, 4p/q, 9p), suggesting that LOH there was caused by chromosome loss. However, some chromosomal arms exhibited high LOH without high CN loss (i.e., 9q, 13q, 17p/q, 21q), indicative of mitotic recombination or a nondisjunction event, and are target regions of interest for CN neutral LOH exploration. In contrast, chromosomal arms with high CN gain were distinct, and did not overlap with chromosomal arms that exhibited either high LOH or high CN loss (Fig. 1).

Subchromosomal region genomic instability. A total of 22 LOH regions were identified based on our criteria (see Materials and Methods), including 11 on chromosome 9, 7 on chromosome 17, and 4 on chromosome 13. These regions represent 395 SNPs from 21 genes (Table 2A). Three CN loss regions were found—3p12.3, 4p15.1, and 9p21.3—which contain 765 SNPs from seven genes (ie, CTNT3, LOC33897, LOC28529, LOC645716, MTAP, CDKN2A, and CDKN2B; Table 2B). Twelve CN gain regions were recognized as follows: six on chromosome 3q, one on 7q, four on 8q, and one on 11q; these 12 regions represented 13,343 SNPs from 482 genes (Table 2C). Among the 37 regions of genomic instability found, a single cytoband, 9p21.3, which contains CDKN2A and CDKN2B, showed both LOH and CN loss (Table 2A and B).

Comparison of genome-wide LOH with the two previous studies. Genome-wide LOH results from this study are compared with two previous studies of ESCC from the same population in China in Table 3 (21, 23). Of note, the cases studied in the three studies described here were all different and nonoverlapping. Our first genome-wide LOH study was performed on 11 ESCC cases and used 366 microsatellite makers with an average heterozygosity of 76%. High LOH frequencies (≥50%) were observed on 14 chromosomal arms (3p, 4p/q, 5q, 8p/q, 9p/q, 11p/q, 13q, 17p/q, 18p; ref. 19). A second study performed on 26 ESCC cases used the Affymetrix 10K SNP array and found high frequency of LOH (≥50%) on 10 chromosomal arms (3p, 4p/q, 5q, 9p/q, 13q, 15q, 17p/q). The present study showed high LOH (≥10%) on 9 chromosomal arms (3p, 4p/q, 9p/q, 13q, 17p/q). Table 3 shows that LOH frequency decreased as the number of markers increased, suggesting that the frequency of LOH in previous studies was overestimated by the low density marker platforms, and that the presence of LOH in ESCC tumors was not as high as previously thought. Alternative explanations for this finding include that the high density array increased the probability of markers disrupting the requirement for continguous LOH regions, or that differences in marker polymorphic information content or marker placement may have played a role. Table 3 also shows that eight chromosomal arms (3p, 4p/q, 9p/q, 13q, 17p/q) consistently showed the highest LOH frequencies across all three studies, suggesting that these regions likely contain tumor suppressor genes involved in ESCC development and disease progression.

Relation between CN alterations and mRNA expression. To investigate the relationship between CN alterations and mRNA expression level, we focused our analysis on chromosome 11q13.1-q13.4, one of the CN gain regions, to evaluate the association between gene expression and CN gain. 11q13 is a highly gene-rich region strongly conserved across zebrafish, mice, and humans, and it has exhibited multiple amplification peaks in several tumors (29). The high CN gain region 11q13.1-q13.4 identified in the present study contained 777 SNPs from 157 genes. We found that 37 of the 157 genes (203 SNPs) had CN gain in >30% of cases (range, 33–63%). We measured RNA expression in 17 cases with sufficient frozen material available for RNA isolation from both tumor and matched normal tissues, and determined that 26 of the 37 genes represented on the U133A chip probesets were suitable for analysis. Among these 26 genes, there were strong (significant) positive correlations between CN gain and expression for 21 of 26 genes, with Pearson correlation coefficients between 0.51 to 0.87 (P < 0.05), including PSCA1, CCND1, CTTN, PPFIA1, and SHANK2 (Table 4; Fig. 2). This result suggests that high CN gain is associated with up-regulated gene expression. To identify specific genes in the focal CN gain regions, we will need to characterize both gene expression and other cellular activities in a larger population of cases in the future.

Our study characterized ESCC tumors for three types of genome-wide instability—LOH, CN loss, and CN gain—in germ-line DNA and matched microdissected tumor DNA using the Affymetrix GeneChip Human Mapping 500K array, and is the first study to report the integration of high density LOH and CN alteration data in ESCC with gene expression analyses on a genome-wide scale.

We observed that chromosomal arms with high LOH and CN loss frequently overlapped (i.e., 3p, 4p/q, 9p), suggesting that these areas potentially harbor tumor suppressor genes, as illustrated in a previous study that showed high LOH on 9p concurrent with frequent CDKN2A mutations and intragenic allelic losses (30). The present study had analogous findings: subchromosomal region 9p21.3 exhibited both high LOH (Table 2A) and high CN loss (Table 2B), whereas one-half the ESCC cases studied had CN loss on all five SNPs in CDKN2A, and nine cases showed biallelic loss for CDKN2A (data not shown). These data indicate that numerous types of alterations occur in CDKN2A—LOH, CN loss, biallelic loss, germ-line and somatic mutations, and intragenic allelic loss—and that these alterations collectively contribute to the inactivation of CDKN2A.

Using criteria stricter than traditionally applied, we identified 22 focal subchromosomal regions of LOH. These regions were found on 9p/q, 13q, and 17p/q, similar to a previous LOH study in which we used the Affymetrix 10K SNP array (24). Both studies used microdissected tumor and matched germ-line DNA, which we believe are essential approaches to obtaining consistent LOH results. The concordance of findings between these two LOH studies using similarly rigorous methods increases our confidence in the importance of emphasizing these 30 regions in future studies to identify ESCC tumor suppressor genes.

Six chromosomal arms showed high CN loss and six arms showed high CN gain, but over all 22 autosomal chromosomes, CN gain occurred twice as often as CN loss. There was no overlap in the chromosomal arms which showed high CN gain compared with arms which showed high CN loss, nor did high CN gain arms overlap with high LOH arms. This is consistent with observations that chromosome-wide CN gain is limited to one or two copies (trisomy or tetrasomy), whereas high-level amplification tends to be a focal event (e.g., gene amplification at 11q13 encompassing CCND1). Although there have been several previous reports of CN alterations in ESCC, most have used low-resolution comparative genomic hybridization methods (6, 24, 31–36). Two have used SNP-based array methods (24, 36), comparable with the present study, but these were also low-resolution arrays.

We previously reported that CDC25B overexpressed both mRNA (37, 38) and protein in ESCC tumors and dysplasias (38). CDC25B is located on 20p13, a chromosomal arm with high CN gain in the present study, where we also observed that 30% of ESCC cases showed CN gain in CDC25B (data not shown), suggesting that CN gain is responsible for CDC25B overexpression in at least a subset of persons with premalignant and invasive ESCC.

Among the 12 CN gain regions identified, six were on 3q, four on 8q, one on 7q22.1, and one on 11q13.1. We focused on 11q13.1, the most gene-rich of these gain regions. Several genes on chromosomal region 11q13 that were amplified in this study are known oncogenes whose amplification has been associated with poor prognosis (39). For example, PPFIA1 is a member of a family of leukocyte common antigen-related trans-membrane tyrosine phosphatase-interacting proteins. Recently, Tan and colleagues (29) reported that PPFIA1 was overexpressed in head and neck squamous cell carcinoma, and suggested that this gene may act as an invasion inhibitor in head and neck squamous cell carcinoma. PPFIA1 was highly (7-fold) overexpressed in our data, and expression was highly correlated with CN (Table 4).

In summary, our study showed that genomic instability varied widely among ESCCs and included cases with both high and low frequencies; the high-frequency instability group may harbor germ-line variants or acquired somatic mutations in genes that maintain genomic stability. Genome-wide studies from this high-risk population show a consistent pattern of high LOH on selected chromosome arms, which are targets in searching for loss-of-function genes involved in ESCC. Our findings also show the potential utility of combining CN and gene expression data to identify genes involved in esophageal carcinogenesis. Future studies should combine results in tumors with germ-line genotypes to find functional changes, and determine if these changes are associated with genetic susceptibility to ESCC or might serve as early detection markers.

No potential conflicts of interest were disclosed.

Grant support: Intramural Research Program of the NIH, the National Cancer Institute, the Division of Cancer Epidemiology and Genetics, and the Center for Cancer Research.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

We thank Dr. Stephen Hewitt of NCI for reviewing the slides used for microdissection.