A flow diagram of the study is presented in Fig. 1. The raw data on patients with gastric cancer containing RNA sequencing (RNA_seq) and clinical information were obtained from The Cancer Genome Atlas (TCGA) repository website (https://cancergenome.nih.gov/). The R ‘limma’ Bioconductor package (12) was adopted to screen the DEGs between normal gastric and tumour tissues based on the following criteria: Fold change (FC), |log₂(FC)|>1; and false discovery rate (FDR) <0.05. To minimize noise in the DEG dataset, in strict accordance with the analysis method on the WGCNA official website (https://labs.genetics.ucla.edu/horvath/htdocs/CoexpressionNetwork/Rpackages/WGCNA/), genes with too many missing values and unknown names were checked for and excluded, and genes for which the names corresponded to multiple probes were deleted. The samples were clustered to exclude the obvious outliers (7,13). The microarray dataset of filtered probe sets was constructed using the remaining genes and samples. In addition, clinical variables, including age and gender, pathological variables [including drug response, tumour (T) stage and tumour grade] and survival information were compiled for the WGCNA analysis. The DEGs were normalized by log2 transformation prior to the subsequent analyses. A total of two investigators separately collated and analysed the data, and a joint decision was reached in cases of disagreement.

Construction of a weighted co-expression network represents an effective method for identifying modules and for defining the intra-module connectivity. In the current study, WGCNA was performed on the microarray data of filtered probe sets using the R ‘wgcna’ Bioconductor package (13). Each paired gene-gene association was correlated utilizing the absolute value of the Pearson product moment correlation (a gene co-expression similarity measure); the absolute value represents the co-expression similarity. An adjacency matrix was constructed utilizing a ‘soft’ power adjacency function, where the ‘soft’ power indicates that the resulting adjacency measures the connection strength. Clusters of co-expressed genes were designated by hierarchical cluster analysis following subtraction of the topological overlap measure of similarity of 1 (7,10).

The modules of WGCNA are groups of highly correlated genes. In network terminology, modules are sets of genes with high topological overlap. To obtain the potential associations between the co-expressed gene clusters and clinical variables, a single column of vectors termed the module eigengenes (MEs) (7) was used. The MEs were generated by reserving the first principal component of a given module and were representative of the gene expression profiles in a module. Since each ME contains the majority of the variance in the original data, it represents a summary measure for the overall co-expression network. The consistency between the expression of a particular gene and the ME expression is termed the module membership. This measure of co-expression network centrality is decided by calculating the Pearson correlation coefficient between each individual gene and the ME. Further details on WGCNA theory and algorithm have been previously published (7,10).

To define the pivotal module, the correlations between the MEs and clinical features were calculated. Furthermore, functional enrichment analysis and protein-protein interaction (PPI) analysis was conducted on the genes in the yellow module to examine the probable mechanism underlying the impact of the genes on chemotherapy resistance and the prognosis of gastric cancer. The Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david-d.ncifcrf.gov/) was used to analyse the enriched biological process (BP) terms and pathways (14). The top 15 most significant BP terms (P≤0.05) are presented. Only those KEGG pathways with P≤0.05 and ≥10 enriched genes were considered to be significant. In addition, the PPI networks were constructed using STRING (https://string-db.org/).

The human gastric cancer cell lines AGS and MKN45 were purchased from the American Type Culture Collection (Manassas, VA, USA). All the cell lines were maintained in RPMI 1640 medium (Gibco; Thermo Fisher Scientific, Inc., Waltham, MA, USA) supplemented with 10% foetal bovine serum (ScienCell Research Laboratories, Inc., San Diego, CA, USA) and 2 mM penicillin and streptomycin, and were cultured in a humidified 37°C incubator at 5% CO₂.

The GLIS2 overexpression plasmids were constructed and cloned into the CV061 vector (Shanghai GeneChem Co., Ltd., Shanghai, China) between the EcoRI and HindIII sites. The primer sequence used for GLIS2 was: 5′-ACGGGCCCTCTAGACTCGAGATGCACTCCCTGGACGAGCCGCTCG-3′. AGS and MKN45 cells were transfected with GLIS2 overexpression plasmids using Lipofectamine^® 3000 (Invitrogen; Thermo Fisher Scientific, Inc.). According to the manufacturer's instructions, 1×10⁵ gastric cancer cells were inoculated onto a 24-well plate and 1 µg/ml of the plasmid was added, and the transfected cells were incubated in a 5% CO₂ incubator at 37°C for 48 h. Subsequent treatments were performed when the cells had grown to 80–90% confluence.

For the western blotting assay, cells were lysed using radioimmunoprecipitation assay and phenylmethylsulfonyl fluoride buffers (100:1; Sigma-Aldrich; Merck KGaA, Darmstadt, Germany), and the lysates were collected by centrifugation at 12,000 × g for 12 min at 4°C. The protein concentration was determined using a bicinchoninic acid protein assay kit (Sigma-Aldrich; Merck KGaA). Equal amounts of protein lysates (30 µg) were electrophoresed on 10% SDS-polyacrylamide gels and transferred to methanol-activated polyvinylidene fluoride (PVDF) membranes. The PVDF membranes were blocked with 5% non-fat dried milk in Tris-buffered saline (pH 7.4) containing 0.1% Tween-20 (TBST) for 1 h, and incubated with the primary GLIS2 (Abcam, Cambridge, UK; cat. no. ab28462; diluted 1:750) and GAPDH antibodies (Abcam; cat. no. ab8245; diluted 1:750) in TBST at 4°C overnight. The membranes were subsequently washed three times with TBST and incubated with a secondary Alexa Fluor-conjugated anti-rabbit antibody (Abcam; cat. no. ab150077; goat anti-rabbit IgG H&L; diluted 1:5,000) for 1 h at room temperature. Signal detection was performed for 5 min using an enhanced chemiluminescence reaction (Dalian Meilun Biotechnology Co., Ltd., Dalian, China). The GAPDH antibody was used to normalize protein expression.

Total RNA was isolated from cells using TRIzol^® reagent (Invitrogen; Thermo Fisher Scientific, Inc.). RT was performed (on 500 ng of RNA) using the PrimeScript RT Reagent kit (Takara Biotechnology Co., Ltd., Dalian, China), and the RT-qPCRs were performed using SYBR Premix Ex Taq (Takara Biotechnology Co., Ltd.), according to the manufacturer's protocol. The following primers were used: GLIS2 forward, 5′-CTTCGGGGAGGCTGGATT-3′ and reverse, 5′-GGTGATACTCAGCTTCAGGTCG-3′; and GAPDH forward, 5′-AATCCCATCACCATCTTCCAG-3′ and reverse, 5′-GAGCCCCAGCCTTCTCCAT-3′. The relative expression of the genes (GLIS2 and GAPDH) was calculated using the 2^−ΔΔCq method (15).

Fluorouracil (5-FU), cisplatin and doxorubicin were purchased from Selleck Chemicals (Houston, TX, USA). Cells (6×10³ per well) were seeded in 96-well plates and incubated for 48 h. 5-FU, cisplatin and doxorubicin were added at final concentrations of 0, 2, 4, 8, 16 and 32 µmol/l with dimethyl sulfoxide (DMSO) alone serving as a control. The cells were incubated in a 5% CO₂ incubator at 37°C for 24 h. The medium was removed, and the cells were cultured in fresh growth medium containing 0.5 mg/ml MTT for 4 h under the above conditions. The supernatant was removed, and 100 µl DMSO was added to dissolve the formazan crystals. The absorption was measured at 490 nm, and the results were expressed as the half-maximal inhibitory concentration (IC50), the concentration of chemotherapeutic drug required for 50% inhibition in vitro.

The majority of the visualizations were generated using R version 3.3.1 (https://www.r-project.org/) except for the KEGG network and PPI visualizations, for which the ClueGO (16) and STRING online tools were used. Survival analysis and clinical correlation analysis was performed using TCGA expression and clinical data on gastric cancer, setting the median expression value of GLIS2 as a cut-off value between high and low expression. The Kaplan-Meier method was used to estimate survival, and a log-rank test was used to assess differences between the survival curves. Furthermore, the prognostic role of GLIS2, including its role in predicting overall survival and disease-free survival, was validated using online tools (http://kmplot.com/analysis/) (17). Statistical analyses were implemented using GraphPad Prism 6.0 (GraphPad Software, Inc., La Jolla, CA, USA). All data are presented as the mean ± standard deviation. Student's t-tests were used for comparisons between groups, except for the analysis of clinicopathological features, for which the two-sample t-test was used. P<0.05 was considered to indicate a statistically significant difference.

As presented in the flow chart (Fig. 1), public RNA_seq datasets were downloaded from TCGA repository website, containing 238 gastric cancer samples and 33 normal tissue samples. Based on the criteria of |log₂(FC)|>1 and FDR<0.05, a total of 5,630 DEGs were screened out. Prior to constructing the microarray dataset of filtered probe sets, 204 low-quality DEGs and 26 cancer samples were excluded (data not shown). To facilitate the subsequent WGCNA analysis, the clinical information was compiled and converted. To obtain an overview of the clinical information, hierarchical cluster analysis was performed (data not shown).

The remaining 5,426 most varying and most connected genes from 212 gastric cancer samples were obtained for construction of the co-expression network. As presented in Fig. 2, 16 distinct gene modules with high topological overlap were identified when the soft threshold power β was set to 6. The number of genes included in the modules ranged between 36 (light cyan) and 3,068 (turquoise). Each module was assigned a different colour to distinguish between them, and the grey module was reserved for genes that were identified as not being co-expressed.

The yellow module had the highest correlation with drug response (r=0.31; P<0.001) and days to mortality (r=−0.3; P<0.001) among the module-trait relationships (Fig. 3A), and was selected for study in the subsequent analyses due to its clinical importance. The yellow module was also associated with new tumour type (including locoregional recurrence, distant metastasis and new primary tumour; r=0.39; P<0.001). It was additionally identified that the red module was associated with histological type (r=−0.24; P<0.001), the blue module was associated with tumour histological grade (r=0.23; P<0.001), and the tan module was associated with microsatellite instability (r=0.37; P<0.001). In addition, the intra-modular connectivity was calculated for each gene based on its Pearson correlation with all other genes in the module (Fig. 3B), to help prove the importance of these modules.

To investigate the biological importance of the yellow module, 174 genes in the yellow module were functionally enriched using DAVID for BP analysis and KEGG pathways analyses. The BPs of the yellow module were primarily enriched in a number of aspects, including ‘cell adhesion’ (P<0.001), ‘vasculature development’ (P<0.001), and ‘regulation of cell proliferation’ (P<0.001. Fig. 4A). In addition, the DAVID-based enrichment analysis identified KEGG pathways for the 174 genes in the yellow module (Fig. 4B and C). Among them, certain tumorigenesis-associated signalling pathways were present, including the TGF-β (P=0.006) and p53 (P=0.036) signalling pathways. In particular, cell adhesion-associated signalling pathways, including extracellular matrix (ECM)-receptor interactions (P<0.001) and focal adhesion (P<0.0001), were observed to serve a key role in the regulation of drug resistance in gastric cancer.

The present study proceeded to develop the PPI network and investigate the hub genes of the yellow module (Table I). To identify novel proteins and pathways modulated by the pivotal genes, GLIS2 was selected as a candidate gene for further analysis and verification subsequent to comparing the data and reviewing the literature. Analysis of the TCGA dataset demonstrated that high expression of GLIS2 was significantly associated with tumour chemotherapeutic resistance (P=0.031; Fig. 5A). Furthermore, high GLIS2 positivity in gastric cancer samples corresponded to a significantly worse prognosis (P=0.003; Fig. 5B). Similar results were obtained using the online tools to verify the prognostic role of GLIS2. Patients with higher GLIS2 positivity had significantly shorter overall survival (P<0.001; Fig. 5C) and disease-free survival (P<0.001; Fig. 5D) times. In addition, the association between GLIS2 and other clinicopathological characteristics was evaluated in patients with gastric cancer (Table II). The expression of GLIS2 was significantly associated with histological type (P=0.039), histological grade (P=0.02), pathological stage (P=0.013), T stage (P=0.002) and microsatellite instability (P<0.001).

To examine whether the overexpression of GLIS2 enhanced chemotherapy tolerance in gastric cancer, MTT assays were performed. Fig. 6A-D demonstrates that the relative mRNA and protein expression levels of GLIS2 were markedly increased following transfection of GLIS2 plasmids in the AGS and MKN45 cell lines, respectively (P<0.05). Following treatment with different concentrations of three chemotherapeutic agents for 24 h, the results indicated that gastric cancer cells that overexpressed GLIS2 had increased survival in response to cisplatin, fluorouracil and doxorubicin compared with parental cells with lower GLIS2 expression. Specifically, the growth activity of the AGS experimental groups were higher compared with the control group (Fig. 6E; P<0.05), and consistent results were also obtained in the MKN45 cell line (Fig. 6F; P<0.05).