This article is mentioned in:

Introduction

Breast cancer (BC), is a type of malignancy that poses the greatest threat to female health worldwide, despite its relatively long 5-year survival rate (~58.75% between 1972 and 2011 in China) (1). BC is one of the most common types of cancer in Asia, where ~39% of all worldwide cases are diagnosed (2). The treatment and care for patients with BC is a large economic burden in developed and developing countries (3–5). Considering the health and economic effects of BC, numerous studies have been conducted to investigate the pathogenesis, and to identify therapeutic targets and potential factors that may contribute its progression. However, the mechanisms behind the development and progression of BC remain poorly understood, and additional studies are required.

G9A [also known as euchromatic histone lysine methyltransferase 2 (EHMT2)], encodes a methyltransferase that methylates lysine residues of histone H3. The dysregulation of G9A is associated with numerous abnormal biological processes and the emergence of diseases, including cancer. For example, the downregulation of G9A could trigger DNA damage responses and inhibit the progression of colorectal cancer (6). In head and neck squamous cell carcinoma, G9A dysregulation has been associated with epithelial-mesenchymal transition (EMT)-mediated metastasis and the maintenance of cancer stem cell-like characteristics (7). Furthermore, it has been identified that G9A performs important roles in the transition of acute-to-chronic pain (8). In BC, the aberrant expression of G9A could affect the stability of SRY-box 2 (SOX2), a tumor suppressor gene, by inducing its epigenetic silencing (9). Furthermore, a study by Si et al (10) revealed that G9A performed important roles in the dysfunction of the reciprocal feedback loop between GATA binding protein 3 and zinc finger E-box binding homeobox 2, and that it contributed to the progression of BC. The combined analysis of its genome-wide profiles and the gene expression data following G9A-knockdown (KO) or overexpression in BC could improve our understanding of its regulatory functions, and has not previously been investigated.

MicroRNAs (miRNAs/miRs) are small non-coding RNA molecules, 22–25 nucleotides in length. miRNAs are formed from the precursor miRNA, a double-stranded hairpin RNA or even the intergenic region, and perform important roles in the regulation of gene expression at the post translation level by silencing mRNA translation or inducing protein degradation (11). Numerous miRNAs have been shown to contribute to the progression of a number of cancer types. For example, the downregulation of miR-711 could inhibit the cellular proliferation of BC, and miR-711 has been shown to be an independent prognostic factor for BC (12). miR-27a/b and miR-494 may regulate the expression of tissue factor pathway inhibitor α (TFPIα) in its estrogen-mediated downregulation in BC (13).

In the present study, through the combined analysis of genome-wide profiles of G9A obtained through chromatin immunoprecipitation-DNA sequencing (ChIP-Seq) and the gene expression dataset [from the Gene Expression Omnibus (GEO)] following G9A-KO in BC, potential biomarkers were identified for the treatment of BC. miRNAs that may regulate those biomarkers were screened out and the core miRNAs and genes were obtained. The findings of the present study may be valuable for the diagnosis and treatment of BC.

Materials and methods

ChIP-Seq and gene expression dataset

In the present study, the ChIP-Seq and gene expression datasets were comprehensively obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo). The ChIP-Seq dataset (GSM1642516) was deposited from the study by Si et al (10), which used the G9A-specific antibody to capture the DNA fragments in MCF-7 cells and sequenced those fragments based on the GPL11154 Illumina HiSeq 2000 (Homo sapiens). Furthermore, an input control sample, which extracted all DNA in the human breast adenocarcinoma MCF-7 cell line, was adopted for the background correction in the identification of G9A binding sites (peaks). The gene expression dataset GSE34925 (14) contained 3 G9A-KO and 3 control BC samples, which were based on GPL6244 (HuGene-1_0-st) Affymetrix Human Gene 1.0 ST Array [transcript (gene) version] (Affymetrix, Inc., Santa Clara, CA, USA).

Differential expression analysis

The raw CEL data (data storage format) were imported to R based on the affy package (15). Background correction and expression normalization was conducted using the oligo package (16). Probe IDs were converted to gene symbols via the annotation package of the microarray platform and the expression values were summarized that corresponding to multiple probe sets. A paired Student's t-test was conducted between logarithmic transformed expression values of G9A-KO samples and control samples. Differentially-expressed genes (DEGs) with a fold-change expression value of >1.41 or <0.71, and a false discovery rate corrected P-value of <0.05 were screened out.

Functional enrichment analysis

The Database for Annotation, Visualization and Integrated Discovery (https://david.ncifcrf.gov/) (17) was used for the functional analysis of the DEGs. Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.kegg.jp/) pathways that satisfied the threshold of P<0.05 were screened out.

Genome-wide binding analysis of G9A

Bowtie 2 (18), an ultrafast genome aligning tool, was used for the mapping of the raw reads to UCSC hg19 genome with a maximum of 2 mismatches in every read. Model-based analysis of ChIP-Seq version 14 (macs14) (19) was used for the identification of the binding sites of G9A (also known as call peak) with a P-value of <1×10−5. Peaks located in −3000 to 3000 bp centered in the transcription start sites (TSS) were considered to be conpeaks. Based on the ChIPseeker package (20), the nearby genes, genome features (including promoter and 5′-untranslated region) and distances to the nearest TSS were assigned to the peaks.

Construction of miRNA-gene network

The overlapping genes between DEGs and genes located in the conpeaks were screened out. miRNAs that may regulate those overlaps were identified using the TargetScan database (http://www.targetscan.org/) (21). The miRNA-gene regulation network was visualized through Cytoscape (http://www.cytoscape.org/).

Results

DEGs

A total of 217 DEGs were identified in the G9A-KO BC samples, and the distribution of these DEGs is demonstrated in Fig. 1. The supervised clustering (Fig. 2) based on those DEGs was able to separate the G9A-KO and control samples well.

Figure 1. Distribution of DEGs in euchromatic histone lysine methyltransferase 2-knockout samples. The red plus signs and blue triangles represent upregulated and downregulated genes, respectively. The black circles represent non DEGs. FC, fold-change.

Figure 2. Two-way supervised clustering of differentially-expressed genes and samples using R software. G9A, euchromatic histone lysine methyltransferase 2; KO, knockout.

Enriched functions

A total of 172 GO terms and 3 KEGG pathways were identified to be enriched in DEGs, and were mainly involved in the cell process, the inflammatory response and cancer progression, for example, the extracellular region part. The top 10 GO terms according to P-value and KEGG pathways are presented in Tables I and II, respectively.

Table I. Top 10 GO terms according to the P-value.

Table I.

Top 10 GO terms according to the P-value.

Category	GO term	Count	P-value
CC	Extracellular space	27	5.22×10−8
CC	Extracellular region part	32	9.50×10−8
BP	Organ development	44	1.04×10−6
BP	Anatomical structure development	53	1.44×10−5
BP	System development	50	1.59×10−5
BP	Integrin-mediated signaling pathway	8	1.70×10−5
CC	Extracellular region	44	1.95×10−5
BP	Inflammatory response	15	2.81×10−5
BP	Blood vessel development	13	3.09×10−5
BP	Defense response	21	3.32×10−5

[i] GO, gene ontology; CC, cellular component; BP, biological process.

Table II. Enriched Kyoto Encyclopedia of Genes and Genomes pathways.

Table II.

Enriched Kyoto Encyclopedia of Genes and Genomes pathways.

Pathway name	Count	P-value	Genes
Cytokine-cytokine receptor interaction	13	2.99×10−4	CSF2, TNFSF4, IL24, IL7R, CXCL11, TNFSF18, CXCL10, KDR, INHBB, VEGFC, CXCL16, IL1A, GHR
Focal adhesion	11	5.24×10−4	VEGFC, ITGB8, LAMA5, ITGB4, ITGA2, RELN, BIRC3, VAV1, MYLK, KDR, SPP1
ECM-receptor interaction	6	6.80×10−3	ITGB8, LAMA5, ITGB4, ITGA2, RELN, SPP1

[i] ECM, extracellular matrix.

Genome-wide binding profile of G9A

A total of 10,422 binding sites (peaks) of G9A were identified, including 1,210 conpeaks located in −3000 to 3000 bp centered in the TSS. The distribution of all peaks across the genome is shown in Fig. 3. In addition, 1,138 genes were identified to be located in the 1,210 conpeaks and 15 of these genes overlapped with DEGs (Table III).

Figure 3. Distribution of peaks of euchromatic histone lysine methyltransferase 2 across the genome. UTR, untranslated region.

Table III. Overlapping genes between differentially-expressed genes and genes located in conpeaks.

Table III.

Overlapping genes between differentially-expressed genes and genes located in conpeaks.

Gene	logFC
SLCO4C1	−0.65617
ARHGAP24	−0.66706
ANK2	−0.55143
TMX4	0.54535
LARP4	−0.50197
TC2N	−0.78340
TMTC1	0.52521
SLFN5	−0.62773
CYTIP	0.67537
PLCB4	−0.55811
ATF7IP2	−1.19761
CAMK2D	−0.50487
TLR4	0.85905
INSIG1	−1.18712
GHR	−0.60367

[i] FC, fold-change.

miRNA-gene regulation network

Based on TargetScan, a total of 35 miRNAs were identified to regulate the 15 overlapping genes, and 65 miRNA-gene pairs were obtained among them. The miRNA-gene regulation network is shown in Fig. 4. The top 10 genes according to their connectivity in the miRNA-gene network are listed in Table IV. In the network, insulin-induced gene 1 (INSIG1) was regulated by 9 miRNAs and had the highest connectivity, which indicated that it was a hub gene. Furthermore, two other public datasets (GSE29044 and GSE36774) were downloaded from the GEO database, which including data from breast cancer and normal tissue samples. These two datasets were background corrected and quantile normalized using the preprocessCore package of R (http://www.bioconductor.org/packages/release/bioc/html/preprocessCore.html). The Pearson and Spearman's correlation indices were respectively calculated using R between relative INSIG1 and G9A mRNA levels. The results indicated a significantly positive correlation (Fig. 5) between G9A and INSIG1, which may demonstrate their important roles in the progression of BC. The top 10 genes according to their connectivity in miRNA-gene network are presented in Fig. 4.

Figure 4. miRNA-gene regulation network. Rectangles show genes, whilst miRNAs are shown in ovals. Arrows represent the interactions between these genes and miRNAs.

Figure 5. Expression values of G9A were plotted against that of INSIG1 based on the other two public datasets (GSE29044 and GSE36774). G9A, euchromatic histone lysine methyltransferase 2; INSIG1, insulin-induced gene 1.

Table IV. Top 10 genes according to their connectivity in the microRNA-gene network.

Table IV.

Top 10 genes according to their connectivity in the microRNA-gene network.

Gene	Connectivity
INSIG1	9
ARHGAP24	8
CAMK2D	8
GHR	8
ANK2	7
TLR4	6
TMX4	6
CYTIP	4
miR-141	4
miR-145	4

Discussion

BC is a complicated disease, and numerous factors may contribute to its initiation and progression, including occupational exposure (22), variations in expression level or mutations of specific genes (23,24). The traditional treatment options for BC mainly involve surgical resection, chemotherapy and radiotherapy, which have been shown to be effective but rarely eradicate the disease. With the development of next-generation sequencing technology, numerous biomarkers of BC have been identified and certain novel therapeutics have been proposed (25). Immunotherapy is one of the most popular therapeutic treatment options (26,27). The identification of novel biomarkers may therefore improve the prognosis for patients with BC.

In the present study, the genome-wide profiles of G9A in genome wide were identified using ChIP-Seq. DEGs in G9A-KO BC samples compared with the control samples were screened out. The combined analysis obtained potential direct targets of G9A, the majority of which were downregulated in the G9A-KO BC samples, which may indicate its upregulated functions in BC.

DEGs in G9A-KO samples are primarily involved in biological processes or pathways associated with cellular processes, inflammatory responses and cancer progression, including extracellular region part, cytokine-cytokine receptor interaction and extracellular membrane (ECM) receptor interaction, all of which are closely associated with the development of BC. In particular, inflammatory triple-negative BC, a type of rare and aggressive BC with a significantly poorer 5-year survival rate compared with other types of BC, is believed to mainly be attributed to an abnormal inflammatory response (28,29). Previous studies have considered that the ECM-receptor interaction pathway may be associated with the progression of BC, which could indicate the repeatability of the present study (30,31).

Among the 1,138 genes distributed in the 1,210 peaks located between −3000 to 3000 bp centered in the TSS, 15 were identified to overlapped with the DEGs in the G9A-KO BC samples. Furthermore, 11 of those overlaps were revealed to be downregulated in G9A-KO samples compared with the control samples, which may indicate that G9A is mainly involved in the upregulation of its target genes. The positive correlation between the expression values of G9A and INSIG1 was verified using two public datasets (Fig. 5). Furthermore, INSIG1 possessed an intimate connectivity pattern in the miRNA-gene regulation network (Fig. 4). Thus, this suggests that G9A and INSIG1 serve important roles in the progression of BC. INSIG1, which is downregulated in G9A-KO BC, was closely associated with the progression of several types of cancer via the regulation of glucose metabolism and hypoxia-induced EMT (32,33). BC had a strong association with insulin resistance, which was mediated by INSIG1 (34). In previous studies, INSIG1 was also demonstrated to be differentially expressed in BC (35,36). Combined with the results of the present study, this indicates the important roles of INSIG1 in BC.

miR-203 regulated 3 genes in the miRNA-gene regulation network, including INSIG1, which has been reported as a tumor suppressor in numerous types of cancer, including BC (37–39). The activity of miR-203 was also affected by a number of other factors, including kallistatin (40). Therefore, it would be useful to study the upstream and downstream regulation loops of miR-203 for the diagnosis and treatment of BC, as well as other types of cancer. Other nodes (miRNAs or genes) with high connectivity in the miRNA-gene regulation network may also be potential biomarkers for BC, which should be additionally verified through molecular biological experiments.

In conclusion, the present study conducted the combined analysis of ChIP-Seq and gene expression profiles in BC, and identified potential biomarkers via screening. The present results may aid the development of novel diagnostic or treatment methods, including immunotherapy, and therefore improve the prognosis of BC.

Acknowledgements

The present study was supported by the Municipal Science and Technology Commission of Tianjin (grant nos. 15ZLZLZF00440 and 16ZLZXZF00120) and the Health Bureau Science and Technology Foundation of Tianjin (grant no. 2014KZ102).

References