An integrative analysis reveals functional targets of GATA6 transcriptional regulation in gastric cancer.

Lineage-restricted transcription factors (TFs) are frequently mutated or overexpressed in cancer and contribute toward malignant behaviors; however, the molecular bases of their oncogenic properties are largely unknown. As TF activities are difficult to inhibit directly with small molecules, the genes and pathways they regulate might represent more tractable targets for drug therapy. We studied GATA6, a TF gene that is frequently amplified or overexpressed in gastric, esophageal and pancreatic adenocarcinomas. GATA6-overexpressing gastric cancer cell lines cluster in gene expression space, separate from non-overexpressing lines. This expression clustering signifies a shared pathogenic group of genes that GATA6 may regulate through direct cis-element binding. We used chromatin immunoprecipitation and sequencing (ChIP-seq) to identify GATA6-bound genes and considered TF occupancy in relation to genes that respond to GATA6 depletion in cell lines and track with GATA6 mRNA (synexpression groups) in primary gastric cancers. Among other cellular functions, GATA6-occupied genes control apoptosis and govern the M-phase of the cell cycle. Depletion of GATA6 reduced the levels of the latter transcripts and arrested cells in G2 and M phases of the cell cycle. Synexpression in human tumor samples identified likely direct transcriptional targets substantially better than consideration only of transcripts that respond to GATA6 loss in cultured cells. Candidate target genes responded to the loss of GATA6 or its homolog GATA4 and even more to the depletion of both proteins. Many GATA6-dependent genes lacked nearby binding sites but several strongly dependent, synexpressed and GATA6-bound genes encode TFs such as MYC, HES1, RARB and CDX2. Thus, many downstream effects occur indirectly through other TFs and GATA6 activity in gastric cancer is partially redundant with GATA4. This integrative analysis of locus occupancy, gene dependency and synexpression provides a functional signature of GATA6-overexpressing gastric cancers, revealing both limits and new therapeutic directions for a challenging and frequently fatal disease.

INTRODUCTION

Some lineage-restricted transcription factors (TFs) specify developing tissues and regulate cell-specific genes in adults. Cancers often amplify such TF genes, including MITF in melanoma (1), NKX2-1 in lung adenocarcinoma (2), SOX2 in squamous esophageal cancer (3), and AR in prostate cancer (4). As tumors may depend on amplified TF genes (2, 3), they are potential targets for cancer therapy. However, TFs other than nuclear hormone receptors are notoriously difficult drug targets (5, 6). Therefore, core downstream genes and pathways might suggest alternative targets that are more sensitive to small molecules. GATA4 and GATA6 are amplified in up to 30% of gastric and esophageal adenocarcinomas (7) and GATA6 depletion in the latter specifically impairs anchorage-independent cell growth (8, 9). We studied the dependencies and transcriptional functions of this TF.

Worldwide, stomach cancer is the second leading cause of cancer death (10, 11). Somatic copy number amplifications (SCNAs) or mutations of ERBB2, EGFR, MET, and FGFR2 offer avenues for targeted therapy in few patients (12–14). Esophageal adenocarcinomas, which are closely related, frequently amplify GATA6 and GATA4 (7), TF gene loci that show especially high expression in gastric and duodenal epithelia (15, 16). In mouse intestine, GATA6 levels are highest in the crypts, where cell proliferation is reduced in conditional Gata6 mice (17). TF co-occupancy, determined by chromatin immunoprecipitation (ChIP), further suggests that GATA6 mediates crypt functions together with CDX2, a master intestinal regulator (18). Because gastric cancer frequently arises in a background of intestinal metaplasia (19), this partnership suggests that GATA gene amplifications may promote proliferative, crypt progenitor-like properties in stomach epithelial cells.

GATA6 is also amplified in pancreas cancer (20, 21), but interference with its functions is hampered by limited information about the targets of transcriptional control. To delineate core downstream genes, pathways, and functions in gastric cancer, we examined genome-wide GATA6 occupancy in relation to GATA6-dependent gene expression in cell lines and GATA6-associated gene expression (synexpression) in human tumor samples. This approach revealed features, consequences and core transcriptional targets of GATA6 in gastric cancer.

RESULTS

Amplification and expression of GATA genes in upper digestive tract cancers

Small regions on chromosomes 8p and 18q are focally amplified in 17% to 22% of stomach and gastro-esophageal junction (G-EJ) adenocarcinomas (7). GISTIC analysis (22) of these cases and public SCNA data from 321 additional primary stomach cancers identified GATA6 and GATA4 as the only genes within the minimal common areas of amplification (Fig. 1A). Among hundreds of diverse cancers, high-level GATA4 amplifications were largely confined to gastric cancer and GATA6 amplifications to stomach and pancreas adenocarcinomas (Suppl. Fig. S1A). GATA4 and GATA6 are homologous TFs that recognize the same DNA sequence and have overlapping functions in some mouse tissues (23, 24), suggesting that they may serve similar roles in gastric cancer. Because GATA6 amplifications are more common and GATA6 antibodies (Ab) perform well in tissue and chromatin studies, we concentrated on this TF.

Figure 1

Somatic copy number alterations (SCNAs) in adenocarcinomas of the upper digestive tract

A) Summary of SCNA data from SNP array analysis in studies to date: Dulak et al. (7), GSE31168 in the Gene Expression Omnibus, and stomach adenocarcinoma STAD in The Cancer Genome Atlas (TCGA). Results of GISTIC analysis are represented by chromosome and show that the minimal common regions of amplification on chromosomes 8p and 18q encompass the single gene loci GATA4 and GATA6, respectively. Also see Suppl. Fig. S1A. B) Representative immunohistochemical analysis of GATA6 and CDX2 in human tissue, showing abundance of both proteins in a gastric cancer with GATA6 amplification. GATA6 is expressed in normal stomach epithelium and often increased in cancer, whereas CDX2 is expressed only in intestinal metaplasia and cancer; additional examples appear in Suppl. Fig. S1B. The right panels display GATA6 and CDX2 expression in CRC. Original magnification 10X, inserts show 6X additional digital magnification. C) Collated GATA gene expression from RNA-seq data in the TCGA_stomach adenocarcinoma database, showing GATA6 overexpression in a significant fraction of cases and that, although GATA4 mRNA loss is common, primary tumors rarely lose GATA6 expression. Color scales interpret the degree of GATA factor mRNA overexpression.

GATA6 is expressed in normal human stomach epithelium, intestinal metaplasia, and carcinoma (Fig. 1B, Suppl. Fig. S1B). Both GATA6-amplified and unamplified gastric cancers showed strong immunostaining in most tumor cells (Suppl. Fig. S1C), consistently stronger than in colorectal cancer (CRC, Fig. 1B, Suppl. Fig. S1B), where GATA6 amplification is uncommon (Suppl. Fig. S1A). Many gastric cancers also express the intestine-restricted factor CDX2, consistent with their likely origin in areas of intestinal metaplasia (19), but GATA6 levels showed no association with CDX2 expression or tumor cell differentiation (Suppl. Fig. S1C). mRNA analysis of 290 gastric cancers in The Cancer Genome Atlas (TCGA, https://tcga-data.nci.nih.gov/tcga/) showed high GATA6 expression in a significant fraction of cases and GATA4 overexpression in fewer cases (Fig. 1C). Tumors with high GATA6 amplification showed significantly higher mRNA levels than diploid samples and GATA6 mRNA or protein (Suppl. Fig. S1C) were rarely lost. In summary, GATA6 amplification is common in gastric cancers and elevates expression, although tumors lacking amplification also may overexpress GATA6.

Expression of GATA6 and possible downstream genes in gastric cancer cell lines

Focal GATA gene amplifications were uncommon in gastric cancer cell lines (Suppl. Figs. S1A, S2A). Only 2 of 26 lines, HUG1N and GC1Y, showed GATA6, and none showed GATA4, amplification (Suppl. Fig. S2B). Protein expression among 9 gastric cancer lines was highest in HUG1N, at levels similar to those found in the CRC cell line Caco-2 (Fig. 2A). Immunoblotting detected no GATA6 in GC1Y cells and two lines lacking amplification, AGS and SNU16, expressed almost the same GATA6 levels as HUG1N. Thus, gene amplification does not invariably produce excess protein and is only one basis for overexpression. Indeed, despite the paucity of SCNAs at chromosome 18p, gastric cancer cell lines commonly express abundant GATA6 mRNA, and hierarchical clustering of gene expression similarities (Euclidean distance, 1000 most variant probe sets) across 36 lines revealed GATA-overexpressing lines as a distinct sub-class (Fig. 2B). The 15 cell lines in this group included all GATA6-expressing lines we detected by immunoblotting: AGS, HUG1N, NCI-N87, and SNU-16. GATA4 mRNA was overexpressed in fewer lines, many of which showed high GATA6. The presence of a characteristic mRNA profile in GATA6-overexpressing gastric cancers suggests that some of these genes represent its transcriptional targets. We chose HUG1N and AGS cells to study GATA6 functions because these lines express more GATA6 than GATA4 (Fig. 2A) and HUG1N, in particular, carries multiple gene copies.

Figure 2

GATA factor expression and synexpression in gastric cancer cell lines

A) Immunoblot analysis of GATA6 and GATA4 in Caco2 CRC cells and 9 gastric cancer lines chosen for presence (HUG1N, GC1Y) or absence of GATA6 gene amplification. B) Hierarchical clustering of mRNA expression in 36 gastric cancer cell lines, based on the 1000 probe sets that vary the most across this cell line collection (variance calculated by Euclidean distance). Bars with a red-blue scale represent relative expression levels of GATA6 (vertical bar) and GATA4 (horizontal bar) mRNAs in each cell line. This analysis revealed 2 distinct tumor classes. One class encompasses GATA6 (and/or GATA4) mRNA-overexpressing cell lines (bottom left), including all GATA6-expressing lines from the immunoblot analysis (A, which are highlighted by asterisks). Gene expression in these 17 cell lines clusters separately from all others, revealing that GATA factors associate specifically with a sizable body of genes and putative transcriptional targets in gastric cancer. The data also indicate that GATA6 is often (but not always) highly expressed, even without gene amplification.

Effects of GATA6 depletion on cell replication and gene expression

Using lentivirus-delivered shRNAs to deplete GATA6, we confirmed that cells selected for viral integration were efficiently depleted compared to those receiving a control shRNA (Fig. 3A). AGS, SNU16, and Caco2 cells all proliferated slower after GATA6 depletion (Fig. 3A). Proliferation of HUG1N was unaffected by short-term depletion (≤6 cell doublings), suggesting possible adaptation in culture to escape dependency. By contrast, in a pooled-shRNA screen of 89 cancer cell lines (25), GATA6 dependency was strongly correlated with its expression level: over ≥16 cell doublings, HUG1N and other high GATA6-expressing gastric and CRC cell lines were among the most susceptible to GATA6 depletion (box #3 in Fig. 3B). Moreover, transcripts affected by GATA6 depletion in HUG1N and Caco2 cells were significantly correlated with GATA6-synexpressed mRNAs in primary tumors and cell lines (Fig. 3C). Thus, transcripts that respond to GATA6 deficiency closely match in vivo gene expression and are likely enriched for bona fide target genes. Transcripts altered in GATA6-deficient HUG1N cells were correlated significantly with synexpressed genes in CRC and even better with synexpressed genes in gastric cell lines and tumors (Fig. 3C), suggesting that GATA6 may control some genes specifically in the stomach or colon and other genes in both organs.

Figure 3

Loss of GATA6 impairs growth and affects gene expression in gastrointestinal cancer cells

A) GATA6-overexpressing cell lines were infected with lentiviral vectors carrying GATA6-specific or a non-specific (NS) 21-bp shRNA, selected in puromycin for viral integration, and assessed daily in triplicate for cell viability using MTT assays. Immunoblots verified efficient GATA6 depletion in all cell lines for each replicate, using actin as a loading control, and each plotted value represents the mean (±SD) optical density. B) Correlation of GATA6 expression levels with GATA6 dependency in an unbiased, pooled-shRNA screen of >11,000 genes in human cancer cell lines (Ref. 25). 89 cell lines were grouped into 4 bins according to the relative level of GATA6 mRNA, 0=absence, 3=high. Relative response in growth of these lines to GATA6 shRNA is shown in box plots with the centers indicating median sensitivity, edges at the 25th and 75th percentiles, and whiskers extending to the extreme non-outlier data points. Red and blue dots represent HUG1N and AGS, respectively, and many high-expressing, shRNA-sensitive cell lines were from CRCs. C) Concordance of GATA6 expression correlations in primary tumors and cell lines and transcripts differentially expressed in GATA6-deficient cells. Pearson correlations with GATA6 expression were computed from primary tumor (125 stomach, 1,732 CRC) and cell line (36 gastric, 59 CRC) datasets. Correlation coefficients for every gene were then correlated with fold-changes (Δ exp) from shRNA experiments in HUG1N and Caco2 cells. The resulting Pearson correlations are displayed in a heat map, with zeroing of the diagonal. Asterisks denote significant concordance based on a t-test and P <0.001.

Compared to cells treated with a control shRNA, 805 transcripts dropped in levels in GATA6-deficient HUG1N cells and 595 transcripts were increased. Like its homolog GATA1 (26, 27), GATA6 may both activate and repress genes directly or one function could be direct and the other indirect. A non-parametric test of all log2 fold-changes in transcript levels in GATA6-depleted HUG1N cells revealed that mean changes were significantly enriched for reduced expression (one-sided Wilcoxon signed rank test, P=2 × 10−8), suggesting that in unperturbed cells GATA6 activates many more genes than it might repress. In line with this idea, no functional Gene Ontology category was enriched among 595 transcripts that increased in HUG1N cells or 144 transcripts that increased in both Caco2 and HUG1N. In contrast, functional categories were readily apparent among transcripts that declined upon GATA6 loss (Suppl. Fig. S3A and B), providing clues to its cellular functions.

The 212 transcripts reduced in both Caco2 and HUG1N cells depleted of GATA6 represent a highly significant, non-random overlap (Fig. 4A; odds ratio 9.0, p<4 × 10−103 by Fisher exact test) and were enriched for Gene Ontology terms related to the cell cycle, particularly M-phase (Fig. 4B). GATA6 knockdown in Caco2, AGS and SNU16 cells resulted in smaller G1- and S-phase fractions and higher G2/M-phase fractions than control shRNA-treated cells, indicating M-phase dysfunction (Fig. 4C). Phospho-histone H3 immunofluorescence and immunoblot analyses verified this G2/M-phase arrest (Figs. 4D–E). Although the cell cycle seemed overtly undisturbed in HUG1N cells, qRT-PCR analysis showed reduced M-phase related transcripts in GATA6-depleted cell lines, including HUG1N (Fig. 4F). These data implicate late cell cycle control as a core GATA6 function.

Figure 4

Depletion of GATA6 impairs cell growth and cell cycle progression in gastrointestinal cancer cells

A) Venn diagram showing the overlap in transcripts reduced in GATA6-depleted Caco2 and HUG1N cells. B) Table of functional categories enriched among the 212 genes reduced in both cell lines and statistical significance of this enrichment over the genome background. Shaded Gene Ontology terms highlight the cell cycle, particularly M-phase. C) Flow cytometry analysis of the cell cycle in GATA6-depleted AGS, Caco2 and SNU16 cells. Propidium iodide staining of cell cycle phases is shown; each value represents the mean ± SD of 3 replicates. D–E) Immunofluorescence (D, merged images of green PH3 immunostaining and blue DAPI) and immunoblot (E, 3 replicate lanes, Noc = 100 ng/mL nocodazole positive control) analysis of phospho-histone H3. GATA6-depleted AGS cells show G2/M-phase arrest, compared to cells treated with a control shRNA. F) Relative expression in AGS, HUG1N and Caco2 cells of M phase-implicated transcripts (derived from the Gene Ontology classification in B) that showed reduced expression in GATA6-depleted cells. qRT-PCR results are expressed as the average fold-change ± SD relative to cells treated with a non-specific shRNA in 3 independent samples.

Delineation of GATA6 cistromes

To identify primary transcriptional targets, we used ChIP-seq to localize GATA6 binding in HUG1N cells. After adjusting for results from mock IgG ChIP, we mapped nearly 7,000 high-confidence GATA6 binding sites in HUG1N cells and considered the data together with occupancy previously mapped in Caco2 cells (18). High sequence conservation within ChIP fragments and greatest enrichment for the consensus WGATAA motif, present in at least 70% of bound regions, implied direct occupancy at most sites (Fig. 5A). Although GATA6 occupied many regions within 3 kb of transcriptional start sites (TSSs), it predominantly bound DNA far from promoters (Fig. 5A). For example, binding near FGFR2, which is implicated and occasionally amplified in gastric cancer (28, 29), occurs 8 kb and 13 kb from the promoter (Fig. 5B). Moreover, GATA6 sometimes occupied exactly the same sites in HUG1N and Caco2 cells (e.g., ANKRD30BL locus, Fig. 5C) or different sites near the same locus (e.g., FGFR2, Fig. 5B), but often bound sites in only 1 of the 2 cell lines (Fig. 5D). Binding unique to each line did not trivially reflect deletion of the corresponding genomic region in the other. To determine if these differences reflect tissue-specific binding in the stomach and colon, we mapped GATA6 occupancy in overexpressing AGS gastric cancer cells and observed modest overlap among the 3 cell lines (Fig. 5E). Sites confined to 1 or 2 lines seemed not to reflect technical artifacts in the others but bona fide absence of GATA6 at those regions (Figs. 5B,D,E). In aggregate, ChIP-seq revealed diverse GATA6 binding to cellular DNA, with cell line-specific, tissue-specific, and shared binding sites. We reasoned that common binding sites and transcriptional targets might represent the most pertinent outcome of frequent gene amplification.

Figure 5

Genome-wide GATA6 occupancy in gastrointestinal cancer cell lines

A) High sequence conservation (left) and significant enrichment of the canonical WGATAR sequence motif (middle) among GATA6 ChIP sites, whose pie-chart distribution across the genome (right) reveals localization mainly in introns and intergenic regions far from promoters. B–D) Wiggle tracks from GATA6 ChIP-seq in HUG1N and Caco2 cells illustrate diverse binding patterns: same locus, different sites (B, FGFR2); same locus, overlapping sites (C, ANKR030BL); and exclusive occupancy in 1 line (D, KCNJ5). E) Venn diagram representing the overlap of GATA6 occupancy in HUG1N, AGS and Caco2 cells. Representative wiggle traces from the 1007 sites common to the 2 stomach cell lines demonstrate co-occupancy in all 3 lines at the FDPS promoter (top) and exclusive binding in the 2 gastric cancer cell lines in a BTRC intron (bottom).

GATA6-dependent genes, synexpression, and DNA occupancy together impute target genes

To determine roles of GATA6 binding in gene regulation, we measured concordance of fold-changes in gene expression in GATA6-deficient HUG1N cells with its occupancy near affected genes. Fig. 6A plots the Pearson correlations between binding sites and genes affected by GATA6 depletion; negative correlations denote higher concordance. Altered gene expression in HUG1N cells correlated with GATA6 binding in HUG1N cells better than binding in Caco2 cells, and ChIP-seq in AGS gastric cancer cells did not improve this correlation. Fig. 6B shows heat maps of ordered average fold-changes per gene, from most repressed (left, blue) to most activated (right, red), and the corresponding average GATA6 binding near those genes (<20 kb, brighter yellow denotes higher occupancy). GATA6 binding in HUG1N cells alone (Fig. 6B) or in both HUG1N and AGS cells (data not shown) correlated with altered gene expression. To identify direct transcriptional targets, we first considered genes that bind GATA6 within 20 kb of the TSS and also reduce expression in GATA6-depleted cells. Gene Ontology analysis of such genes in HUG1N (Suppl. Fig. S3C) suggests that GATA6 directly regulates 75 genes that control cell replication; 41 of these genes promote cell proliferation, and GATA6 also binds DNA near dependent genes implicated in cell death. Among genes that bind GATA6 and depend on it in Caco2 cells, 86 genes regulate cell proliferation (Suppl. Fig. S3D). “Digestion” and “Response to hormone stimulus”, functions strongly linked to gastrointestinal physiology, were also highly enriched among genes near GATA6-binding sites in HUG1N and Caco2 cells, respectively, suggesting a breadth of cellular roles.

Figure 6

Identification of GATA6 target genes from analysis of DNA binding and effects on gene expression

A) Pearson correlation coefficients of GATA6 binding (ChIP-seq data) and genes affected by GATA6 depletion (microarray analysis following shRNA treatment). GATA6 binding was encoded as 0 or 1; negative correlations denote higher concordance. B) Top: heat map of sorted transcript fold-changes after GATA6 depletion, averaged over bins of 500 genes and ordered from most reduced (left, blue) to most increased (right, red). Bottom: heat map of average GATA6 occupancy within 20 kb of the corresponding 500-gene bins. C) Pearson correlation coefficients of nearby GATA6 binding and the PCA ranks of genes based on reduced expression in GATA6-deficient cells and synexpression in human cancer specimens and cell lines; negative correlations denote higher concordance. D) Top: heat map of combined PCA ranking score, averaged over windows of 500 genes and ordered from best correlated (left, blue, high negative z-scores) to the least correlated (right, red, high positive z-scores). Bottom: heat map of GATA6 occupancy within 20 kb of the corresponding bins of 500 PCA-ranked genes. E) The 20 best-scoring genes by PCA rank, extracted from Suppl. Table S1; genes showing nearby GATA6 binding in at least 2 cell lines are highlighted. F) GATA6 occupancy in cell lines within 20 kb of the TSS of PCA-ranked gene groups. G) Representative wiggle trace showing GATA6 occupancy at the SPINK4 locus, representing a high-ranking gene from panel E.

Reasoning that bona fide transcriptional targets might additionally co-express with GATA6 RNA in human cancer tissues, we integrated these cell line data with information on synexpressed genes, i.e., correlation scores for every gene’s expression with GATA6 expression. We used Principal Component Analysis (PCA) (30) to determine directions of maximal variation in the combined space of altered gene expression in GATA6-deficient HUG1N and Caco2 cells (axes 1 and 2) and GATA6 synexpression across 4 large datasets of colon and gastric cancer tissues and cell lines (axes 3–6); the first 2 PCA axes accounted for 53% of the variation. A linear combination of normalized PCA scores for each gene, incorporating cell line and tumor sample data, yielded a ranked list of candidate GATA6 transcriptional targets (Suppl. Table 1). This integrated PCA ranking (Fig. 6C) gave far stronger concordance with individual binding sites than did GATA6-dependent gene expression in cell lines (Fig. 6A). HUG1N binding showed the highest concordance, but GATA6 occupancy in Caco2 and AGS cells also was significantly correlated. PCA ranking (Fig. 6D) also performed significantly better than GATA6-dependent gene expression (Fig. 6B) at assigning probable regulatory functions to GATA6 binding sites.

We compared the genes best correlated in PCA analysis (using z <−4 as a stringent measure; Fig. 6E shows the 20 highest scoring genes) with the most anti-correlated genes (z >+4, Suppl. Table 1). Among the 109 genes highly associated with GATA6 in expression space, 11 genes showed GATA6 occupancy within 20 kb of the TSS in all 3 cell lines (shaded yellow in Suppl. Table 1) and 31 genes bound GATA6 in 2 or more lines (highlighted in Figs. 6E and 6G, shaded green in Suppl. Table 1). Thus, within a reasonable distance for TF regulation of nearby genes (31, 32) the data reveal a high likelihood that GATA6 activates synexpressed genes. If GATA6 also represses transcription, we might expect to detect its binding near the 97 genes most anti-correlated with GATA6 expression (z>+4). However, none of these genes showed nearby GATA6 binding in all 3 lines (P=9e−4, Fisher exact test) and only 1 gene bound GATA6 in 2 cell lines (P=6.5e−9, Fig. 6F). The disparity in GATA6 binding near genes that are synexpressed and those that are anti-correlated in expression persisted for weaker associations (z <−3 or >+3, Fig. 6F, Suppl. Table 1). Together, these data argue for a predominantly activating function for GATA6. Moreover, as most GATA6-dependent, co-expressed genes do not bind GATA6 within 20 kb of the promoter, they may represent mainly indirect transcriptional targets.

Integrative analysis of GATA6 occupancy, gene dependencies, and tumor synexpression

At least 30% of GATA6-bound regions detected in only 1 cell line lacked canonical GATA motifs and may therefore represent indirect binding or spurious signals. By contrast, all of the hundreds of sites common to 2 or 3 cell lines carried a consensus GATA motif. Thus, shared sites reflect occupancy best, and although HUG1N and Caco2 cells shared only a fraction of binding sites (Fig. 5E), these might denote authentic, common target genes. Genes that control hormonal responses, cell proliferation, and lipid metabolism were significantly enriched near such shared binding sites (Fig. 7A). Indeed, 18 genes associated with cell proliferation showed similar GATA6 dependency and GATA6 binding in both cell lines (e.g., HES1 promoter, Fig. 7A). Eight of these likely bona fide GATA6 transcriptional targets encode TFs.

Figure 7

Integrative analysis of GATA6 occupancy, dependent genes, and synexpression groups to delineate primary transcriptional targets

A) Biological process Gene Ontology (GO) terms enriched among genes that bind GATA6 and require it for expression in both HUG1N and Caco2 cells. A Venn diagram illustrates the 18-gene overlap between 2 cell lines of GATA6-regulated and GATA6-occupied genes associated with cell proliferation and a representative data trace shows GATA6 binding at the HES1 promoter. B) Genes with a significant z-score in Principal Component Analysis (PCA) and nearby GATA6 occupancy in gastric cell lines. Yellow shading marks genes that also bind GATA6 in Caco2 CRC cells. C) Transcripts ascribed to the GO term “Transcriptional regulation” represent the most enriched category of genes with common GATA6 binding sites and levels that respond to GATA6 depletion in AGS, HUG1N and Caco2 cells. D) qRT-PCR for these 9 TF genes confirmed reduced expression in HUG1N and Caco2, and showed similarly reduced levels in GATA6-depleted AGS gastric cancer cells. Values represent the average fold-change ±SD in 3 independent experiments. PCA rank and z-score for each gene are listed, and genes showing extreme association with GATA6 in expression space are marked in bold. E) GATA6 occupancy at one of these 9 TF loci, CEPBG.

To identify targets most relevant to gastric cancer, we captured genes with high PCA rank and nearby GATA6 binding in HUG1N and AGS cells, irrespective of occupancy in Caco2 (Fig. 7B). Target genes determined by integrative analysis of expression and GATA6 binding encode diverse transcriptional, secretory, survival, cytoskeletal and metabolic factors, and as such, fail to deliver a unifying view of GATA6 function. Because binding in gastric cancer cells generally matched that in CRC cells (yellow shading in Fig. 7B, illustrated in Fig. 7D; see Venn diagram in Fig. 5E), we considered 733 sites present in all 3 cell lines, noting that 49 of the sites lie within 20 kb of genes that depend on GATA6 in both Caco2 and HUG1N cells. The Gene Ontology biological process associated with the highest statistical significance was Regulation of Transcription (Fig. 7C), as represented by TF genes CDX2, CEBPG, HES1, IRF8, LRRFIP1, MLXIPL, MYC, NR5A2, and RARB. All 9 of these transcripts were reduced in GATA6-depleted AGS, HUG1N, and Caco2 cells (Fig. 7D) and GATA6 bound the same sites in each cell line (e.g., 8 kb upstream of CEBPG, Fig. 7E). Five of these 9 TF genes showed very little expression in gastric cancer cell lines that lack both GATA6 and GATA4 (Suppl. Fig. S3E), and 4 of them fell within the highest percentile of PCA ranks, with z-scores <−3.36 (Fig. 7D). These data reveal a strong correlation with GATA6 in expression space. Other TFs (NR0B2, MYCN, KLF5) also were strongly correlated in expression but bound GATA6 in 1 or 2 cell lines, and the gene with the highest PCA rank, FOXE1 (z=−23), also encodes a TF. Identification of these candidate target genes might explain indirect GATA6 effects on gene regulation and suggests that it governs many genes through these TFs.

Functional redundancy among GATA factors

Because GATA6 and GATA4 recognizes the same DNA sequence, are often co-expressed, and have overlapping functions in some tissues (23, 24), we asked next if the two homologous TFs might share transcriptional targets in human gastric cancer. Lentiviral delivery of specific shRNAs depleted GATA4 efficiently in AGS and HUG1N gastric cancer cells (Fig. 8A), which we selected because they express both proteins (Fig. 2A). The same representative cell-cycle and TF genes that respond to GATA6 deficiency and bind GATA6 showed significantly reduced levels in GATA4-depleted cells (Fig. 8B). To test if the two factors are truly redundant, we examined expression of the same genes following shRNA-mediated depletion of GATA4, GATA6, or both in HUG1N cells. Some target genes showed similar effects upon loss of one or both TFs, but most transcripts showed at least an additive effect (Fig. 8C). The heightened sensitivity to dual TF depletion indicates that both GATA4 and GATA6 contribute toward regulation of a wide array of gastric cancer genes (Fig. 8D).

Figure 8

Additive and redundant effects of GATA6 and GATA4

A) Immunoblot evidence that GATA4 shRNA depletes GATA4 but not GATA6 protein in AGS (top) and HUG1N (bottom) gastric cancer cells. The lanes represent triplicate samples from 1 of 3 experiments. B) qRT-PCR analysis showing GATA4 dependency of most GATA6 target genes in HUG1N and AGS cells. GATA4 shRNA reduced GATA4 mRNA levels 17- (AGS) to 19- (HUG1N) fold and relative levels of candidate cell-cycle (green labels) and TF (brown labels) target transcripts of GATA regulation are shown from 3 independent replicates. C) qRT-PCR in HUG1N cells, showing that most target genes respond to GATA4 or GATA6 deficiency and more significantly to loss of both, revealing redundant activities. shRNAs reduced GATA6 and GATA4 mRNA levels about 30-fold and relative mRNA levels of target genes were measured in 3 replicates. D) Model for GATA function in gastric cancer. GATA6, in partial redundancy with GATA4, controls diverse genes, including those involved in cell differentiation and late phases of the cell cycle. Although some genes are direct transcriptional targets, others are controlled indirectly through GATA-factor regulation of intermediary TF genes, such as CDX2, HES1, and NR5A2.

DISCUSSION

GATA4 and GATA6 regulate diverse digestive epithelia (17, 23, 33, 34). Their genes are amplified in a high fraction of gastric and G-EJ adenocarcinomas and frequently overexpressed, even in the absence of gene amplification, and GATA6 is also amplified in pancreatic cancer, a related neoplasm (20, 21). Cell replication and colony assays suggest GATA6 dependence in esophageal and pancreas cancers (8, 9, 21) but the transcriptional basis of this dependence is unknown. Among 39 gastric cancer cell lines, GATA6 overexpression delineates a disease subset that expresses many genes in common, likely reflecting GATA6 regulatory functions. Genome-wide GATA6 occupancy, considered in relation to GATA6-dependent gene expression, provided one layer of insight into target genes. Transcripts that increased in GATA6-depleted cells, possibly owing to relief of GATA6-mediated repression, rarely showed GATA6 occupancy higher than the genome background. In contrast, genes with reduced expression in GATA6-deficient cells were enriched for nearby GATA6 binding and for functions related to cell cycle M-phase control, digestion, metabolism, and hormonal responses. Thus, in gastrointestinal cells, GATA6 coordinately activates genes of related function, with perhaps little direct transcriptional repression, distinct from its homolog GATA1, which activates and represses genes in blood cells (27). There is no obvious mechanism for this apparent difference, other than the transcriptional contexts in epithelial and blood cells.

Our consideration of gene synexpression across primary gastric and colorectal cancers significantly improved imputation of likely direct targets and functions. Most GATA6-dependent, synexpressed genes, however, lack GATA6 occupancy within 20 kb or more, implying indirect regulation through other TFs. Sequence-specific TFs that may fulfill this role are readily found among transcriptional targets that bind GATA6 in multiple cell lines and show extreme associations in expression space: CEBPG, IRF8, HES1, NR5A2 and CDX2; we propose that these TFs execute many cellular functions downstream of GATA6 (Fig. 8D). As TFs generally make poor drug targets, our observations have implications for the selective treatment of GATA6-overexpressing gastric cancers. Suitable targets for new drugs might include enzymes that modify components of this TF network or nuclear hormone receptors such as NR5A2.

GATA6 overexpression begins early in esophageal neoplasia (8) and the spectrum of transcriptional targets we identify supports early biologic effects. Control of M-phase genes and arrest of GATA6-depleted cells late in the cell cycle, for example, may explain accelerated cell replication. Target genes include CDX2, a TF that controls intestinal differentiation (35, 36) and partners with GATA6 to regulate intestinal cell growth (18). Presence of both CDX2 and GATA6 in many gastric cancers reflects the prevalence of intestinal differentiation in this disease (37) and their partnership may promote tissue-specific proliferation. Among other intestinal gene targets (e.g., FABP1, SPINK4, GCNT1, VILL), LGR5 in particular marks intestinal (38) and some gastric (39) stem cells. At least two TF target genes, MYC and NR5A2, are components of self-renewal and pluripotency networks (40) and we implicate GATA6 in regulating HES1, a transcriptional target of Notch signals (41). GATA6 activity at the HES1 promoter is noteworthy because Notch signaling plays an important role in mouse stomach epithelial stem and progenitor cell activity (42) and because HES1 and other GATA factors regulate each other’s expression in other tissues (43, 44). Thus, GATA6 targets in gastric cancer place its activity within pathways for cell replication and self-renewal, differentiation, and transcriptional control.

Epigenetic and signaling determinants of the GATA6 cistrome in each tissue and tumor are likely complex. Although mutations and other features may influence TF binding, we reason that a shared transcriptional program mediates GATA6-dependent oncogenic functions. We propose a core set of target genes that GATA6 occupies in multiple cell lines and that associate tightly with GATA6 expression in primary gastric cancers. The homologous protein GATA4 seems to regulate the same target genes redundantly with GATA6.

MATERIALS AND METHODS

Tissue samples and copy number analysis

SCNAs in gastric and G-E junction adenocarcinomas were reported recently (7) and supplanted with extensive public data from The Cancer Genome Atlas and GSE31168 series in the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/). Immunohistochemistry was performed on tumor specimens archived after fixation in formalin and embedding in paraffin. Gastric cancer samples (n=125) used for synexpression analysis were derived from 4 GEO series, GSE19826, GSE2109 and GSE13911, GSE22377, and 24 additional tumors assayed on Affymetrix HG U133 Plus 2.0 arrays. Colon cancer samples (n=1732) are from GEO series GSE10714, GSE13059, GSE13067, GSE13294, GSE13471, GSE14333, GSE17536, GSE17537, GSE17538, GSE18088, GSE18105, GSE20916, GSE2109, GSE21510, GSE23878, GSE26682, GSE26906, GSE28702, GSE31595, GSE33113, GSE33114, GSE4107, GSE4183 and GSE9348, and 11 additional samples. Transcript data in gastric and colon cancer cell lines are from the Cancer Cell Line Encyclopedia (http://www.broadinstitute.org/ccle/home), normalized by the GCRMA method (45).

Cell lines

All cell lines were obtained from the American Type Culture Collection. Caco-2, AGS, HTB135, NCIN87, and GC1Y were cultured in Dulbecco’s Modified Eagle Medium (DMEM) containing 10% fetal bovine serum (FBS). SNU5 and HGC27 cells were maintained in Iscove’s Modified Dulbecco Medium (IMDM) containing 20% FBS. HUG1N, SNU16 and ECC12 were cultured in Roswell Park Memorial Institute-1640 medium containing 15% FBS. All media were purchased from Life Technologies and supplemented with penicillin and streptomycin (Life Technologies).

RNA and protein expression

Total RNA was isolated using Trizol (Invitrogen), treated with the RNeasy Mini Kit (Qiagen), and DNA was removed using Turbo DNA Free (Ambion). For qRT-PCR analysis, 2 μg of total RNA was reverse transcribed with Superscript III First Strand Synthesis System (Invitrogen) and cDNA was amplified using SYBRGreen PCR Master Mix (Applied Biosystems). For transcript profiling, 1 μg of total RNA was processed for hybridization to U133A 2.0 oligonucleotide arrays for human genes (Affymetrix). To prepare nuclear lysates, cells were suspended in lysis buffer (20 mM Tris pH 7.9, 1.5 mM MgCl2, 10 mM KCl), followed by nuclear lysis buffer (20mM Tris pH 7.9, 1.5 mM MgCl2, 0.42M NaCl, 0.2 mM EDTA) in the presence of protease and phosphatase inhibitors. Whole-cell extracts were prepared in parallel to assess equal inputs for immunoblotting. Proteins resolved by SDS-PAGE were transferred to nitrocellulose membranes and probed with GATA6 (Cell Signaling, # 4253), GATA4 (Millipore, # AB4132), and Actin (Santa Cruz, # 47778) antibodies (Ab). Binding was detected by chemiluminescence after incubation with horseradish peroxidase-coupled secondary Ab (Santa Cruz, # 2048). Immunohistochemistry was performed as described previously (46) using CDX2 (Biogenex, # 10M MU392A-UC) and GATA6 Ab.

GATA6 and GATA4 knock-down, proliferation, and cell cycle analyses

Cells were infected with lentiviruses in the pLKO.1 vector (Open Biosystems) carrying GATA6 (5′-AGAACAGCGAGCTCAAGTATT-3′), GATA4 (5′-CCAGAGATTCTGCAACACGAA-3′), or a control shRNA not complementary to any human gene and known to lack cytotoxicity (NS, 5′-CCTAAGGTTAAGTCGCCCTCG-3′). Stable clones were selected in 2 μg/ml puromycin and GATA6 or GATA4 depletion assessed 10 days later by immunoblotting. Triplicate samples of 5,000 cells maintained in puromycin-free medium for 48h were seeded in 96-well plates and proliferation was assessed using CellTiter 96 (Promega G4000) at 570 nm absorbance using a Synergy-2 multi-detection plate reader (BioTek). For cell cycle analysis, cells were fixed in 70% ethanol, incubated overnight at −20°C, washed in PBS, stained with 10 μg/ml propidium iodide for 30 min, and analyzed on a FACScan instrument (Becton-Dickinson) at 488 nm.

Phospho-histone H3 (PH3) analyses

Histones were extracted as described (47), with substitution of 0.2 N HCl for 0.4 N H2SO4, neutralized with 3 volumes of Trizma base (pH 11), resolved by SDS-PAGE, and immunoblotted as described above with Histone H3 (Cell Signaling, cat. #9715) and PH3 (Cell Signaling, 3377) Ab. shRNA-treated AGS cells were also seeded on glass cover slips, fixed, incubated sequentially with PH3 Ab and FITC-conjugated goat anti-rabbit IgG (Jackson Immunoresearch, 111-096-144), and mounted in medium containing DAPI (Vectashield, H-1200).

Chromatin immunoprecipitation (ChIP) and ChIP-seq

ChIP and ChIP-seq were performed as described previously (18), using 3 μg GATA6 Ab (Cell Signaling 4253S). We used Cistrome tools (www.cistrome.org/) to call and annotate peaks, generate wiggle files and conservation plots, identify enriched sequence motifs and linked genes, and compare data across ChIP-seq libraries (48). Wiggle traces were generated using the Integrative Genome Viewer (IGV) (49).

Data analysis

Integrative expression analysis was performed in the R statistical computing environment (http://cran.r-project.org). Hierarchical clustering of the 1,000 most varying probesets in mRNA expression across gastric cancer cell lines was generated with a Euclidean distance metric and complete linkage clustering. A distinct cluster containing the HUG1N and AGS cell lines corresponded to those that express high levels of GATA6 (probeset 210002_at) and GATA4 (probeset 205517_at) mRNAs, independent of presence of GATA4 or GATA6 probesets in the distance matrix calculation. Functional categories enriched among groups of genes were determined using DAVID software (http://david.abcc.ncifcrf.gov/). Principal Component Analysis (PCA) was applied to the following 6 dimensions of measurements across all genes: average fold-changes in gene expression in Caco2 and HUG1N cells following shRA-mediated GATA6 depletion (measured by dChip, http://biosun1.harvard.edu/complab/dchip), and Pearson correlations of GATA6 mRNA levels with all mRNAs in primary stomach and colorectal cancers and cell lines. Fold-changes (FC) in transcript levels have a positive or negative sign to denote increases or decreases; for PCA analysis we centered fold-change values at 0 by the transformation FC -> sign(FC) (abs(FC) - 1.0). PCA axes are ordered by decreasing amount of variance explained in the 6-dimensional space (30) and the first 2 principal axes were used to derive scores for each gene, by projection. PCA scores for each axis were converted into z-scores and added to yield a single ranking that integrated shRNA and synexpression data.

Supplemental Figure S1. GATA4 and GATA6 amplifications and overexpression in human cancers

A) GISTIC analysis of public data on 16 tumor types shows that frequent, focal amplifications of GATA4 and GATA6 are most prominent in gastric adenocarcinoma. Chromosomes 8 and 18 are shown in detail. B) Additional representative immunohistochemical analysis of GATA6 and CDX2. Expression of both proteins is shown in a gastric tumor that lacked GATA6 amplification and in a CRC. Stained tissue sections were photographed at 10X magnification and insets are digitally magnified 6X further. C) Aggregate analysis of immunohistochemical staining for GATA6 and CDX2 in 22 gastric cancers, including 4 tumors that carried ≥4 copies of GATA6. Amplification was estimated conservatively from qPCR analysis of genomic DNA using 2 separate primer pairs each for GATA6 and control TUBB loci. Mod, moderate. D) DNA heat map of cases from The Cancer Genome Atlas with the highest amplification of GATA6. As these copy-number changes tend to result from highly focal amplifications, PCR estimations in Fig. S1C are unlikely to represent an artifact of aneuploidy.

Supplemental Figure S2. Infrequent GATA6 amplification in gastric cancer cell lines

A) GISTIC analysis of SCNA data from 42 gastric cancer cell lines (top) and 522 tumor samples (bottom) shows the low frequency of amplification in cell lines compared to primary cancers. B) Table of GATA6 and GATA4 gene copy numbers in SNP array analysis of 27 gastric cancer cell lines. CN, copy number.

Supplemental Figure S3. Biological functions enriched among candidate GATA6 target genes. A–B)

Functional categorization of genes that depend on GATA6 for their expression in HUG1N gastric (A) and Caco2 CRC (B) cells. C–D) Top 8 biological processes (Gene Ontology terms) enriched among genes that depend on GATA6 for expression in HUG1N (C) or Caco2 (D) cells and bind GATA6 within 20 kb of the TSS in that cell line. GATA6 binds near genes implicated in cell proliferation and in digestive processes. E) qRT-PCR analysis of candidate TF target genes in HTB135, HGC27 and GC1Y, gastric cell lines that lack GATA6 and express little to no GATA4 (see Fig. 2A). Low basal expression of these genes in GATA-deficient cell lines supports their position as direct transcriptional targets.

Suppl. Table 1 ( Integrative analysis of genes with nearby GATA6 occupancy in Caco2, HUG1N and AGS cells; expression changes (Δ exp) in Caco2 and HUG1N cells upon shRNA-mediated GATA6 depletion; and synexpression (Synexp) with GATA6 in gastric and colorectal primary tumors and cell lines. DNA occupancy is denoted as present (1) or absent (0); Δ exp values represent log2 fold-changes in cells treated with GATA6 shRNA compared to cells treated with a control shRNA; and Synexp values represent Pearson correlations of GATA6 mRNA levels with the corresponding gene across all examined samples. PCA scores were calculated only for Δ exp, which accounted for 53% of the variation, and for the combination of Δ exp and Synexp. Cumulative scores in the far right column yield this ranked list of candidate GATA6 transcriptional targets. The bold blue boxes contain highly correlated (PCA score <−4, top of table) and highly anti-correlated (PCA score >4, bottom of table) genes. Those contained within the thin black box show moderate GATA6 correlations (PCA score <−3). Yellow shading highlights nearby GATA6 binding in all 3 cell lines and green shading highlights binding in any 2 of the 3 cell lines. Because some genes are represented by multiple Affymetrix probe sets (here designated as Gene X_probe 2, Gene X_probe 3, etc.), gene numbers indicated in the text do not correspond exactly to the numbers of rows in this table.

Suppl. Table 2. Primers used for quantitative PCR and RT-PCR analysis