Ethical approval for the present study was granted by the Stockholm Regional Ethical Review Board (Stockholm, Sweden; approval no. 2009/366-31). Patients eligible for inclusion were those >18 years, accepted and planned for surgery involving resections of the pancreas or liver, pre-operatively estimated to last >4 h. Patients requiring long-term GC treatment pre-operatively or with any known disease of the glucocortid system were excluded. The patients recruited for the present study were enrolled after informed consent was obtained, orally and in writing, in accordance with The Declaration of Helsinki. Whole blood (4 ml) was collected from 24 patients who were to have major abdominal (pancreatic or liver) surgery. Samples were drawn during a pre-operative visit at the Department of Surgery, Karolinska University Hospital Huddinge (Stockholm, Sweden) in 2014. The present study did not allow for a power calculation of the required number of participants, as the distribution of SNPs in the patient population was not known.

Complications were graded according to Clavien-Dindo classification (21), which is a well-established grading system in the surgical community. It is an ordinal scale, ranging from 1 to 5, with 1 denoting only a minor deviance, from the expected postoperative course, i.e., additional antiemetic treatment. Grade 2 denotes need for pharmacological therapy not considered routine or standard post-operative treatment. Grade 4 denotes organ failure in any organ system and grade 5 denotes mortality of the patient due to any cause. Grade 3 signifies complications, such as anastomic leakage or wound dehiscence, necessitating re-intervention and is divided into two sub-classes depending on whether this re-intervention is performed without (3a) or with (3b) general anaesthesia. Complications graded ≥3b are regarded as major.

Pre-operative samples of whole blood were collected in EDTA tubes (BD Vacutainer; BD Diagnostics) and the erythrocytes were immediately treated with cold (4–6°C) ammonium lysis buffer (153 mM NH₄Cl, 10 mM KHCO₃ and 0.01 mM EDTA adjusted to pH 7.4-7.5) until tranluscent. After lysis and washing in PBS, leukocytes were resuspended in 350 µl of RLT Buffer from the Qiagen RNEasy kit (Qiagen AB) and were subsequently frozen. DNA was extracted using the Qiagen RNEasy kit in a QIAcube (Qiagen AB) according to the manufacturer's protocol. Concentration and purity assessments were performed using a NanoDrop 1000 (Thermo Fisher Scientific, Inc.). The 260 nm/280 nm absorbance ratios were in the range of 1.75-1.85.

The CpG island is located at chromosome (chr)5:143402505-143405805 [Human Dec 2013 (GRCh38/hg38) Assembly], i.e., GRCh38/hg38, Human Genome Reference Consortium, human genome build 38 (human genome 38) (12). Targeted BiS sequencing was performed by Zymo Research Corp. Samples containing approximately 400–1400 µg DNA were diluted and BiS-converted using the EZ DNA Methylation-Lightning kit (Zymo Research Corp.) at 98°C for 8 min and 54°C for 60 min. BiS-converted DNA was then PCR-amplified using ZymoTaq PreMix (Zymo Research Corp.) for library construction, using a microfluidics system (Fluidigm Access Array; Fludigm Corporation), the specific cycling conditions and final concentrations are proprietary to Zymo Research Corp. The primer sequences used are presented in Table SI. BiS-converted, amplified DNA underwent next generation sequencing using the Illumina MiSeq platform (Illumina, Inc.) and corresponding kits by Zymo Research Corp. A paired-end 300 bp configuration was used and libraries were loaded at 8 pM concentrations. The raw data was aligned for CpG methylation calling. Sequence reads were identified using standard Illumina base-calling software. Sequence alignment was made if the coverage exceeded 10 reads. Furthermore, in addition to the CpG methylation analysis, the full sequences were analysed for SNPs to the extent allowed for following BiS conversion using Bioconductor software package R (version 3.7, http://bioconductor.org/news/bioc_3_7_release/) with the packages Rqc (version 1.10.2), Rsamtools (version 1.35.2), msa version 1.12.0 and GenomicRanges (1.32.0).

Multifactorial analysis was done in R (version 3.5.0) (22). Consensus sequences from each patient were extracted from the sequence data obtained via BiS-sequencing, which used a threshold of 0.5 as minimum probability threshold in the multiple sequence alignment package (23). All biological parameters were assessed using multiple factor analysis in Factominer (version 2.4) (24). This multifactorial approach was selected as it enabled the analysis of individual and combined SNP effects, methylation patterns, patient characteristics and clinical data. No patient outliers skewed the model. No individual factors, except for the Clavien-Dindo classification, were able to explain the significant portion of the variance observed. Therefore, all other factors were excluded from further analysis.

Parameters with minimal risk of non-biological bias were used to construct the statistical model, including liver or pancreas surgery, age, sex, and CRP levels. Factors omitted at this stage were those including some degree of subjective assessment, for example ASA-classification or pre-operative chemotherapy. Those variables were instead subsequently fitted into the model as supplementary variables. Patients with postoperative complications (defined as Clavien-Dindo ≥1; n=14) were included for partial least squares (PLS) analysis which was performed using the PLS package. PLS was chosen because it can accept ordinal variables, including Clavien-Dindo classification. Moreover, Clavien-Dindo classifies postoperative complications and therefore the non-differentiation between patients without complications presented a lower boundary of measurement. This truncation of the outcome variable, complications (graded according to the Clavien-Dindo classification), produced a region with zero variance that disqualified it from appropriate factor analysis. Therefore, it was only possible to analyse the severity of complications and patients without complications were excluded from further analysis.

The Clavien-Dindo classification system was used as the response variable and all investigated nucleotide positions with more than zero variance were used as the predictor matrix. ‘Leave one out’ cross-validation in Factominer (version 2.4) was used to avoid overfitting and the number of components with the smallest mean square error of prediction was retained (24). Two components explaining 71% of the variance were retained. Scores were used as an estimate of contribution of each SNP. P<0.001 was considered to indicate a statistically significant difference and was used as the cut-off value.

Data are presented as medians with ranges or interquartile ranges, as appropriate. Analysis of differential allele comparisons was performed using an independent-samples Mann-Whitney U Test (SPSS version 26; IBM Corp.). P<0.05 was considered to indicate a statistically significant difference.

Genomic positions determined to be significant were acquired from the statistical analysis and were investigated using the University of California, Santa Cruz (UCSC) Genome Browser Human Dec 2013 (GRCh38/hg38) Assembly (12). A VISTA plot was used to visualise the gene (25). Allele frequencies (MAFs) were examined in the 1000 Genomes Phase 3, ESP and gnomAD (v2.1) cohorts (26–28)

Patient characteristics and clinical data are presented in Table I. The median age distribution was 65.5 (range 42–81) years and sex distribution was 11/13 male/female. Overall, 14/24 patients had a postoperative complication. The specimen histopathology demonstrated different clinical conditions of the liver or pancreas, as listed in Table I.

Of the parameters available for factor analysis, sex and benign disease were statistically significant (P<0.05) in dimension one and two with 7.1 and 6.7% of the variance, respectively. These effects were not, however, assessed to skew the model and consequently not considered in subsequent analyses. Hence, a multifactorial PLS regression-based approach was employed to assess the association of methylation with postoperative complications, using the full BiS-converted sequencing data and the Clavien-Dindo classification system only.

When examining all 319 CpG islands in the BiS-converted sequence, a homogeneous demethylated pattern was demonstrated. The whole region was demethylated, except for the last CpG island (position, chr5:143405794). There was no association between CpG methylation at this site and postoperative complications (Fig. 2).

The combined effect of multiple contributing nucleotide sites was assessed, taking the clinically relevant heterogeneity into account. The results demonstrated that 4 positions were significantly associated with postoperative complications (Fig. 2). These four sites were previously reported to be SNPs. The properties of these SNPs are presented in Table II, along with the summarised data from these aforementioned studies [obtained using the UCSC genome browser (Fig. S1)]. For ease, the SNPs are numbered 1–4. Of the four SNPs identified in our analysis, two were previously identified common SNPs; NCBI (National Center for Biotechnology Information) dbSNP (data base for Single Nucleotide Polymophisms) accession number rs10482614 (chr5:143402837; SNP1) and NCBI dbSNP accession number rs3806855 (chr5:143404564; SNP4). GRCh38/hg38 demonstrated that the highest population minor allele frequencies (MAFs; observed in any population, including the 1000 Genomes Phase 3, ESP and gnomAD cohorts (26–28) were 0.26 and 0.25, respectively. SNP1 consists of a minor allele (A), which deletes a possible methylation site (CpG) (16). SNP4 has previously been described to affect promotor activity (16). The same nucleotide as in SNP4 can also be part of a three-nucleotide deletion, (NCBI dbSNP accession number rs796817133), highest population MAF <0.01), which was not encountered in our cohort (29). NCBI dbSNP accession number rs1039242888 (chr5:143404507; SNP2) has not been thoroughly investigated but has previously been reported to have been encountered (highest population MAF 0.02). This SNP, SNP2, was only present in one patient in the present study and should therefore be interpreted with caution. NCBI dbSNP accession number rs904759782 (chr 5:143404514; SNP3) was previously reported to have a highest population MAF <0.01 but was not further investigated and therefore its clinical significance is unknown. SNP3 was fairly common in the present study with over 40% of the patients being carriers of this minor allele. The minor alleles of SNPs 1–3 all result in depletions of CpG sites. Bioinformatics data on transcription factor (TF) binding sites at these four specific positions were obtained, revealing five different TFs: EGR1, SMARCA4, TFAP2C, STAT1, and RBL2. These are listed in Table III along with their mechanisms.

The PLS analysis was repeated containing each allele (major or minor) for each SNP, one at a time and the findings from the pooled PLS analysis, as presented above were confirmed. Both homozygotes and heterozygotes for the minor alleles were included in the minor allele group. The influence of major vs. minor allele carriership on postoperative complications was compared and of the four SNPs which emerged significant in our PLS (above), only SNP4 was significantly different between the major and minor alleles (p=0.016) when analysing one factor, i.e. major vs minor allele carriership for this one SNP, at a time (Fig. 3).