Systems-level effects of ectopic galectin-7 reconstitution in cervical cancer and its microenvironment

Abstract

Background

Galectin-7 (Gal-7) is negatively regulated in cervical cancer, and appears to be a link between the apoptotic response triggered by cancer and the anti-tumoral activity of the immune system. Our understanding of how cervical cancer cells and their molecular networks adapt in response to the expression of Gal-7 remains limited.

Methods

Meta-analysis of Gal-7 expression was conducted in three cervical cancer cohort studies and TCGA. In silico prediction and bisulfite sequencing were performed to inquire epigenetic alterations. To study the effect of Gal-7 on cervical cancer, we ectopically re-expressed it in the HeLa and SiHa cervical cancer cell lines, and analyzed their transcriptome and SILAC-based proteome. We also examined the tumor and microenvironment host cell transcriptomes after xenotransplantation into immunocompromised mice. Differences between samples were assessed with the Kruskall-Wallis, Dunn’s Multiple Comparison and T tests. Kaplan–Meier and log-rank tests were used to determine overall survival.

Results

Gal-7 was constantly downregulated in our meta-analysis (p < 0.0001). Tumors with combined high Gal-7 and low galectin-1 expression (p = 0.0001) presented significantly better prognoses (p = 0.005). In silico and bisulfite sequencing assays showed de novo methylation in the Gal-7 promoter and first intron. Cells re-expressing Gal-7 showed a high apoptosis ratio (p < 0.05) and their xenografts displayed strong growth retardation (p < 0.001). Multiple gene modules and transcriptional regulators were modulated in response to Gal-7 reconstitution, both in cervical cancer cells and their microenvironments (FDR < 0.05 %). Most of these genes and modules were associated with tissue morphogenesis, metabolism, transport, chemokine activity, and immune response. These functional modules could exert the same effects in vitro and in vivo, even despite different compositions between HeLa and SiHa samples.

Conclusions

Gal-7 re-expression affects the regulation of molecular networks in cervical cancer that are involved in diverse cancer hallmarks, such as metabolism, growth control, invasion and evasion of apoptosis. The effect of Gal-7 extends to the microenvironment, where networks involved in its configuration and in immune surveillance are particularly affected.

Background

Galectins are a family of carbohydrate-binding proteins involved in immune response, angiogenesis and cancer development [1]. Different members of the galectin family, such as galectin-1 (Gal-1) and galectin-7 (Gal-7), have been found to influence the pathogenesis of cervical cancer (CaCx). The expression of Gal-1 is functionally linked to histopathological grading in cervical cancer patients, namely by affecting the rate of proliferation, lymph node metastasis and tumor invasion [2]. Gal-1 is also the receptor for Trichomonas vaginalis [3], a sexually transmitted protozoan parasite and risk factor for cervical cancer [4]. In contrast, Gal-7 is preferentially expressed in stratified squamous epithelia from skin, genital and upper digestive track [5]. Gal-7 is a p53-inducible gene, which is upregulated in response to UVB radiation in normal human keratinocytes [6].

As an endogenous lectin, a fraction of Gal-7 is constitutively localized at the mitochondria. It has been found to interact with the anti-apoptotic protein Bcl-2, suggesting its regulatory role in apoptotic processes [7]. Importantly, increased Gal-7 expression has been shown as a positive predictive biomarker for clinical responses after adjuvant radiation therapy in cervical cancer patients [8]. While a plethora of distinct properties of Gal-7 are known, an integrative analysis of the molecular mechanisms with which Gal-7 expression shapes the tumorigenic process has not yet been performed.

In the present study, we performed an integrative analysis of the impact of Gal-7 reconstitution in cervical cancer cells and their microenvironment at the systems level in silico, in vitro, and in a mouse model. For that purpose, we conducted a meta-analysis of a whole spectrum of clinical data in which cervical cancer patients showed a significant longer life span when tumors had simultaneous high Gal-7 and low Gal-1 expression. To validate these observations in the biological system, we ectopically expressed Gal-7 in CaCx cell lines and evaluated them through transcriptomics, SILAC-based proteomics, gene methylation profiling, and network analysis. We identified numerous circuits implicated in cancer hallmarks that were affected by Gal-7 re-expression (Fig. 1a). These results suggest a bi-directional regulation between the tumor and its microenvironment where Gal-7 could be a critical mediator.

Fig. 1

Study design and galectin expression in cervical cancer. a Pipeline of the complete experimental approach. b Gene expression profiles of galectin family members in clinical samples (normal cervix and squamous cell carcinoma, SCC) and CaCx cell lines obtained from the Scotto cohort. c Gal-7 transcription in normal cervical tissue, high grade squamous intraepithelial lesions (HSIL) and squamous cell carcinomas (SCC) from the Zhai cohort. d qPCR and Western blot analysis of Gal-7 in primary keratinocytes (PK), HPV 16 E6, E7 and E6/E7 immortalized human keratinocytes, and CaCx cells (CaSki, SiHa, HeLa). e Analysis of Gal-7 expression in clinical samples (derived from HPV16 positive SCCs) and normal tissue. (^BKruskall-Wallis Test, P < 0.0001, Dunn’s Multiple Comparison Test P < 0.05; ^ET Test, P < 0.05, *** means highly significant, ** moderately significant)

Methods

Galectin expression profiles in independent cohorts of cervical cancer

To analyze the expression profiles of the galectin genes in cervical epithelium, we used the data from the differential expression profiles established by Scotto et al. [9] (Gene Expression Omnibus NCBI, GEO accession number GSE9750). The cohort consists of three groups of 24 normal cervical epithelium samples, a panel of nine CaCx cell lines and 28 squamous cervical cancer samples (SCC). Clinical samples represent the different cancer stages as suggested by the International Federation of Gynecology and Obstetrics (FIGO) [10]. The differential expression analyses were held using the reported signal intensity information for each galectin gene. Comparisons among groups were performed using a Kruskall-Wallis test followed by a post-hoc Dunn’s multiple comparison test. To test the Gal-7 profile through the CaCx progression, we used the cohort analysis established by Zhai et al. [11], (GEO accession number GDS3292) composed of 10 normal cervix, 7 high-grade squamous intraepithelial lesion (HSIL) samples and 21 SCCs. The Zhai cohort was used to study the modulation of Gal-7 along the process of transformation towards malignancy.

Hierarchical clustering and survival analysis

To study the expression profiles of Gal-7 and Gal-1 genes, the data compiled in a Provisional cohort (November 2014) for Cervical Squamous Cell Carcinoma and Endocervical Adenocarcinoma from The Cancer Genome Atlas (CESC-TCGA) [12] were used. This consisted of 185 RNA-Seq samples from patients in different tumor stage and grading. RNA-Seq Expectation-Maximization (RSEM) [13] data normalized by TCGA were used directly. The comparisons among groups were performed using a Kruskall-Wallis test followed by a post-hoc Dunn’s multiple comparison tests.

Unsupervised Hierarchical clustering was performed on a centered Pearson correlation coefficient algorithm and a complete linkage method to compare the expression patterns of the Gal-7 and Gal-1 genes across the CESC-TCGA data set. Overall survival (OS) was defined from the day of the sample intake to the patient’s death. Data of the patients who had survived until the end of the observation period were censored at their last follow-up visit. The OS curve was plotted using the clusters obtained from the hierarchical clustering analysis, using the Kaplan-Meier method. A log-rank test was used to compare the survival curves. All the statistical analysis and graphics were performed using R environment, version 3.0.1 (2013-05-16) “Good Sport”.

CESC-TCGA data methylation analysis

The relationship between Gal-7 methylation and mRNA expression was analyzed using the CESC-TCGA Methylation data versus mRNA expression. To confirm the correlations, the Spearman’s rank correlation coefficient was used, with a two-tailed P-value and an alpha = 0.05.

Analysis of methylation

Bisulfite-PCR-pyrosequencing focused on 13 CpG sites between nt −294 and nt +132 of the Gal-7 gene (the target region encompasses the promoter, the first exon and the first intron of the gene; the Eukaryotic Promoter Database ID: LGALS7 (chr19:39,262,157–39,266,157) was employed. PCR primers were designed on the in silico-converted target sequence, with the correspondent sequencing primers, to generate two products (Additional file 1: Table S2). Total DNA (500 ng) extracted from cells underwent bisulfite conversion by using the EZ DNA Methylation^TM Kit (Zymo reasearch), following manufacturer’s instruction. Two PCR reactions were performed in a final volume of 50 μl with 200 mM dNTPs mix, 500 nM each primer, 2.5 mM MgCl₂, 0.62 U Hot Start GoTaq® DNA polymerase (Promega) and 5 μl of bisulfite-converted DNA, under the following conditions: 95 °C for 2 min, 50 cycles of 95 °C for 30 s, 56 °C for 30 s, 72 °C for 30 s and final extension at 72 °C for 10 min. For the sequencing reactions of the PCR products, the PyroMark™ Q24 System (Qiagen) was used following manufacturer’s instructions (400 nM of one sequencing primer/reaction). The final pyrograms were analyzed using the allele quantification (AQ) mode in the PyroMark™ Q24 Software, to determine the proportion of C/T (forward sequencing) and, hence, the methylation frequency of the target CpG sites.

Cell culture, azacytidine treatment

HPV16 immortalized human keratinocytes, HeLa, SiHa and CaSki Gal-7 positive and control cells as well as 293 T cells were maintained under standard conditions in Dulbecco’s modified Eagle medium supplemented with 10 % Fetal Calf Serum and 1 % penicillin/streptomycin (P/S). Azacytidine (5′-Aza; Cayman) was dissolved in DMSO (Merck). Treatment was done for four days using 1 (HeLa) to 10 (SiHa and CaSki) μmol/L.

Gal-7 cloning and retrovirus construction

The Human Gal-7 gene cloned in the pEF1 vector [14] was amplified by PCR using primers containing the restriction sites XhoI and BamHI: Gal-7-XhoI-F: 5′-GAGCTCGAGCCGCCATGTCCAACGTCCCCCACAA-3′ and Gal-7-BamHI-R: 5′-GCGCGGATCCTCAGAAGATCCTCACG-3′ under the following parameters: denaturation at 98 °C for 3 min, followed by 30 cycles at 94 °C 30 s, 58 °C 45 s, 72 °C 1 min and a final extension at 72 °C for 10 min. Subsequently, 2 μg of Gal-7 insert DNA was digested with XhoI and BamHI enzymes (Fermentas, Sankt Leon-Rot) for 1 h at 37 °C and purified using the QIAquick® PCR Purification Kit (Qiagen, Hilden). Afterwards, the digested product was ligated into pLXSN viral vector previously digested with XhoI and BamHI and dephosphorylated. Subcloning was further validated by sequencing of the resulting pLXSN-Gal-7 construct (GATC Biotech, Germany). For production of viral particles, 293 T cells were transfected with the construct pLXSN-Gal7 together with the corresponding retroviral packaging vectors pVSV.G [15] and pSVΨ [16]. Viral particles containing empty pLXSN vector were also produced as control. 293 T cell supernatant was harvested and filtered through a 0.45 μm filter (Minisart Plus, Sigma-Aldrich). Filtered supernatant was mixed with 5 μg/ml of the cationic polymer Polybrene® and used for infection of CaSki, HeLa and SiHa cells. After infection, the cells were subjected to antibiotic selection with G418 for two weeks (Life Technologies). Gal-7 expression was finally analyzed by Western blot.

Protein extraction and Western blotting

Gal-7+ cervical tumor cell lines and the respective control cells were collected, washed in 1 × PBS and resuspended in RIPA buffer (20 mM Tris pH 7.5; 150 mM NaCl; 1 mM Na₂EDTA; 1 mM EGTA; 1 % NP-40; 1 % sodium deoxycholate) containing 1× complete protease inhibitor cocktail (Roche Diagnostics, Mannheim). Protein extraction was performed by incubating for 30 min on ice and subsequently centrifuging for 30 min at 4 °C and 13,000 rpm. Supernatants were quantified using the Bio-Rad Protein Assay Dye Reagent Concentrate (Bio-Rad, Munich). 80 μg protein/lane was used for Western blotting. After transfer, the PDVF filters were incubated with the following antibodies: anti-human Galectin-7, epr4287 (Genetex, USA), anti-Tubulin (G8), sc-55,529 (Santa Cruz Biotechnology, USA). Blots were developed using the Western lightning plus-ECL system (Perkin Elmer, Rodgau).

Colony formation assay and mitochondrial membrane potential measurements

The colony formation assay was essentially performed as described by Rotem et al. [17]. HeLa cells were cultured at 500, CaSki and SiHa at 1×10³ cells in 6-well plates. Cells were cultured for 7 days or until overlapping colonies started to appear. Cells were then briefly washed with DPBS and fixed for 30 min with a 10 % glutaraldehyde solution. Colonies were stained with 6 % crystal violet for 30 min. Each well in the plate was photographed with a Stereomicroscope (Olympus, Hamburg) and the colonies were automatically counted using the image processor software ImageJ.

To measure mitochondrial membrane potential, 15,000 cells from each cell line were plated in 100 μl normal DMEM in 96-well μClear black plates. The next day, the cells were treated with different concentrations of HA14-1 (10 μM, 25 μM and 40 μM) for 24 h. Afterwards, 100 μL/well of 30nM JC-1 Solution (Life Technologies, Darmstadt) were added and the cells were incubated for 30 min at 37 °C in the dark. After washing the plates twice with 1× Dilution Buffer solution (PBS + 2 % FCS) fluorescence measurement was performed by reading the plate in a Synergy 2 Multi-Mode microplate reader with the following settings: excitation: 475 nm, emission at 530 (monomers emission) and 590 nm (aggregate emission). A decrease in the ratio between the fluorescence of monomers and aggregates indicates mitochondrial depolarization and cell death.

Animal tumor model

7–8 week-old female nude Balb/c mice (Janvier Labs, St Berthevin, France) were maintained under pathogen-free conditions. 5×10⁵ (HeLa and HeLa Gal-7) or 5×10⁶ (SiHa and SiHa Gal-7) cells were suspended in 100 μl or 200 μl ice-cold PBS and subcutaneously implanted into the flanks of mice (5 mice/group). Tumor sizes were measured with an electronic digital caliper (Farnell, Germany) three times a week and the tumor volume was calculated according to the formula: V = 1/2 × length × width² (mm³). Animals were killed according to the animal welfare act when a tumor volume reached a volume of 1500 mm³. The animals at the German Cancer Research Center were maintained in compliance with German and European statutes and all animal experiments were undertaken with the approval of the responsible Animal Ethics Committee.

RNA extraction, reverse transcription and quantitative PCR analysis

RNA was extracted from cells using the RNeasy Mini Kit (Qiagen) according to the manufacturer’s instructions. One μg of RNA was reverse transcribed using RevertAid Reverse Transcriptase (Thermo Scientific, USA) and dT22primers according to the manufacturer’s protocol. The resulting cDNA was used for quantitative PCR analyses using the CFX96 Touch™ Real-Time PCR Detection System (BioRad, USA), the iTaq Universal SYBR Green (BioRad, USA) and the primers described in Additional file 1: Table S3.

Gene expression profiling analysis

RNA was extracted from tissue culture cells using the RNeasy Mini Kit (Qiagen) according to the manufacturer’s instructions. A total of 500 ng of RNA from every sample was labeled and hybridized in the HumanHT-12v4 expression BeadChip (Illumina) following standard procedures of the Genomic Core facility at the DKFZ. For the tumors, fresh-frozen xenotransplants isolated from nude mice were homogenized in a Precellys® 24 tissue homogenizer (Bertin Tech. USA). Subsequently, RNA from tumors was isolated using standard TRIZOL procedure (Invitrogen) and Direct-zol RNA (Zymo Research) following the manufacturer’s instructions. Samples within the groups were pooled. A total of 500 ng of RNA from every pooled group was labeled and hybridized in the HumanHT-12v4 expression BeadChip (for human genes) or the mouseWG-6v2 expression BeadChip (for mouse genes) following standard procedures of the Genomic Core facility at the DKFZ. Bioinformatic data analysis was performed using Chipster [18] software version 3.1. mRNAs with a fold change of at least 2 between Gal-7 negative and Gal-7 positive cells or tumors were considered significant and were used for further analysis.

Cell lysis and protein digestion for mass spectrometry analysis

Cell pellets were resuspended in 1 % (w/v) sodium deoxycholate (SDC) and 50 mM ammonium bicarbonate, boiled at 99 °C for 5 min and cooled to 4 °C. The lysates were then treated with benzonase (Sigma-Aldrich, Germany) in a final concentration of 150 U/mL for 60 min on ice. Insoluble material was removed by centrifugation at 14,000 G, 15 min, 4 °C. Protein concentration was determined by bicinchoninic acid protein assay kit (Sigma-Aldrich, Germany). Corresponding light and heavy lysates were mixed in 1:1 protein ratio. Lysates were then reduced with 10 mM dithiothreitol at 37 °C for 60 min, alkylated with 20 mM iodoacetamide at RT for 30 min in the dark and the unreacted iodoacetamide was quenched with further 10 mM dithiothreitol at RT for 15 min. The samples were diluted with 50 mM ammonium bicarbonate to decrease the concentration of SDC to 0.5 % and digested with sequencing grade trypsin (Promega, USA) overnight at 37 C. SDC was removed by the modified phase transfer protocol [19]. Briefly, ethyl-acetate was added and the digested product was acidified by trifluoroacetic acid (TFA) to a final concentration of ca. 2 % (v/v). The mixtures were vortexed vigorously for 1 min, centrifuged at 14,000 G for 5 min and the upper organic layer was removed. The extraction was repeated with fresh portion of ethyl-acetate. The aqueous phases were desalted on Empore™ extraction cartridges (Sigma-Aldrich, Germany) and dried in vacuum.

SILAC labeling

HeLa and SiHa cells as well as their Gal-7 reconstituted counterparts were SILAC labeled by cultivating the cells for 5 passages in DMEM media (without glutamine, arginine and lysine; Silantes 282,986,444) containing 10 % (v/v) dialyzed FBS. Isotopically labeled L-lysine [¹³C₆ ¹⁵N₂ - labeled] 0.798 mM and L-arginine [¹³C₆ ¹⁵N₄ - labeled] 0.398 mM (Silantes, Germany) were added to the DMEM media (Fig. 4b). Unlabeled L-proline was added to a final concentration of 2.61 mM (300 mg/L) to prevent arginine-proline conversion [20].

Peptide separation and mass spectrometry analysis

Peptides were separated by two-dimensional liquid chromatography (LC). In the first dimension, reversed phase LC under high-pH mobile phase conditions was performed using Alliance 2695 LC system (Waters, UK). The mobile phases were (A) water, (B) acetonitrile (ACN) and (C) 200 mM ammonium formate pH 10. Dried peptide mixtures were dissolved in 25 % C and 4 % ACN, and an aliquot of 200 μg was loaded on trap column (Gemini C18, 2 × 4 mm) and column (Gemini C18, 3 μm, 110 Å, 2 × 150 mm; Phenomenex, USA). Peptide separation was performed by linear gradient from 5 to 55 % of B in 62 min with constant 10 % of C. Flow rate was 0.16 mL/min, the column was kept at 40 °C and the separation was monitored at 215 nm. Fractions were collected manually in 2-min intervals over the sample elution window from 10th to 52nd min, acidified with TFA and dried in vacuum. The last 2 fractions were combined with the first 2 fractions in sequential order. Second dimension of the separation was performed on an Ultimate 3000 RSLCnano system (Dionex, USA) coupled on-line through Nanospray Flex ion source with Q-Exactive mass spectrometer (Thermo Scientific, Germany). Fractions were dissolved in 2 % ACN/0.05 % TFA and loaded on capillary trap column (C18 PepMap100, 3 μm, 100 Å, 0.075 × 20 mm; Dionex) by 5 μL/min of 2 % ACN/0.05 % TFA for 5 min. Then they were separated on capillary column (C18 PepMap RSLC, 2 μm, 100 Å, 0.075 × 150 mm; Dionex) by step linear gradient of mobile phase B (80 % ACN/0.1 % FA) over mobile phase A (0.1 % FA) from 4 to 34 % B in 48 min and from 34 to 55 % B in 10 min at flow rate of 300 nL/min. The column was kept at 40 °C and the eluent was monitored at 215 nm. Spraying voltage was 1.75 kV and heated capillary temperature was 275 °C. The mass spectrometer operated in the positive ion mode performing survey MS (at 350–1650 m/z) and data-dependent MS/MS scans on 10 most intense precursors with dynamic exclusion window of 30 s. MS scans were acquired with the resolution of 70,000 from 10⁶ accumulated charges; maximum fill time was 100 ms. The intensity threshold for triggering MS/MS was set at 5×10⁴ for ions with z ≥ 2 and the isolation window was 1.6 Da. Normalized collision energy for HCD fragmentation was 27 units. MS/MS spectra were acquired with the resolution of 17,500 from 10⁵ accumulated charges; maximum fill time was 100 ms.

Protein identification and quantification

Database search and quantification were performed by Proteome Discoverer v.1.4 software (Thermo Scientific). The reference proteome set of Homo sapiens containing canonic and isoform sequences was downloaded from UniProt [21] (http://www.uniprot.org/) on Aug 18th 2014 and merged with the common contaminants file downloaded from the MaxQuant web page (http://www.coxdocs.org/doku.php?id=maxquant:common:download_and_installation); the merged database contained 89,252 sequences. The search parameters were as follows: digestion with trypsin, max. 2 missed cleavages, allowed peptide mass tolerance of 10 ppm, fragment mass tolerance of 0.02 Da, fixed carbamidomethylation of cysteine, variable modifications: oxidation of methionine, acetylation of protein N-term and SILAC labels Arg10 a Lys8. The strict target value of FDR for a decoy database search of 0.01 was applied (high confidence). Only unique peptides were considered for quantification and the heavy to light ratios were normalized on protein median for each replicate (Additional file 2: Supplementary File 1).

For relative protein quantification only protein groups with a minimum of 2 identified peptides in all three replicates and a minimum of 1 quantified peptide per each replicate were considered. Log2 values of protein ratios from each replicate were then subjected to the ranking test to find the most significantly regulated proteins [22] taking the protein groups found in top (T) or bottom (B) groups in all three replicates (TTT or BBB). The false discovery rate (FDR) was evaluated by non-parametric estimate as an average number of proteins in the “false” groups (TTB, TBT, BTT, BBT, BTB and TBB). The significance cut-off of the FDR was set around 5 %. For analysis, only proteins with a log2 fold change value of 1.2 or higher were considered (Additional file 3: Supplementary File 2). The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium [23] via the PRIDE partner repository with the dataset identifier PXD001806.

Pathway and GO enrichment analysis

We performed an enrichment analysis of pathway-based sets of proteins considering all the nodes of our extended network. Enrichment was done employing ConsensusPathDB, of the Max Planck Institute for Molecular Genetics, by using the overrepresentation analysis online tool. As input, we uploaded the UniProt protein identifiers of all the elements of the extended network. We searched against pathways as defined by Reactome [24] and KEGG [25], with a minimal overlap with the input list of 5 and a p-value cutoff of 0.001.

Also, employing the same website and the same analysis tool, we performed an enrichment analysis based on Gene Ontology (GO) [26] level 3 category of biological processes. For this analysis, we considered only the identified core proteins and set the p-value cutoff on 0.001.

Network construction

Network reconstruction was performed with the aid of the Cytoscape Plugin, BisoGenet [27], using the identified proteins as bait nodes and adding edges with the following parameters: Organism > Homo sapiens, protein identifiers only; Data Settings > protein-protein interactions; all data sources and all experimental methods; method > By adding edges connecting input nodes and as Output > Proteins.

The iRegulon plugin was used to predict their transcriptional regulators using the default setting. Only predicted transcriptional regulators with normalized enrichment scores (NES) >3 were used.

Results

Gal-7 is downregulated in squamous cervical cancer, high-grade squamous intraepithelial lesions and cervical cancer cell lines

To get an insight into the galectin status of high-risk human-papillomavirus (HPV)-induced tumors, we first examined the expression of nine different members of this gene family on microarrays from the Scotto cohort [9] (Additional file 4: Figure S1). A significant down-regulation of Gal-7 expression was observed in fresh samples derived from squamous cell carcinoma (SCC), as well as established HPV16/18-positive CaCx cell lines (Fig. 1b). Considering the signal intensity of the microarrays, Gal-2 and Gal-3 transcription was also affected to some extent in SCC/CaCx samples compared to normal cervical tissue (Fig. 1b). However, in contrast to Gal-7, these two galectins are not directly related to tumorigenesis [28]. Only Gal-7 negative regulation showed a high statistical significance (Fig. 1b). This could be further supported by the Zhai cohort [11] (Additional file 4: Figure S1. B), where Gal-7 expression was also reduced in high-grade squamous intraepithelial lesions (HSIL) and in cervical cancer samples (Fig. 1c). Hence, the gradual downregulation of Gal-7 in premalignant lesions and a marked reduction in SCCs suggest that it might play a role in the development of cervical cancer. We validated our observations from the meta-analysis by qPCR and Western blot (Fig. 1d). For this purpose, we also included human keratinocytes that were separately immortalized by the E6-, E7- and E6/E7 oncoproteins of HPV 16 [29]. While CaSki, SiHa, and HeLa cells showed almost a complete suppression of their mRNA steady-state levels, Gal-7 was increased in immortalized cells when compared with primary keratinocytes (Fig. 1d). However, no significant differences among the same cells could be discerned at the protein level (Fig. 1d, below). Moreover, consistent with the microarray data (Fig. 1b), clinical specimens obtained from SCC patients also demonstrated a significant decrease (p < 0.05) of the Gal-7 mRNA in biopsies when compared to normal control samples (Fig. 1e). Altogether, these data imply that Gal-7 downregulation correlates with cervical cancer progression.

Mutually exclusive expression of Gal-7 and Gal-1 determines clinical outcome and overall survival

In order to confirm the biological significance of the negative regulation of Gal-7 in a clinical context, we analyzed a cohort of patients to determine the correlation between their survival rate and the absence or presence of Gal-7 and Gal-1. Since Gal-1 is considered as a pro-tumorigenic galectin [30], we anticipated a mutually exclusive expression with respect to Gal-7 as an indication of positive clinical outcome in cervical cancer. In order to prove this assumption, we used Illumina 450 k data from the Cervical Squamous Cell Carcinoma project of The Cancer Genome Atlas (CESC-TCGA; n = 185) [31] and performed an Unsupervised Hierarchical Clustering analysis (UHC) (Fig. 2a). Here, two significantly differentiated clusters (p = 0.0001) were obtained (Fig. 2b). Group A was integrated by 104 patients expressing high levels of Gal-7 and low levels of Gal-1. In contrast, group B revealed just the opposite composition (low Gal-7 and high Gal-1 expression, data derived from 81 patients). Figure 2e shows a representative immunohistochemical staining of Gal-1/Gal-7 in normal versus tumor tissue (obtained from the Human Protein Atlas) [32]. To address the clinical significance of the UHC analysis, we also assessed the overall survival in this cohort through a Kaplan-Meier survival analysis. Intriguingly, Group A had a significantly higher overall survival rate when compared to group B (p < 0.0001) (Fig. 2c). These results suggest that a mutually exclusive expression of Gal-7 and Gal-1 has beneficial prognostic value in CESC patients.

Fig. 2

Analysis of Gal-1/-7 expression and survival in the TCGA-CESC cohort. a Hierarchical clustering of Gal-7 and Gal-1 expression in a 185-patient panel of CESC from the TCGA. b Gal-7 and Gal-1 expression in TCGA-CESC patients. c Kaplan-Meier curve for 5000 days of overall survival in the CESC panel of TCGA. Censored events were marked with vertical black lines. d Inverse correlation of Gal-7 expression and methylation. e Immunohistological staining of Gal-1 and Gal-7 in normal and cervical cancer tissue sections (immunohistochemistry images were taken from the Human Protein Atlas Project). Gal-1 was detected using the HPA000646 antibody. Gal-7 was detected using the HPA001549 antibody. (^BMann Whitney test P < 0.0001 Gaussian Approximation two-tailed P-value; *** means highly significant. ^Clog-rank test; p = 0.005. ^DP < 0.0001, Spearman r, two tailed P value, alpha = 0.05)

Tumors are known to downregulate incompatible genes through epigenetic mechanisms, such as gene methylation [33]. Therefore, we analyzed if this was the mechanism behind Gal-7 repression in the cohort. We found a high correlation (p = 0.001) between the clinical outcome and the degree of expression of the Gal-7 gene, which in turn inversely correlated with methylation (Fig. 2d) between positions −499 to +100 relative to the initiation site (Additional file 1: Table S1). Our data suggest that Gal-7 negative regulation is a biological phenomenon with a strong impact in the outcome of cervical cancer patients.

Hypermethylation is responsible for reduced Gal-7 expression in CaCx cells

As is indicated by the meta-analysis, we examined whether the Gal-7 gene became de novo methylated during the multi-step progression to cervical cancer. For this purpose, cells were treated with 5′-azacytidine as a demethylating agent. Re-expression of the Gal-7 gene could be observed after demethylation (Fig. 3a). Protein recovery was not as strong as in HPV16-immortalized keratinocytes (Fig. 3a, right panel), implying that additional mechanisms control endogenous Gal-7 expression or protein stability in CaCx cells [34].

Fig. 3

Methylation analysis and re-expression of Gal-7 in CaCx cell lines. a qPCR and Western blot of Gal-7 four days after 5′-Azacyditine (Aza) treatment. Actin was used as the loading and transfer control. b Above: schematic overview of the Gal-7 gene: the locations of CpG sites, the transcriptional start site (TSS) and the first exon are indicated. Below: percentage of methylation after bisulfite pyrosequencing analysis. c Western blot after reconstitution of Gal-7 in CaCx cells and mock infected control cells. Tubulin was used as the loading and transfer control. d Colony formation assay of Gal-7+ CaCx cells (Gal-7+, red) versus mock-infected cells (Gal-7-, black). e Induction of apoptosis due to loss of the mitochondrial transmembrane potential after incubating Gal-7+ and Gal-7- CaCx cells with increasing amounts of H14-1. (^CKruskall-Wallis P < 0.0001, Dunn’s Multiple Comparison Test P < 0.05; ^{E, F G} T Test, P < 0.05, ** means moderated significant)

In silico analysis using the Eukaryotic Promoter Database [35] revealed several CpG sites at the 5′end as well as inside the gene (Fig. 3b). After bisulfite-sequencing, a methylation pattern from positions −300 to −12 was almost identical in all cell lines (Fig. 3b). Cervical carcinoma cells revealed strong hypermethylation of CpG sites localized within the first intron (at positions +2 to +132), whereas the same region was almost methylation-free in HPV16 Gal-7+ immortalized keratinocytes. These data suggest a gradual de novo methylation of Gal-7 during HPV-induced carcinogenesis.

Reconstitution of Gal-7 in CaCx cells conferred sensitivity to apoptosis and anchorage-independence

We next examined the phenotypic changes of CaCx cells after retrovirus-mediated Gal-7 reconstitution in vitro (Fig. 3c). Low-attachment and anchorage-independence are well-known and stringent in vitro parameters for transformation [17]. In a tumor formation assay in vitro, colony formation capacity was found to be significantly lower in Gal-7+ cells than in mock-transduced controls (Fig. 3d). Moreover, re-expression also enhances the susceptibility to apoptotic stimuli, which is consistent with the pro-apoptotic function of Gal-7 in keratinocytes [36]. This was further confirmed by treating cells with the specific Bcl-2 inhibitor and apoptosis inductor HA14-1 [37]. Subsequent staining with the lipophilic cationic dye JC-1 revealed lower cell viability as determined by the loss of the mitochondrial transmembrane potential (ΔΨm) (Fig. 3e). Our data shows that Gal-7-expressing cells are more susceptible not only to intrinsic, but also to extrinsic apoptotic signals [14]. We show that the expression of Gal-7 decreases cell viability and induces a apoptotic response in transfected cells.

Impact of Gal-7 re-expression on cellular networks of CaCx cells

Having shown that Gal-7 re-expression negatively affects cell growth and impairs colony formation by means of apoptotic signals, we hypothesized that there are deep changes in the cellular networks in response to Gal-7 re-introduction. In order to analyze this question at system-level, we combined the results from microarray expression profiling (Fig. 4a) with a SILAC-based proteomic approach (Fig. 4b). In total, 213 candidate genes were identified in vitro as being differentially expressed either at the RNA or protein level (38 in HeLa Gal-7+ and 185 in SiHa Gal-7+ cells versus mock-transduced Gal-7- control cells). Surprisingly, only three genes were differentially regulated in HeLa Gal-7+ cells at the transcriptome level (Fig. 4c). The proteomic analysis of HeLa/HeLa Gal-7+ revealed 35 proteins that were differentially expressed in three biological replicates (FDR 5 %, Fig. 4d). In the case of SiHa/SiHa Gal-7+ cells, 60 mRNAs were found to be differentially transcribed (Fig. 4e). The proteomic analysis showed 125 differentially regulated proteins (FDR 5 %, Fig. 4f). The difference between transcriptome and proteome suggests that there are post-transcriptional regulation mechanisms affecting protein expression levels.

Fig. 4

Changes in gene and protein expression in Gal-7+ CaCx cell lines. a Experimental strategy for differential gene expression analysis in Gal-7+ CaCx cell lines. b SILAC proteome experiments for the profiling of Gal-7+ cell lines. c-f Differentially expressed transcripts in HeLa (c) and SiHa (e) Gal-7+ cell lines in vitro with respect to Gal-7- controls. The X axis is the P value in log10 scale; the Y axis is the fold change. Differentially expressed proteins between HeLa (d) and SiHa (f) Gal-7+ cell lines in vitro are shown aside. The X axis represents the score in log10 scale (the identification score was obtained from the average score reported by triplicate identification experiments) and the Y axis is the fold change (FC). g, h Functional modules affected in Gal-7+ HeLa (g) and SiHa (h) cells (red: up-regulation, blue: down-regulation)

To study the biological processes in which the identified molecules are involved, we performed a functional analysis of proteins and transcripts using levels 3 and 4 of the “Biological processes” Gene Ontology (GO) [38] through which 120 GOs were obtained (q-value = 0.001). We condensed the redundant GOs to obtain “functional modules” by pooling the GO domains according to their participation in pathways and cancer hallmarks [39]. As summarized in Fig. 4g and h (Additional file 1: Table S4 and S7), HeLa/HeLa Gal-7+ and SiHa/SiHa Gal-7+ shared eight functional modules, but their composition and the extent of regulation between them was different. Moreover, two additional modules were recognized in SiHa, namely “Signal transduction” and “Post-translational modifications” (Fig. 4h). Notably, although we studied the same cancer entity, our proteo-transcriptomic data show that re-expression of Gal-7 can trigger different cell-context dependent responses but leads to a convergent phenotype, as reflected in specific biological processes and cancer hallmarks.

Mouse microenvironment pressure exerts differential gene expression in Gal-7+ tumors

To investigate if Gal-7 reconstitution impairs the tumorigenicity of xenografts, HeLa and SiHa Gal-7 + cells, as well as their mock Gal-7- control cells, were subcutaneously injected into athymic nude mice and tumor growth was monitored. The tumors were excised when they reached a volume between 800–1000 mm³ and were used for transcriptome analysis (Fig. 5a). Consistent with the lower colony formation capacity and higher sensitivity to apoptotic signals under in vitro conditions, both HeLa Gal-7+ (Fig. 5b, d) and SiHa Gal-7+ cells (Fig. 5c, e) showed a significant impairment in their tumor formation capacity.

Fig. 5

The microenvironment induces changes in Gal-7+ tumors. a Experimental strategy. b and c Growth of Gal-7+ and Gal-7- Hela and SiHa tumors. d and e Incidence of tumors in injected mice. f, g Gene profiling analyses in Gal-7+ and Gal-7- HeLa and SiHa tumors. X axis: P value in log10 scale, Y axis: fold change (FC). h, i Functional modules and differentially expressed transcripts in Gal-7+ HeLa. h and SiHa (i) tumors in vivo. (^{B C D E} ***p < 0.001; **p < 0.01; *p < 0.05 Student’s t test). (red: up-regulation, blue: down-regulation)

To understand this phenomenon at the molecular level, we performed microarray gene profiling to monitor changes that could be linked to decreased tumorigenicity and animal survival. A total of 42 differentially expressed genes were found in HeLa Gal-7+ cells (Fig. 5f) that could be divided in 8 functional modules (Fig. 5h, Additional file 1: Table S4). In SiHa Gal-7+ cells, 58 genes were differentially regulated (Fig. 5g), and these could be split into 10 functional modules (Fig. 5i, Additional file 1: Table S7), where only PLAC8 and TCN1 were in common. On the other hand, Gal-7 re-expression led to the identification of 8 genes that were regulated in both CaCx cell lines, either at transcriptional or translational level (OASL, COX19, FAM129A/NIBAN, CDKN2A/p14arf, CRYAB, BCL3, FST and TACSTD2) (Fig. 6b). We validated their expression by qPCR, and observed whether the same genes were affected when Gal-7 reconstituted HPV16-positive CaSki cells were examined (Additional file 5: Figure S2). Our data indicate that Gal-7 has a general impact on these network components, which are shared among the three different cervical carcinoma cell lines. These results highlight the complex influence of the microenvironment on the modulation of gene expression in comparison to the in vitro growth.

Fig. 6

Differential network analysis of Gal-7+ CaCx cells and tumors. a Integrative analysis strategy. b Common candidates that were differentially regulated at the transcriptional or translational level compared against their Gal-7- controls. Numbers indicate fold changes. c, d Differential network analysis in Gal-7+ HeLa and SiHa cells and tumors. Color depth represents the fold change with respect to the control. Red: up-regulation, blue: down-regulation. Major transcription factors (MTF) are shown in the center. The colored lines derived from MTFs show the transcriptional relationship. The shape of the symbols indicates the following: square for transcript, circle for protein, hexagon for TF, diamond for miRNA, and triangle for lncRNA. Wide, green shape edges represent that the information for this molecule was obtained in vitro, while the thin, black edges mean in vivo molecule information

Living systems are dynamic and complex. Their behavior may be hard to predict from the properties of individual parts, since their interactions produce emergent properties [40, 41]. To get a more accurate insight into the emergent properties of the Gal-7-dependent regulation, we performed an integrative analysis using the transcriptome and proteome biological layers (Fig. 6a). We used the Cytoscape software [42] to build a network for each cell line, and we made a prediction of transcriptional factors (TFs) with iRegulon [43]. The entities were grouped according to their participation in each synthetic module. In HeLa Gal-7+ cells (Fig. 6c), five major TFs are able to control about 60 % of all the genes in the network; these are SP9, ARID3A, MYF5, IKZF2 and SCRT2 (Additional file 1: Table S5). When the function of each individual module was dissected, we found an up-regulation of several candidates involved in the impairment of growth and proliferation (cell cycle and proliferation/regulation of programmed cell death modules), but also in cell adhesion, migration, metabolism, tissue organization and integrity.

The integrative analysis of SiHa Gal-7+ cells (Fig. 6d) showed 6 major TFs that are able to regulate up to 65 % of the genes within this network: E2F1, E2F6, EP300, SMAD1, NFIC and FOXJ3 (Additional file 1: Table S6). Despite having different members, we found a regulation of the same modules that were present in HeLa Gal-7. This could account for the similar phenotypes showed by our xenotransplantation analysis. Notably, there were also changes in genes involved in post-translational modifications and signaling transduction, two modules that were not identified in the Hela Gal-7 network (Additional file 1: Table S7).

Gal-7 re-expression in tumors modulates transcriptional changes in their associated microenvironment

Xenografting of human cancer cells into mice gives us a unique opportunity to study the crosstalk between the tumors and their microenvironment. To follow that purpose, we analyzed how Gal-7 reconstitution shapes the behavior of the microenvironment and the infiltration of immune cells into tumor tissues. To selectively obtain the gene expression profile of infiltrates (e.g. tumor-associated-epithelial, stromal, vascular and immune cells), we followed the competitive cross-species hybridization of microarray experiment (CHME) strategy [44] of the xenografts (Fig. 7a). In the HeLa-Gal-7+ Tumor associated microenvironment (TAM), a total of 87 differentially regulated mouse genes were revealed (Fig. 7b). Subsequent functional module analysis assembled these genes into 9 modules (Fig. 7d), where “tissue morphogenesis”, “metabolism and transport”, and “chemokine activity and immune response” comprised most of the differentially expressed genes. A module related to “hypoxia response” appears only in this TAM. In the SiHa Gal-7+ TAM, 91 differentially regulated genes were found (Fig. 7c). 11 functional modules could be clustered (Fig. 7e), mostly showing a strong down-regulation in their differentially expressed genes. Notably, while in the HeLa-Gal-7+ TAM the genes related to “chemokine activity and immune response” and “proteases” were over-regulated, in the corresponding modules of the SiHa-Gal-7+ TAM all the genes were down-regulated. Also remarkable is the inclusion of two additional modules in the SiHa-Gal-7+ microenvironment: “calcium signaling” and “degradation of macromolecules”. These results suggest that the mechanisms that determine the growth impairment phenotype in response to Gal-7 reconstitution depend on the cellular context, supporting the notion of the uniqueness of cancer cells.

Fig. 7

Differential gene profiling of the Gal-7+ tumor associated microenvironment. a Experimental strategy. b, c Differential gene profiles of the of Gal-7+ HeLa and SiHa TAM with respect to controls. X axis: P value in log10 scale; Y axis: fold change (FC). d, e Modules and differentially expressed transcripts in Gal-7+ HeLa (d) and SiHa (e) TAM (up-regulation in red, down-regulation in blue)

We made a network and TF analysis of the TAM’s of the Gal-7+ tumors. In the HeLa-Gal7+ TAM, regulation of all 9 modules could be attributed mainly to three mouse TFs: Arfgap1, Egr1 and Srf (Fig. 8a and Additional file 1: Table S8). Interestingly, in the “Regulation of transcription” module we found the down-regulation of the Abra gene, which activates Srf and Egr1 [45, 46]. This confirmed our TF prediction and could provide an explanation for the generalized down-regulation of genes. In the “chemokine activity and immune response” module, an up-regulation of several genes (Additional file 1: Table S9) could be discerned, arguing for a wide immune host response against the Gal-7 expressing tumor.

Fig. 8

Differential network analysis of the Gal-7+ tumor associated microenvironment. Differential network analysis in Gal-7+ HeLa (a) and SiHa (b) TAM. Color depth represents the fold change. MTFs are shown in the center. The colored lines derived from MTFs show the transcriptional relationship. The shape of the symbols indicates the following: square for transcript, hexagon for TF, and triangle for lncRNA (up-regulation in red, down-regulation in blue)

To further support the potential immune mechanism in this system, differentially expressed genes were also analyzed using the immunological genome project (IGP) database [47]. The HeLa Gal-7+ mouse microenvironment was enriched with the expression of Stfa1, Stfa2, Retnlg, S100a8 and S100a9. The two latter marker genes are constitutively expressed by neutrophils, activated monocytes and macrophages. We also identified IL1β, CXCL2, CCL3, CCL4 and Arg1, characteristic genes expressed in macrophages, NK cells, and activated polymorph-nuclear neutrophils, respectively. Moreover, Galnt9 is present only in macrophages and monocytes (Additional file 6: Figure S3B), while Zmynd19 expression is found within hematopoietic stem and multi-lineage progenitor cells (Additional file 6: Figure S3A). The SiHa Gal-7+ microenvironment network analysis showed strong repression in all the functional modules. The major TFs that drove the network were Trp53, Nr2c2 and Mecom (Fig. 8b and Additional file 1: Table S10). A detailed description of some important genes is provided in Additional file 1: Table S11.

Using the IGP database, we mapped the up-regulated genes in the SiHa Gal-7+ mouse microenvironment. Here, we found that the expression of Taf15 (TATA Box Binding Protein (TBP)-Associated Factor), amongst others, was restricted to B-cells, dendritic and NK cells, but not to stromal cells. Whether SiHa Gal-7+ cells induce an environment that suppresses the infiltration of tumor-associated stromal and immune cells awaits further elucidation. Altogether, these data suggest a possible role of the innate cytotoxicity and an adaptive antibody and cytokine secretion mechanism in the growth impairment of SiHa Gal-7+ tumors.

Discussion

Cervical cancer is a highly robust and flexible network [48] that is shaped by particular microenvironments [49]. Consequently, the combination of clinical data and integrative system-level approaches [50–53] is necessary to tackle the inherent complexity of the disease [54, 55]. In fact, the network of CaCx cells is deeply influenced by the expression of the E6 and E7 high-risk-HPV oncoproteins, which always attack hubs to get selective advantage [56, 57]. These events occur not only via interfering with the p53 and the retinoblastoma (pRB) pathways [58], but also through a direct evasion of anti-tumorigenic responses by 1) impairing central immunomodulatory molecules including type I interferons [59], IL-1β [60], and chemokines [61], and 2) acquiring resistance to immune (CTL and NK)-mediated apoptosis [62, 63].

Gal-7 is incompatible with the malignant phenotype in cervical cancer

Gal-7 possesses well known functions in pro-apoptotic response [7, 64, 65], regulation of cellular proliferation [66], and epithelial homeostasis [67]. In our study, several lines of evidence pinpointed Gal-7 as incompatible with the selection process during cervical carcinogenesis. Such phenomenon is clearly visible in the negative regulation of the gene at the transcript and protein level, in clinical samples as well as in cell lines. As shown by bisulfite sequencing in our cell lines and meta-analysis of methylation arrays from clinical samples, the absence of Gal-7 could be attributed to an intragenic de novo methylation event, a mode of gene silencing also observed in gastric cancer [68]. Gal-7 expression, however, could be restored by treating cell lines with a demethylating agent (Fig. 3a).

As depicted in our results, Gal-7 reconstitution does not lead to the suppression of the tumorigenesis, but rather results in a significant reduction of colony formation capacity and an increased sensitivity towards apoptosis in vitro, consistent with observations in prostate cancer cells [65] and in colon carcinoma cells [69]. Xenotransplantation in nu/nu mice also showed a strong tumor growth retardation, a phenotype in common with Gal-7+ colon carcinoma cells in SCID mice [69]. Interestingly, other cancer models display an opposite behavior. Gal-7 is over-expressed in lymphoma cells and is related to promoting tumorigenesis by induction of MMP9 [70] and in ovarian cancer it is associated with cell proliferation and poor prognosis [71].

Gal-7 reconstitution triggers changes in proteins and mRNA in the cell lines and tumor

Our results demonstrated that re-expression of Gal-7 in the HeLa and SiHa cell lines triggers a cascade of changes at both transcript and protein levels. In vitro, HeLa Gal-7+ only yielded 3 differentially expressed genes and 35 proteins, while in SiHa Gal-7+ we found 60 differentially expressed transcripts and 125 proteins. Several studies show significant differences between the abundance ratio of the messenger RNA transcript and the corresponding protein product [72, 73], a phenomenon that could be mainly due to post-transcriptional and translational regulation, as well as protein and mRNA turnover [74].

Within tumors, a continuous selection exists for cells with the highest survival and proliferative advantage; moreover, microenvironmental conditions provide selective pressure for tumor evolution and the progression to invasion [75]. One of the most probable mechanisms for inhibition of the tumor growth could be the expression of genes with “asymmetric activity” capable of suppressing tumorigenicity, without affecting in vitro growth [76].

Three central points can be extracted from the present work, the first of which is related to the robustness of the cancer network. Despite the gene composition, there is a powerful convergence in the functional modules towards the hallmarks of cancer, as suggested by our previous work [48]. A factor to be discussed is if these functional modules are a response to the imminent cell destruction, a desperate response to survive, or a combination of both. An example in HeLa Gal-7+ cells is the case of over-expression in vitro (protein) and in vivo (transcript) of SERPIN3B, a protease inhibitor that might act in protecting cancer cells through inhibition of apoptotic signals [77, 78] and is capable of inducing epithelial-mesenchymal transition [79]. Another one is the over-expression of GAGE genes in SiHa Gal-7+ cells, whereas their expression in normal tissues is limited to the germ cells of immune-privileged organs [80]. It is possible to speculate whether this response is an attempt to escape the immune modulation.

System-level response in Gal-7 network is context dependent

We selected the SiHa and HeLa cell lines because they are transformed by the most frequent high-risk HPV types (HPV16 and HPV 18, respectively). The E6/E7 oncoproteins of both HPV types possess subtle differences in their interaction partners [81], in the affinity of these interactions [82] and, consequently, in the mechanisms of transformation [83, 84]. These differences could potentially generate a chain reaction in the whole cell network, and each cell line could give a particular response due to its specific background. Although tumors of similar types and clinical outcomes can have patterns of gene expression or mutations that are strikingly different, these behaviors recurrently target the same pathway or network [55]. Which leads us to central point number two: while each cancer is as unique as its patient, there is convergence towards a very similar phenotype.

Although Gal-7 affects the specific cervical cancer networks in different ways, we identified a group of eight molecules that are regulated in both cell lines: CRYAB, TACSTD2, BCL-3, COX19, OASL, CDKN2A/p14ARF, NIBAN/FAM129A, and FST. All these proteins have been previously described as key regulators of tumorigenesis in some cancers (Table 1).

Table 1 Common proteins between Hela and SiHa Gal-7+ cell networks

In our differential network analysis, we predict the TFs that can potentially regulate the respective modules. The HeLa Gal-7+ network can be controlled by five major TFs: SP9, ARID3A, MYF5, IKZF2 and SCRT2, while in SiHa Gal-7+ cells, six major TFs could be identified: E2F1, E2F6, EP300, SMAD1, NFIC and FOXJ3 (Additional file 1: Tables S5 and S6). Surprisingly, we found a deep connection involving E2F TFs in both cell lines. In HeLa-Gal-7+, ARID3A, also known as E2FBP1, is a TF involved in the control of cell cycle progression by the RB1/E2F1 pathway, and is induced in response to p53 [85, 86]. In the SiHa-Gal-7+ network, E2F1, E2F6, and FOXJ3 are involved in pathways of cell cycle control, where FOXJ3 may be involved in regulating a network of zinc finger–binding proteins that may affect gene expression themselves in response to the E2F pathway [87, 88]. E2F TFs are part of a wide network involved in differentiation, apoptosis and immune response [89–91].

Gal-7 plays a role in the crosstalk between cancer cells and the tumor associated microenvironment

Since a tumor is always composed of heterogeneous cells in a context-dependent manner [92], expression profiling of infiltrates provides an insight into the bi-directional communication between transformed cells and the microenvironment. Our CHME strategy allowed us to monitor the expression profile of genes that were indicative of the infiltration of potential tumor-associated epithelial, stromal and immune cells [93, 94].

The integrative analysis of Gal-7+ TAM cells revealed deep changes, presumably activated by the interaction of the tumor and microenvironment networks. In the HeLa Gal-7+ TAM, the chemokine and immune responses stood heavily stimulated. Meanwhile, 30 genes related to tissue morphogenesis, metabolism and transport, and hypoxia response were down-regulated. Excitingly, we found the interplay between Abra (STARS) and its targets, Srf and Egr1. These three TFs have an important role in the remodeling, homeostasis and metabolism of skeletal muscle [95], and appear to be under control of AKT [96]. In the SiHa Gal-7+ TAM, there is not an obvious interplay between the predicted and regulated TFs. Nevertheless, we observed the negative regulation of almost 97 % of all the modules, and specifically of some genes that promote carcinogenesis, such as CCL5, CXCL9 and CXCL10 (see Additional file 1: Table S11). These genes are characteristic of tumor-associated macrophages and stromal fibroblasts.

It is noteworthy that the role of the immune system is highlighted by the expression of specific genes of NK cells in the HeLa Gal-7+ TAM. Meanwhile, in SiHa Gal-7+ TAM, we found evidences of activated monocytes and the infiltration of macrophages. This would indicate that the innate and adaptive immune system, stroma cells, and fibroblasts play a complex leading role in impairing the growth of Gal-7+ tumors.

These observations sustain the third central point of our work: there is an aggressive negotiation between the tumor and the TAM, where constant changes in both networks shape their overall behavior and fate. We summarize our three central points in Fig. 9. Supporting our empirical observations, formal mathematical modeling demonstrated that tumor mass evolution is the integrated result of the dynamics of two linked complex systems, tumor cell population and tumor microenvironment [97], in addition to emergent properties caused by the genomic background [98] and clinically significant macro-environmental variables [99, 100].

Fig. 9

Bidirectional crosstalk between Gal-7+ tumors and their associated microenvironments. The expression of genes related to stromal cells, cancer associated fibroblasts, and immune cells results in an aggressive crosstalk that could potentially decide the fate of the tumor: proliferation and immune evasion, or immune destruction. The role of the common proteins in the Gal-7 network remains partially unknown. The color of the participant transcripts and proteins represents differential expression (up-regulation in red, down-regulation in blue)

Numerous studies reinforce the hypothesis of a complex control by the stem cells over their niche, allowing for immunologic privilege. This effect allows them to exert their multipotentiality and proliferate at a different rate than the rest of the tissue [85–87]. Excitingly, it is known that secreted galectins are able to exert immunomodulation in the stem cell niche [88]. For instance, Gal-9 has recently been shown to prevent allo-reactivity of immune cells in a non-cancerous model [89], and has been shown to be a major mediator of the anti-proliferative and functional effects of multipotent mesenchymal stromal cells on T and B cells [90]. Moreover, Gal-3 is over-secreted in pancreatic cancer stem-like cells [91]. The interaction with stem cells could explain why Gal-1 and Gal-3 were found to drive a strong immunomodulation in cancer in a previous work [1].

Based on the previous evidence, we speculate that galectins (including Gal-1 and Gal-7) could act like molecular masterminds behind the regulation of immune activity in different cancers, particularly through interactions with stem cells. A delicate equilibrium in galectin expression could regulate immunosurveillance, while its disturbance by tumors and cancer stem cells would lead to immune escape, tumor growth and metastasis. In our Gal-7+ xenografts, the re-expression of Gal-7 could have either restored immunosurveillance or counteracted any immune escape mechanisms, which would explain the immune signature found in the TAM. Future studies should confirm if Gal-7 plays a direct role in the recruitment of immune cells, as well as whether they interact in any way with cancer stem cells to achieve such results.

Conclusions

Gal-7 expression presents an interesting correlation with a significantly better survival rate of patients diagnosed with cervical cancer. This phenomenon highlights the biological importance of Gal-7 repression in CaCx.

The cellular networks reconstructed by Gal-7 expression prove three central hypotheses implicated in the emergency of properties of the Gal-7 network: 1) Despite the gene composition, there is a convergence in functional modules that potentially govern the retardation of growth in CaCx Gal-7+ cells. 2) Regardless of their anatomical and cellular origin, each cancer is as unique as its patient. 3) A tumor-microenvironment cross-talk in response to Gal-7 reconstitution leads to an alteration of regulation and interaction networks, providing a great insight into molecular basis of the Gal-7 network.

Integrating several biological layers of information and comparing the variability of the dynamics and the complexity of the networks, using tools of the emergent field of differential network biology, is essential to gain definitive insights into cancer biology at the molecular level.