Surgical Tumor Resection Deregulates Hallmarks of Cancer in Resected Tissue and the Surrounding Microenvironment Open Access

This article is featured in Selected Articles from This Issue, p. 513

Surgery remains the mainstay of cancer management with curative intent (1). However, there is also some evidence that there could be metastatic dissemination of cancer cells during surgical removal of tumors (2–11). Although several studies have investigated gene expression in various tissues after surgery, changes induced in the primary tumor during surgical removal remain inadequately characterized (12–20). Moreover, most literature on gene expression profiles in tumors is based either on pre-surgical biopsy samples or those obtained from the resected surgical specimen, but not from both. The transient nature of surgical stress and its individual components makes it difficult to study them using single time point (diagnostic biopsy or surgical specimen) tissue samples.

The chief events during surgical resection include stress due to mechanical tissue shearing, hypoxia due to blood vessel ligation, and the systemic effects of anesthesia. Surgery-induced hypoxia is primarily perfusion-limited (acute hypoxia) rather than diffusion-limited (chronic hypoxia followed by cyclical reoxygenation), which is more commonly observed during natural neoplastic growth in vivo (21). The latter has been relatively well studied (22), whereas the former has yet to be systematically investigated (23). Our group has previously reported two randomized trials (24, 25), which suggest that short-term peri-operative interventions could favorably impact long-term survival outcomes in patients with breast cancer. Therefore, it is of interest to characterize the gene expression changes resulting from surgical stress, which could be potentially modulated for improving patient outcomes (26).

In this study, we evaluated the effect of surgical stress and the resulting acute, progressively severe tissue hypoxia on global gene expression changes using intraoperatively sampled biopsies from patients with breast cancer and head and neck squamous carcinoma (HNSC). We subjected RNA from breast cancer and HNSC patient samples obtained at different time points during surgery, to whole-transcriptome sequencing using next-generation sequencing (NGS), and gene expression microarrays, respectively. We studied the transient gene expression changes and biological pathways that were deregulated during surgical intervention.

All the studies from which patient samples were accrued for this analysis were conducted in accordance with the Declaration of Helsinki. All patients were recruited after obtaining written informed consent on The Tata Memorial Center (TMC) Institutional Ethics Committees (IEC)–approved peri-operative hypoxia studies 168 and 254 (ClinicalTrials.gov identifier NCT03797482) for breast cancer and 1259 for head and neck cancer, respectively. All patients were recruited in this study after screening as per the study eligibility criteria. Patients were explained the study protocol, their queries addressed, and subsequently, consented in writing to participate on the corresponding study. All Patients received standard cancer treatment as per institutional practice. They received general anesthesia (propofol, sevoflurane or isoflurane) at the time of surgery followed by standard curative surgical resection and procedures.

Patients with histologically confirmed, treatment-naïve, non-metastatic breast cancer with any receptor phenotype, who were scheduled for upfront curative surgery were eligible for this study. The following tissue samples were obtained during surgery.

• Time point A: Core needle biopsy of the tumor and adjacent normal tissue was performed at the start of surgery, immediately after anesthesia, and when the tumor surface was completely vascularized.

• Time point B: Core needle biopsy tissue was obtained midway during surgical resection, approximately 15–20 minutes after initiating surgery, and when the tumor was partially devascularized.

• Time point C: Tissue was obtained from surgical specimen at the end of surgery when the tumor was completely devascularized.

All tissue samples were preserved in RNAlater solution (Qiagen). RNA for breast and head and neck cancer samples was extracted using TRizol (Invitrogen), followed by purification using the PureLink RNA Mini Kit (Invitrogen). RNA quality control (QC) was performed using the Agilent 2100 Bioanalyzer (Agilent RNA 6000 Nano Kit), and total RNA concentration, RNA integrity number (RIN), 28S/18S ratio, and size were evaluated. The purity of the samples was tested using NanoDrop. Only samples with RIN ≥ 7.0, minimal r-RNA as evaluated by 28S/18S ratio and adequate quantity for library preparation were subjected to RNA-sequencing (RNA-Seq) library preparation.

Whole-transcriptome sequencing libraries were prepared using the manufacturer's (Illumina) recommended protocol (TruSeq RNA Sample Prep Kit v2) and sequenced on Illumina HiSeq 2000. This library preparation and sequencing methodology have been previously reported from our Centre in 31 patients with breast cancer (27). The NGS library preparation and sequencing methods are described in detail in Supplementary Materials and Methods.

Breast cancer RNA-Seq data were analyzed using two distinct and complementary methods. The first was an individual patient time-series analysis wherein sample at one time point was compared with the same patient's sample at another time point, that is, A versus B, B versus C, and A versus C. This enabled a within-patient comparison of gene expression at different time points during surgery. The gene expression changes that met the statistical criteria were selected for further analysis. This method allowed us to identify patient-specific gene expression changes (28–30). The secondary method was a replicate analysis, wherein all samples obtained at the same time point (A, B, or C) were treated as replicates and clubbed together, and differential gene expression was evaluated between time points A, B, and C. This analysis allowed us to identify gene expression changes that occur in the whole cohort with high statistical significance.

The tuxedo suite pipeline comprising the Bowtie (version 2.3, RRID:SCR_016368; ref. 31), TopHat (version 2.06, RRID:SCR_013035; ref. 32), Cufflinks (version 2.02, RRID:SCR_014597; ref. 33), and Cuffdiff (version 2.02, RRID:SCR_001647; ref. 34) was used for analysis as described by Trapnell and colleagues (35). Tophat aligner with default settings was used to align the RNA-Seq raw reads to the University of California Santa Cruz Biotechnology (UCSC) Homo sapiens reference genome (build hg19). A reference-based transcript assembly using Cufflinks was carried out using hg19-based GTF file containing all known transcripts (36) from the UCSC genome browser (RRID:SCR_005780; ref. 37) and was used for further downstream alignments with Cufflinks. The raw RNA-Seq reads for each sample (at time points A, B, and C) were run through this pipeline. For each patient, one aligned BAM file (Tophat output) and one assembled transcript GTF file (Cufflinks output) were generated for each sample A, B, and C. FPKM (Fragments Per Kilobase per Million mapped reads)-based normalization of global gene expression profile with respect to the differing library sizes and transcript length was performed, to generate a global normalized expression profile for each sample. Cuffcompare was used to merge the transcript assembly files for samples A, B, and C into a single GTF file for each patient, which was used for further downstream analysis. Cuffdiff was used to evaluate differential gene expression between the three samples (A vs. C, B vs. C, and A vs. C) of each patient in a time series experiment (35). Statistically significant differentially expressed genes were determined on the basis of the statistical significance criteria for each gene, as specified by Cuffdiff.

We applied a statistical filter using p < q value as a cutoff value for identifying significantly differentially regulated genes, as specified by Cuffdiff's significance criteria. We only analyzed genes with a p value of <0.05 by comparing the p value with the corresponding q value (Benjamini–Hochberg corrected p value) and determining whether the p value was less than the q value. This analytic strategy was followed in all patients. A matrix of statistically significantly differentially expressed genes was prepared. A gene was enriched in our cohort if it was deregulated in at least three patients.

Pathway analysis was performed using the Canonical Pathway module of Molecular Signatures Database (MSigDB, RRID:SCR_016863; ref. 38), an expertly curated, backend database used for gene set enrichment analysis (RRID:SCR_003199; ref. 39). We preferred MSigDB given its ease of usage, up-to-date backend curation and wide acceptance in the community, including pan-cancer analyses from The Cancer Genome Atlas (TCGA; refs. 40–42). Hallmark analysis was performed by submitting the selected gene lists to the H Hallmark module (43) of the MSigDB database. Gene ontology analysis was performed using the C5 module to select for biological process (BP) and molecular function gene sets of MSigDB. The top 10 statistically significant pathways involving these genes, identified by filtering against an FDR-corrected p value (q value) of <0.05, were considered for further analysis.

Transcript-level quantification was performed using Salmon (version 0.8.1, RRID:SCR_017036; ref. 44) on the RNA-Seq data (FASTQ files) using the default settings. A transcriptome index was built using the UCSC Homo sapiens reference genome (build hg19) GTF file followed by transcript quantification.

Low-count genes were filtered before running the DESeq2 (RRID:SCR_015687) function, and class comparison was carried out between groups. R v2.7 (RRID:SCR_001905) for all packages, v3.3 for packages that required newer R visualization libraries, and Bioconductor v3.4 (RRID:SCR_006442), with default settings were used for this analysis. RNA-Seq data were subjected to batch correction using log transformation, r-log transformation and variance stabilization transformation as implemented by DeSeq2 v1.6.3 (45) with default options. The results were visualized as dot plots. The “csDendro” command in CummeRbund R Package (version 2.7.2, RRID:SCR_014568) that implements hierarchical clustering using the Jensen–Shannon divergence distance matrix was used to assess intra-center–sequencing quality and batch effects and plot a genome wide correlation matrix. Differential gene expression analysis was performed using the DEseq2 based on a negative binomial generalized linear model. Transcript-abundance files from Salmon were imported and converted to gene-level information using the tximport (v1.0.3) R Bioconductor package (46), followed by class comparisons between samples obtained at time points, A, B, and C.

Log-fold change and p value (Wald test) were calculated for each gene for each comparison. Genes were considered differentially expressed if the Wald test p value was less than 0.05, and the log fold-change was greater than 2 in either direction.

Patients with treatment-naïve, histopathological-proven squamous carcinoma of head and neck (oral cancer) undergoing upfront curative surgery were eligible for this study. Patients were sampled at four time points during surgical intervention: T1 (fully vascularized tumor, corresponding to time point A in the breast cancer study), T2 (approximately midway during surgery, partially devascularized tumor, corresponding to time point B in the breast cancer study), T3 (immediately before final extirpation of the tumor from the body, nearly totally devascularized tumor, no corresponding time point in breast cancer study), and T4 (after tumor was completely removed from the body, fully devascularized tumor, corresponding to time point C in breast cancer study). Tumor-adjacent normal tissues were collected at T1 from all patients. RNA was extracted as described above.

Samples with low RIN (<7), indicating partial degradation of RNA, were processed using the Illumina WGDASL assay according to the manufacturer's recommendations. RNA samples with RIN ≥ 7 were labeled using the Illumina Total Prep RNA Amplification kit (Ambion) and processed according to the manufacturer's recommendations. Gene expression data were generated using Illumina BeadArray technology on HiScan.

Data were analyzed using GenomeStudio (v2011.1 Gene Expression module 1.9.0, RRID:SCR_010973). All assay controls were verified to ensure assay quality and chip-scanning integrity. Raw signal intensities were exported from GenomeStudio for pre-processing and further analyzed. The encrypted raw data (.idat files) were converted into a human-readable format using the Illumina gs plugin. The probe-level data were converted into gene-centric data and used for further processing. The Bioconductor lumi package was used for data preprocessing, which included quality-control steps, background correction, normalization, and log transformation, after which the data were ready for differential gene expression analysis. Robust spline normalization and variance stabilization transformation methods were used for normalization and transformation, respectively. To limit the analysis to highly variable genes, we removed unexpressed genes from the analysis. Normalized and transformed data for highly variable genes were further analyzed using BRB-ArrayTools (RRID:SCR_010938). Class comparisons between groups of arrays were based on univariate parametric testing, and the inclusion of genes was based on P values less than 0.05. Median absolute deviation analysis was performed to make the samples comparable with each other, and heat maps were generated using the Tigr Multi Experiment Viewer (MeV) module (version 4.9.0) of the TM4 Microarray Software Suite (RRID:SCR_001915; ref. 47).

We analyzed gene expression datasets from TCGA (RRID:SCR_003193) to test the hypothesis that surgical intervention induces distinct gene expression changes in excised human tissues. We did not identify any dataset with a mid-surgery tumor sample. Therefore, we evaluated the A versus C and T1 versus T4 comparison of our study by comparing the gene expression profiles of vascularized tissues (pre-surgery biopsies from tumors in situ) and de-vascularized tissues (tumor from surgically resected specimens). We identified seven cancer types in TCGA gene expression cohort (breast, cervix, diffuse large B-cell lymphoma, melanoma, mesothelioma, esophagus, and sarcoma) for which this information (tissue origin from pre-surgery biopsy or surgically resected specimen) was available. Patients in whom the biomaterial was procured using any of the labels as “incisional biopsy,” “core biopsy,” “segmental biopsy” were considered analogous to samples at time point A or time point T1 of our study. Patients where the biomaterial was procured using any of the labels as “lumpectomy,” “modified radical mastectomy,” “bilateral mastectomy,” “mastectomy,” “modified radical mastectomy and axillary dissection,” “segmental mastectomy,” “simple mastectomy,” “total mastectomy,” or “tumor resection” were considered analogous to samples at time point C or time point T4 of our study. The mRNA data and clinical information acquisition for all samples were accessed from TCGA portal on July 1, 2018. TCGA level 3 transcriptome data for the above were acquired and divided into two classes—biopsy and surgery. We performed class comparisons between the gene expression profiles of these two classes for each cancer type and identified significantly differentially expressed genes based on a parametric p value of <0.05. We then compiled a matrix of genes deregulated across all cancer types to identify a pan-cancer signature of the effects of surgical intervention.

We validated our RNA-Seq results in an independent set of breast cancer patient samples (n = 8) at time points A, B, and C, using a customized real-time qPCR gene array. Twenty genes, which were significantly deregulated in at least 50% of breast tumors in RNA-Seq analysis, were chosen for validation. An additional 23 genes known to be involved in canonical hypoxia, inflammation, and EMT (epithelial-to-mesenchymal transition) pathways were also selected for validation. A custom assay was designed for these 43 genes along with two reference genes, one genomic DNA control, one reverse transcription control, and one positive control for PCR. RNA was isolated using the TRizol reagent (Invitrogen) as per the manufacturer's instructions. The assay was conducted as per the manufacturer's instructions and is described in Supplementary Materials and Methods. Relative gene expression was analyzed using the ΔΔC_t method. A gene was considered validated if it was significantly deregulated in ≥4 of 8 patients with fold-change of ≥2 after applying one-way ANOVAV followed by the Tukey multiple-comparison test.

We designed a custom assay for 780 genes involved in the canonical biological pathways of metastasis, cellular invasion and migration, cancer, hypoxia, inflammation, immunity, drug response, drug resistance, stemness, cell cycle, DNA damage repair, cellular stress, wound healing, and the circadian cycle, which we tested in an independent set of 30 patients with breast cancer (time points A, B, and C) with an additional 13 patients from the RNA-Seq cohort (total n = 43). Of these genes, 156 were significantly deregulated in our RNA-Seq analysis, whereas the remaining 681 genes were not significantly deregulated in the RNA-Seq analysis. 16 housekeeping genes were included in the assay for reference-based normalization and batch correction across the samples. RNA QC was carried out using an Agilent 2100 Bioanalyzer (Agilent RNA 6000 Nano Kit) to measure the total RNA concentration, RIN, DV200, and fragment size. Samples with mRNA DV200 <0.3 were excluded from further assays. Samples passing QC were run on the nCounter Analysis System (NanoString Technologies, RRID:SCR_021712), according to the manufacturer's protocol as described in Supplementary Materials and Methods. NanoString nsolver software (version 4.0, RRID:SCR_003420) was used to analyze the data. The samples were run in multiple batches, which required batch normalization and correction. We performed technical normalization using positive control spikes to correct for any experimental variability between samples. Additional background correction was performed by subtracting the geometric mean of eight negative control probes. The expression values were normalized to the geometric mean of the 16 housekeeping genes. The expression values were log₂ transformed and standardized within each sample. Replicate analyses were performed by grouping samples at time points A, B, and C and evaluating gene expression changes A versus B, B versus C, and A versus C. A p value was computed for pairwise distribution using the unpaired t test, and the fold change for the difference was calculated. Transcripts with statistically significant changes (P < 0.05) were considered deregulated.

The human sequence data generated in this study are not publicly available due to patient privacy requirements but are available upon reasonable request from the corresponding author. Other data generated in this study are available within the article and its Supplementary Data Files.

This study was conducted in accordance with the Declaration of Helsinki. Patients were recruited after obtaining written informed consent under the appropriate protocol for each study as approved by the Tata Memorial Hospital (TMH) and ACTREC Institutional Ethics Committee (IEC) for IEC studies no. 168, 254 (ClinicalTrials.gov identifier NCT03797482), 1259 and TMH National Tumor Tissue Repository Project (NTTR).

Where required, all software versions and relevant settings are mentioned in the article. If required, all computational pipelines used for analysis in this article, will be made available as versioned bash files that can be executed on UNIX platforms upon a request to the corresponding author.

90 patients with breast cancer were enrolled on TMC-IEC approved studies after written informed consent. Patient sampling is shown in Fig. 1A Nine of 90 patients did not have at least one analyzable tissue sample and were screen-failed and excluded from analysis (Fig. 1B). Five out of 81 patients were biopsied 10–15 days apart without a mid-surgery sample, including a diagnostic biopsy at disease presentation (representing time point A) and a resected tumor sample at surgery 10–15 days after the diagnostic biopsy (time point C). Samples from 46 patients were subjected to whole-transcriptome paired-end RNA-Seq on an Illumina platform, those from 43 patients (including some that also underwent RNA-Seq) were subjected to a custom NanoString nCounter gene expression profiling assay, and those from 8 patients were subjected to a custom RT² profiler qPCR assay. The RNA-Seq and NanoString nCounter assays had an overlap of 13 patients, although samples did not overlap because only those that did not pass RNA-seq QC were subjected to nCounter assays provided they passed QC for the latter. The NanoString nCounter and RT² profiler assays had an overlap of 3 patients with an overlap of all samples from these patients. There was no overlap between the patients subjected to RNA-Seq and the custom RT² profiler qPCR assay.

The clinical characteristics of the included patients with breast cancer and availability of various time point samples for the three assays from each patient are shown in Supplementary Table S1.

The 34 patients with analyzable samples in at least two time points were eligible for this analysis (Supplementary Table S2; Fig. 1A and B). However, samples at time point C from 2 patients were mislabeled and/or interchanged and excluded from this analysis because this was a patient-specific analysis. These patients were included in the replicate analysis (Table 1). A time-series RNA-Seq analysis resulted in 67 comparisons: A versus B, 23 comparisons; B versus C, 18 comparisons; and A versus C, 26 comparisons. Four of these 67 comparisons were for tumor-adjacent normal tissue from two patients (patient 51N: three time points A, B, and C, patient 95N: two time points A and B).

We identified 4,947 significantly deregulated genes (P < q value; FDR = 0.05) across 67 comparisons in 32 patients (Fig. 1C). A cutoff value of statistically significant deregulation in at least three patients for each comparison (A vs. B, B vs. C, A vs. C) identified 737 deregulated genes (523 in A vs. B, 423 in B vs. C, and 250 in A vs. C, with overlaps; Fig. 1D). A total of 141 genes were significantly deregulated in all three comparisons (Fig. 1E and F). Notably, 487 of these 737 genes (235 unique to A vs. B, 148 unique to B vs. C, and 104 common to A vs. B and B vs. C) were transiently deregulated during surgery and would have been missed if samples at time point B were excluded from this study.

Genes related to stress responses (FOSB, FOS, JUN, JUNB, DUSP1, CYR61, EGR1, EGR2, EGR3, ATF3, RGS1, RGS2, and ZFP36), inflammation (GABRP, PIGR, NR4A1, NR4A2, CCL19, CCL21, and SLPI), epithelial markers (KRT14, KRT15, KRT17, KRT5, and KRT7), cellular invasion and migration (MMP9, MMP11, and MMP13), and lipid and fat metabolism (CD36, LIPE, CYP4B1, CYP4Z1, FABP4, PLIN1, PLIN4, LPL, and LEP) were among the 523 significantly deregulated genes in A versus B comparison, after partial blood vessel ligation when the tumor surface was semi-vascularized leading to acute hypoxia. Most of these genes were among the 141 genes found to be significantly deregulated in all three comparisons (Fig. 1E). Hallmark module analysis using the molecular signature database (MSigDB; see Supplementary Materials and Methods) revealed the enrichment of TNFα signaling via NF-κB, EMT, and hypoxia as the most enriched hallmarks (Fig. 1G; Supplementary Table S3A).

Pathway analysis of 523 deregulated genes between time points A and B revealed enrichment of genes belonging to the activator protein 1 (AP-1) transcription factor-network and extracellular matrix (ECM) protein-related genes (Fig. 1H; Supplementary Table S3B). The AP-1 transcription factor network intricately interacts with inflammation and hypoxia and mediates cellular hypoxic effects (48–52). Further analysis revealed that most deregulated genes were involved in Gene Ontology pathways such as the regulation of cell death, tissue development, circulatory system development, and response to oxygen (Fig. 1I; Supplementary Table S3C).

We identified 423 significantly dysregulated genes between time points B and C, of which 318 genes were also significantly deregulated in A versus B or A versus C comparisons, whereas 148 genes were uniquely deregulated in only B versus C comparison (Fig. 1E). Hallmark analysis of 423 significantly deregulated genes revealed sustained changes associated with inflammation, hypoxia, angiogenesis, and stress response (Fig. 1G; Supplementary Table S3A). Pathway analysis revealed that six of the top 10 deregulated canonical pathways were related to ECM degradation (Fig. 1H; Supplementary Table S3B). Gene ontology analysis identified genes involved in responses to external stimuli, receptor binding, and hormone responses (Fig. 1I; Supplementary Table S3C).

We also evaluated the gene expression changes that occur between a fully vascularized (time point A) and a completely devascularized (time point C) tumor, and identified 250 genes that were significantly deregulated (Fig. 1E), of which 36 genes were uniquely deregulated only in this comparison. Hallmark and Pathway analyses of the 250 deregulated genes largely recapitulated the changes observed in its two component comparisons, A versus B and B versus C (Fig. 1G and H; Supplementary Table S3A and S3B). Gene ontology analysis identified BPs associated with responses to oxidative stress and previously identified responses to external stimuli (Fig. 1I; Supplementary Table S3C).

We included RNA-Seq data of treatment-naïve, diagnostic biopsies of 17 additional patients with breast cancer from a previous study at our center (27) in this analysis, representing time point A, to increase the statistical power. Thus, a total of 63 patients (46 from this study, and 17 from the previous study) were eligible for this analysis (Table 1). We subjected these data to three batch-correction methods (Supplementary Fig. S1). Post batch-correction samples subjected to unsupervised clustering on the basis of genome-wide batch-corrected FPKM counts showed a clear batch effect, with a sequencing center-wise clustering pattern visible in the dendrogram (Supplementary Fig. S2), indicating that batch correction was unable to adjust for these effects.

Class-comparison analysis of the three groups, A (n = 63), B (n = 28), and C (n = 32) identified 322 unique genes (410 total) that were significantly deregulated in A versus B (n = 246), B versus C (n = 64), and A versus C (n = 100; Supplementary Table S4). Three genes were significantly deregulated in all 3 comparisons, 22 genes in A versus B and B versus C comparisons, 55 in A versus B and A versus C comparisons, 5 in B versus C and A versus C comparisons (Supplementary Fig. S3), 166 uniquely in A versus B, 34 uniquely in B versus C, and 37 uniquely in A versus C, comparisons. Overlap analysis of 322 deregulated genes from replicative analysis with 737 deregulated genes in individual patient-specific analysis identified 85 genes that were commonly significantly deregulated in both analyses, with Benjamini–Hochberg corrected P value 3.469149e-22, computed with DynaVenn (53). Gene ontology analysis of these 85 genes using ClueGo (54) on CytoScape (55) revealed their involvement primarily in the BPs of the p38MAPK cascade, cellular proliferation, cellular differentiation, and response to hormones (Supplementary Fig. S4; Supplementary Table S5).

Ten patients with head and neck cancer who were eligible were recruited to this study after written informed consent. Nine of these 10 patients had tissue with mRNA of sufficient quality available at all four time points (N1, T1, T2, T3, and T4) and were included in this analysis (Table 1). The mean age was 42 years (range, 22–62 years), 8 (88.9%) were of biologically male sex (n = 8), 8 (88.9%) had history of tobacco use, and 2 (22.2%) had history of alcohol use, and all 9 (100%) had oral cavity squamous carcinoma. One patient had tissue with inadequate quality at all time points and was excluded from the analysis. Detailed clinical information is provided in Supplementary Table S6.

We carried out a replicate analysis of HNSC samples, like the one in breast cancer cohort. Tumor samples at four time points (T1–T4) from patients P2–P10 were considered replicates and each time point was considered a class. All samples were normalized with respect to each other's expression profile across the four classes. We performed multiple-class comparisons among classes T1–T4. As a QC, we clustered all samples based on their genome-wide expression profiles. All samples from the same patient clustered together with some exceptions (Fig. 2A). On the basis of a P value cutoff (P < 0.05), we identified 16 genes (FOSB, EGR1, CTGF, DUSP1, SNORD3A, FOS, SNORD3D, CYR61, KLF9, JUN, EGR2, DPYSL3, RFC4, RNY5, GSDMA, and CKS2) as differentially expressed in at least one possible comparison between T1, T2, T3, and T4 (Fig. 2B). Next, we conducted binary comparisons between T1 for all patients (treated as replicates) and T2, T3, and T4, respectively, for all patients (treated as three separate replicates). In this comparison, each time point was normalized only for samples within that time point. This analysis produced binary comparisons between two classes regardless of the effect of other time points on that comparison and resulted in 30 significantly deregulated genes in T1 versus T2, 20 in T1 versus T3, and 50 in T1 versus T4 comparisons, respectively (Fig. 2C–E; Supplementary Table S7). Six genes (KLF9, EGR1, FOS, DUSP1, CTGF, and FOSB) were upregulated in all three comparisons (Fig. 2F).

A class comparison of mRNA profiles of diagnostic biopsies with surgically resected tissues in seven cancer types from TCGA study (Fig. 3A) identified AP-1 transcription factor-network genes (Fig. 3B; FOSB, FOS, EGR1, RGS1, DUSP1, NR4A1, CYR61, CTGF, EGR2, and EGR3) as the top differentially expressed genes, in five of seven cancer types evaluated—seen in breast, esophageal, cervical cancers, mesothelioma, and diffuse large B-cell lymphoma but not in melanoma and soft tissue sarcoma. This observation in multiple cancer types in an independent dataset increased our confidence that AP-1 transcription factor pathway genes are frequently deregulated in human tissues after surgical intervention.

We validated these changes using qPCR in an independent cohort of eight patients (Table 1). Ten of the 17 genes (9/15 upregulated and 1/5 downregulated) from the discovery cohort were validated as significantly differentially expressed (Fig. 3C; Supplementary Table S8) in one or more of the three comparisons (A vs. B, B vs. C, A vs. C). These genes included AP-1 transcription factor-network family genes and those involved in inflammation. Of the 23 candidate genes chosen from literature based on their involvement in hypoxia, inflammation and EMT, three genes, including the hypoxia-related genes HIF1-A and VEGF-A, were significantly differentially expressed (upregulated) in all three comparisons (A vs. B, B vs. C, and A vs. C).

On the basis of these results, we sought to further validate the findings using a custom NanoString nCounter assay of 780 genes. The NanoString validation cohort consisted of samples from 43 patients, of whom, 13 patients were also a part of the discovery cohort (RNA-Seq). Samples from these 13 patients could not be subjected to RNA-Seq due to QC failure, but successfully passed QC for NanoString assays. These samples included 32 at time point A, 36 at time point B and 34 at time point C (Table 1). Of note, no samples were subjected to both RNA-Seq and NanoString assay.

Using this assay, we observed that 82 unique genes (24 from our RNA-Seq data and 58 candidates from literature) were significantly deregulated (P < 0.5; Fig. 3D) in A versus B (n = 47), B versus C (n = 12), and A versus C (n = 62) comparisons, respectively (Supplementary Table S9). The overwhelming majority of the 24 validated genes from our RNA-Seq data were involved in the AP-1 transcription factor-network pathway, with few genes involved in inflammatory response pathways. The literature-based 56 genes were involved in regulating hypoxia (HIF1A and IDH1), ion channels, cell-cycle progression, cellular differentiation, stemness, the NF-κB-TNF-A axis, angiogenesis, fibrinolysis, the angiotensin pathway, cellular invasion and migration, the circadian cycle, DNA damage response (DDR), cell signaling, and epigenetic modulation. Surprisingly, we identified two immune checkpoint genes (CTLA4 and PDL1) that were upregulated during surgical resection.

We performed QC steps, background correction, normalization, log transformation, and class comparison for all nine head and neck cancer patient tissues at time points T1, T2, T3, and T4, along with the tumor-adjacent normal tissue biopsied at time point N1. We compared N1 versus T1, N1 versus T2, N1 versus T3, and N1 versus T4 and identified AP-1–related genes, EGR1, FOS, DUSP1, FOSB, JUN, CTGF, CYR61 to be upregulated in the tumor tissue in all comparisons, except T1 versus N1 (Supplementary Fig. S5; Supplementary Table S10). This suggests that AP-1–related genes were at the basal expression level in normal tissue and tumor at time point T1, with surgical stress subsequently upregulating these genes at various points during surgery. These findings are supported by NanoString assay of 20 tumor-adjacent normal breast tissue samples at time point A, 17 normal breast samples at time point B, and 15 normal breast samples at time point C that were compared with each other (normal breast versus normal breast at time points A, B, and C; Table 1). We identified 91 unique deregulated genes in tumor adjacent normal tissue as surgery progressed in A versus B (n = 50), B versus C (n = 9), and A versus C (n = 75) comparisons, respectively (Fig. 3D; Supplementary Table S10). Of these 91 genes, 40 were also significantly deregulated in breast tumor samples at time points B and C (DynaVenn computed P value 2.575964e-09) in the NanoString assay. These 91 genes included all AP-1 transcription factor-network genes and those involved in Inflammation, EMT, and cellular differentiation.

We evaluated the role of deregulation of hypoxia-related genes during surgical procedure by overlapping the combined list of intraoperatively deregulated genes (n = 1,113 unique genes) obtained from RNA-Seq of breast tumors (n = 975 genes; 737 from individual patient analysis, Supplementary Table S2 and 322 from replicate analysis, Supplementary Table S4), microarray of HNSC (n = 80 genes, Supplementary Table S7), and NanoString assay of breast tumors and normal tissues (n = 135 genes, Supplementary Table S9) with 2,009 proteins known to be deregulated in hypoxia (56). This led to the identification of 164 common genes (Supplementary Table S11) that are intraoperatively deregulated in our study and known to be involved in chronic hypoxia.

We studied gene expression changes in resected breast cancer and HNSC tissues from patients at different time points during surgical resection using genome-wide transcriptomic methods. A comparison of the transcriptional states of tumor tissues immediately before, during, and immediately after surgery revealed deregulation of genes involved in inflammatory responses (TNFα via NF-κB), EMT, and hypoxic stress response, which underlie many of the hallmarks of cancer. This is consistent with previous studies showing that hypoxic response can activate the AP-1–signaling cascade (57) and NF-κB pathway (58), inflammatory responses can activate TNFα signaling (59) and that established mediators of inflammation are deregulated in cancer (60, 61). Hallmark and pathway analysis using MSigDB identified pathways and genes involved in 7 of the 10 hallmarks of cancer, namely, Tumor-Promoting Inflammation (TNF-A, NFK-B, IL18 pathways), Activation of Invasion and Migration (various ECM related pathways, cell migration), Sustained Proliferative Signaling (K-Ras signaling), Evasion of Growth Suppressors (P53 signaling, regulation of cell death), Deregulating Cellular Energetics (response to lipid, secreted factors, adipogenesis), Inducing Angiogenesis (hypoxia signaling, myogenesis), and Avoiding Immune Destruction (CTLA-4 and PD-L1), to be significantly deregulated during surgical resection.

The custom NanoString nCounter assay recapitulated the gene expression changes observed in the RNA-Seq discovery cohort in an independent cohort of breast cancer and tumor-adjacent normal samples collected during surgery. The coherence of intraoperative gene expression changes in tumor and its adjacent normal tissue in two cancer types (breast cancer and HNSC) suggests that these changes could be a common response of tissue subjected to surgery. This is further supported by the finding that the top differentially expressed genes between biopsy and surgery specimens from five out of seven different cancer types in the TCGA dataset were nearly overlapping with those identified in our study.

The upregulation of AP-1 transcription factor-network genes was validated by qPCR in an independent cohort of patients with breast cancer. In addition, the custom NanoString nCounter assay identified several intraoperative gene expression changes central to biological pathways of EMT, stemness and cellular differentiation, inflammation, cell-cycle regulation, and, surprisingly, tumor immune checkpoints (PD-L1 and CTLA-4). The latter is an intriguing finding that needs to be further explored for several reasons, including the fact that these genes could potentially enable tumor cells to evade immune surveillance during the perioperative period. The activation of pro-survival gene expression programs in cancer cells during surgery could adequately prime them for metastatic spread. We identified 164 genes de-regulated during surgery, which are also known to be involved in hypoxia pathways as reported in an expert-curated database, HypoxiaDB. These 164 genes could plausibly constitute a backdrop in transcriptional studies on cancer. The majority of intraoperatively deregulated genes identified in our study have not been previously reported to be involved in hypoxic stress. It is possible that they are either deregulated due to surgical stress or are yet to be associated with acute onset hypoxia, or a combination of both. It would be appropriate to perform a regression analysis to identify a gene expression signature that is able to predict surgical sample class. We have not attempted this in the current analysis because of the relatively small sample size of the validation cohort. This is being planned in an ongoing study (ClinicalTrials.gov identifier NCT03797482) wherein gene expression analysis will be performed on intra-operative samples from 500 patients with breast cancer.

Other investigators have reported similar gene deregulation patterns from surgical stress. In their seminal Nature report, Sorlie, Perou and colleagues (62) mentioned that a cluster of genes, including FOS and JUN, could have deregulated expression due to prolonged sample handling after surgical resection. Sortiriou and colleagues (63) classified breast cancer tumors into basal 1 and basal 2 subtypes based on JUN, FOS, and FOSB expressions. JUN, JUNB, JUND, DUSP1, NR4A1, and CXCR4 were dysregulated in pre-surgical and post-surgical samples from 12 patients with prostate cancer who underwent proctectomy (64). DUSP1, FOSB, and NF-κB were also identified as the top deregulated genes between preoperative and postoperative breast cancer tissues obtained 2–8 weeks apart from 13 Norwegian patients with breast cancer (65). A study investigating the effect of surgery on gene expression in and tumor adjacent normal tissue from three patients with breast cancer undergoing surgical resection identified JUNB, EGR1, EGR3, ATF3, CYR61, and CXCL2 among top upregulated genes in tissues collected immediately post-induction and at end of surgery (66). A study that investigated the effect of cold and warm ischemia in colorectal tissue identified FOSB, EGR1, DUSP1, ATF3, IL8, RGS1, and SOCS3 to be differentially expressed before and after clamping of the hepatic artery (67). FOS, FOSB, and JUN were differentially expressed between core cuts obtained at the end of surgical resection and those obtained after a 90-minute delay from resected breast tumor specimens (68). These three genes were also differentially expressed between diagnostic biopsies and surgical resection specimens in an independent breast cancer dataset in the same study. Another study identified deregulation of JUN, FOS, EGR1, MMP9, PTGS2, and IL6, between percutaneously obtained and surgical specimen-derived normal liver tissue (69). FOS and FOSB are among the most deregulated genes, when investigating the effects of the timing of fine needle-aspiration (before and after surgery) biopsies on gene expression profiles in breast cancers (70). FOS, FOSB, JUN, EGR1, DUSP1, CYR61, NR4A1, NR4A2, NR4A3, PTGS2, and ZFP36 were differentially expressed between diagnostic breast cancer biopsies (from tumors in situ) and resected tumor specimens in 37 patients (71).

There is some evidence to suggest that these genes are upregulated by ischemic changes. FOS and ZFP36 are the most deregulated genes because of ischemic changes induced by cold blood cardioplegia in the human heart (72). Another possibility is that these changes are surgical-wound related. A study of early transcriptional changes in wound healing using a rat model, identified FOS, FOSB, EGR1, EGR2, and EGR3 among the top 15 genes deregulated within a span of 5 minutes after mechanical loading of the rat achilles tendon (73).

The roles of the deregulated AP-1 transcription factor-network genes have been well established in various biological pathways. AP-1 transcription factor is a dimeric complex comprising members of the JUN, FOS, activating transcription factor (ATF), and musculoaponeurotic fibrosarcoma protein families. AP-1 proteins dimerize through a leucine zipper motif and contain a basic domain for interaction with the DNA backbone. Thus, the AP-1 complex can form different heterodimers and homodimers that determine the genes regulated by AP-1 (74). The main AP-1 proteins in mammalian cells are FOS and JUN. We observed upregulation of four primary members of the AP-1 family (FOS, FOSB, JUN, and JUNB) in all our datasets—breast cancer RNA-Seq, HNSC microarray, breast cancer NanoString, breast cancer qPCR, and TCGA five cancers. The AP-1 family is implicated in diverse processes such as inflammation (75, 76), cancer (77, 78), invasion (79, 80), migration (81), angiogenesis (82, 83), chemotherapy resistance (84–86), and targeted-therapy resistance (87–89). AP-1 transcription factor-network genes have also been identified as central to the angiogenic switch in pancreatic cancer (90). Clinically, AP-1 pathway inhibition leads to improved outcomes following neonatal hypoxic-ischemic brain injury (91), and it is a therapeutic target for inhibiting COX2-mediated inflammatory pathways (92).

Our study has some strengths, most notably the unique sample set and analytical techniques. A within-patient analysis of sample obtained at different time points from each patient enabled identification of differentially expressed genes while avoiding inter-tumor heterogeneity (28–30). The secondary replicate analysis, clubbing all samples obtained at each time point together, allowed greater power because of higher sample numbers. These strategies were complimentary and the confidence in our findings is enhanced because of overlapping gene expression changes discovered by two different statistical models. Second, the effect of surgical stress on tissue gene expression is concordant between two cancer types included in our study and also concordant with biopsy-versus-surgical specimen comparison in five cancers from the TCGA study. However, without conducting our study design in all cancer types, it is not possible to definitively conclude that surgical stress will produce the same or similar gene expression changes across all cancers.

Our study also has some limitations. First, our findings represent observations from a relatively small cohort of patients with breast and head and neck cancers. Although these findings were replicated in an independent validation cohort using orthogonal technology and in TCGA cohort, they need to be systematically reproduced in a larger patient population. Second, the underlying mechanisms that are responsible for the molecular effect of surgery remain to be investigated and could be acute hypoxia due to blood vessel ligation during surgery, mechanical shearing stress, hypoglycemia, anesthesia, wound healing or a combination of these. In this context, it has been reported that anesthetic drugs can modify gene expression in tumors (93) and other tissues (94, 95). It is not possible from our study to discern the effect of surgery from that of anesthetic agents on tumor gene expression. Third, a more frequent short-interval tissue sampling during surgery could potentially provide a more complete picture of gene expression changes but would be difficult to conduct given the finite size of operable tumors and requirements of pathology grossing and clinical reporting of surgical outcomes. Fourth, inclusion of a control group of non-surgical samples, obtained at the same time points as our study but without surgery, would definitively attribute gene expression changes to the surgical procedure. However, carrying out such repeated biopsies in such a short period would be very difficult in human patients. Fifth, multiple testing in a limited sample size increases the possibility of false-positive results. However, the fact that we obtained similar gene expression results in multiple experiments on various technology platforms suggests that the findings of our study represent an actual signal rather than noise. Sixth, physiological and pathological implications of AP-1 deregulation need to be better delineated in the context of development of metastasis.

Our results show that gene expression findings based on either biopsy or surgically resected samples alone provide an incomplete understanding of tumor–treatment interactions, which could potentially impact the natural history of tumors. The magnitude of expression profile changes underscores the sensitivity of tumor tissue to its microenvironment. To our knowledge, this is the first systematic study to characterize the effects of surgical stress on gene expression profiles of human breast and head and neck tissues (cancer and normal) using sequentially acquired samples, including an intraoperative one, and our findings could have profound implications for understanding of solid tumor biology.

No disclosures were reported.

The funders had no role in the design and conduct of this study or interpretation of the results.

R. Chaubal: Data curation, software, formal analysis, supervision, validation, investigation, visualization, writing–original draft, project administration, writing–review and editing. N. Gardi: Software, formal analysis, methodology. S. Joshi: Data curation, investigation, methodology, project administration. G. Pantvaidya: Resources, formal analysis, supervision, investigation, methodology, project administration. R. Kadam: Data curation, investigation, methodology, project administration. V. Vanmali: Resources, supervision, investigation, methodology, project administration. R. Hawaldar: Data curation, supervision, methodology, project administration. E. Talker: Data curation, investigation, methodology. J. Chitra: Data curation, investigation, methodology. P. Gera: Validation, investigation, methodology. D. Bhatia: Validation, investigation, methodology. P. Kalkar: Validation, investigation, methodology. M. Gurav: Validation, investigation, methodology. O. Shetty: Validation, investigation, methodology. S. Desai: Validation, investigation, methodology. N.M. Krishnan: Validation, investigation, methodology. N. Nair: Resources, supervision, investigation, methodology, project administration. V. Parmar: Resources, supervision, investigation, methodology, project administration. A. Dutt: Methodology. B. Panda: Supervision, validation, investigation, methodology. S. Gupta: Resources, data curation, formal analysis, supervision, validation, investigation, methodology, writing–original draft, project administration, writing–review and editing. R. Badwe: Conceptualization, resources, data curation, supervision, funding acquisition, validation, investigation, methodology, writing–original draft, project administration, writing–review and editing.

This study was funded by the Department of Atomic Energy, Government of India (GOI). We acknowledge the generous support and funding received from Sunil Gupta for the NanoString nCounter Assays. We also acknowledge funding from Mizuho Bank Limited for research infrastructure at our institute. We are thankful to Akhil Gupta for funding laboratory infrastructure. We acknowledge part research funding for this study from the Women's Cancer Initiative (WCI)—Tata Memorial Hospital. R. Chaubal and N. Gardi were funded by a fellowship from HBNI, Mumbai, and TMC, Mumbai. R. Chaubal, N. Gardi, and R. Kadam are also funded by a fellowship from the Department of Biotechnology (DBT), GOI, through the DBT-Virtual National Cancer Institute (VNCI) Breast Cancer 2015 grant (BT/MED/30/VNCI-Hr-BRCA/2015) awarded (to S. Gupta). N. Gardi is also funded by the Department of Science and Technology (DST)—Scientific Engineering and Research Board (SERB) Prime Minister's Fellowship.