This study was approved by the Central Ethics and Scientific Committee at the National Cancer Institute in Mexico City (INCan) (015/012/ICI, CEI/961/15) and has been conducted in agreement with the ethical standards as laid down in the 1964 Declaration of Helsinki and its later amendments. A total of 92 biopsies from patients with LACC cancer were obtained. Tumor samples were collected from 2014 to 2018 by the Pathology Department, INCan, Mexico City. After confirmed diagnosis, all patients received concurrent chemoradiotherapy using cisplatin [weekly cis-diamminedichloroplatinum (II) at a dose of 40 mg/m²] for a total of five or six cycles and radiation (external radiation and intracavitary brachytherapy, for a total dose of 64–66 Gy over 67 days) (17). The patients’ therapy response was assessed according to RECIST criteria defined as follows: the disappearance of all target lesions was assigned as complete response (CR); meanwhile, patients with partial response, progressive disease, or stable disease were considered as therapy resistant (TR). The biopsies were divided into two cohorts: the first with 22 patients (12 CR and 10 TR) used as a discovery cohort to generate a microarray specific for CpG islands Array-Based Profiling of Reference-Independent Methylation Status (aPRIMES) (18); the second cohort, with 70 biopsies (40 CR and 30 TR), used for molecular data validation. The patient eligibility criteria consisted of (a) confirmed pathological diagnosis of CC stages from II-B to IV-B (LACC), (b) biopsies with a pathology report confirming more than 80% tumorous cells, (c) age range of 29–65 years, (d) high-quality DNA and RNA samples, (e) no other comorbidity, (f) no previous oncological treatment, and (g) patients able to receive the standard therapy based on concurrent chemotherapy and radiotherapy.

The DNA extraction from the 92 biopsies was performed as follows: 20 mg of fresh tissue was placed in a Fisherbrand Bead Mill homogenizer, and 2 ml soft tissue homogenizing Mix Tube was preloaded with lysis buffer [10 mM Tris–HCl, 2 mM ethylenediaminetetraacetic acid (EDTA), 1% sodium dodecyl sulfate (SDS)]. The tissue was homogenized using the MagNA Lyser instrument at 6,000 rpm for 1 min. To purify the genomic DNA, the QIAamp DNA Blood Kit (Qiagen, CA, USA) was used according to the manufacturer’s protocol. Finally, the purified DNA was stored at −20°C.

We employed the 3x720K CpG Island Plus RefSeq Promoter Arrays (Roche, Penzberg, Germany). These arrays cover the annotated CpG islands and the promoters of the RefSeq genes derived from the UCSC RefFlat files (Hg 38). Then, the hybridization probes were synthesized by aPRIMES assay. Briefly, genomic DNAs were digest by MseI, and the fragments obtained were subjected to linker-mediated PCR as described by Klein and coworkers (19); later, through enzymatic digestion by methylated-sensitive and methylated-specific enzymes, we obtained a methylated and unmethylated fraction of DNA, which were labeled with Cy5 and Cy3 fluorophores, respectively, and competitive hybridizing in a Human DNA Methylation 3x720K CpG Island Plus (Roche) as shown in Supplementary Figure 1. Then, arrays were scanned in an MS200 scanner (Roche). Finally, the alignment of the images and the extraction of data were carried out using the software DEVA Project Manager—1.2.1 (Roche). Next, for each region in the array, we obtained a continuous numerical ratio that represents if a region is hyper- or hypomethylated; we termed this ratio as bi-weight (BW) and is represented with the following formula:

Finally, to determine the significative methylated regions between responsive and resistant tumors, we calculated the Student’s t-test for each region between both groups and considered as statistical significance those methylated regions with a p < 0.01. Then, selected regions were ranked in ascending and descending orders accordingly to the difference of the means for both groups. All statistical analysis were executed in R environment.

Differentially methylated regions were analyzed by Pathway enrichment analysis by using Webgestalt (20) and ReactomePA (21); we only considered pathways with a p < 0.05 as subject of regulation.

We obtained from the Genome Browser database (22) a sequence of 2,000 bp (1,000 bp downstream; transcription start site, 1,000 bp upstream) that included the promoter region from each analyzed gene. These sequences were analyzed using MethPrimer web tool from Urogene to determine the CpG density (23).

To determine the methylation status of selected genes (BRD9, CTU1, and DOCK8), genomic DNA from each sample was modified using Methylation-Direct EZ DNA Kit (ZYMO, CA, USA). DNA bisulfite treatment changed unmethylated cytosines to uracil, but the methylated bases remained as cytosines. Then, two PCR reactions were performed per sample using specific primers to determine the methylated (M) or unmethylated (U) DNA status. The list of primers and its characteristics are shown in Supplementary Table S1. The product of each reaction was analyzed in agarose gels and resolved in Minigel OWLTM Easy CastTM B2 system (Thermo Scientific, MA, USA). Later, the gel was stained with ethidium bromide and photo-documented using a Gel Doc EZ Imager transilluminator (Bio Rad, CA, USA).

Chi-squared tests were employed to determine the differences in the distribution of the methylation status of the genes and the clinicopathological characteristics, considering p < 0.05 as statistically significant.

Kaplan–Meier plotter was calculated using the survival package in R, where the significance testing was assessed using the log-rank test. Significance was considered as p < 0.05.

This study was approved by the Central Ethics and Scientific Committee at the National Cancer Institute in Mexico City (approval number 015/01271B/CEI/961/15). The 92 patients who were enrolled accepted and signed the informed consent. All patients received treatment based on cisplatin and radiotherapy as mentioned in Material and Methods. The median age was 48 years. Patients were classified following the last version of International Federation of Gynecology and Obstetrics (FIGO) staging criteria as II (51.8%), III (37%), and IV (11.2%) stages. According to the FIGO’s guidelines, 52 patients (56.52%) showed complete response (CR) to therapy; meanwhile, 40 (43.48%) exhibited therapy resistance (TR). The HPV-genotype of all patients was determined by nested PCR (24). Table 1 shows the clinicopathological characteristics; a supplementary table that compiles all clinical data is available as Supplementary Table 2.

The determination of DNA methylation patterns has been proposed as a prognosis predictor in several types of cancer (25, 26). In this study, we employed aPRIMES arrays to obtain the genome-wide DNA methylation patterns on both groups, namely, responsive (CR) and therapy-resistant (TR) LACC patients. Next, to establish the differentially methylated regions (DMRs) on CR and TR groups, we compared the bi-weight ratio values from each analyzed region, using a Student’s t-test. The results showed a methylated DNA profile between the two groups composed of 16,538 DMRs that corresponded to 7,957 unique regions, where 2,833 of them were hypermethylated, 5,881 were hypomethylated, and 757 regions had both hyper and hypomethylated DMRs regions (p < 0.05) (Figure 1). As expected, a global hypomethylation pattern across the genome was observed, where the distribution of these DMRs varied to a considerable extent depending on the chromosome (Supplementary Figure S2). We noticed that chromosomes 1 and 19 had the higher number of gene promoters with DMRs, 859 and 753, respectively (Supplementary Figures S2, S3). Next, we observed that clustering of these DMRs using the Euclidean distance algorithm could distinguish the TR (black bar) from CR LACC patients (green bar) (Figure 1). Since the bi-weight ratios are continuous variables, we transformed them to a Z-score to visualize the DMRs in a heatmap, which clearly shows a DNA methylation profile that includes 3,533 DMRs hypermethylated and 13,005 DRMs hypomethylated (Figure 1).

Furthermore, we were interested in evaluating the impact of the methylation profile in biological pathways. The gene set enrichment analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database revealed that multiple key carcinogenic pathways such as the PI3K-AKT signaling pathway, nuclear factor (NF)-kappa B pathway, RNA polymerase, and pathways associated with breast cancer were dysregulated as a consequence of differentially methylation profile (Figure 2). As the PI3K-Akt pathway was the most enriched, we focused on analyzing it in more detail. As shown in Figure 2, multiple key genes in this pathway were hypomethylated, including the insulin receptor substrate-1 (IRS-1) and oncogene JAK3. In contrast, genes such as RELA that inhibit the tumor growing were hypermethylated.

To select methylated genes as potential biomarkers of response to chemoradiation, we further narrowed the methylation profile by considering only those DMRs that showed hyper- or hypomethylation status for further analysis. The results showed 4,463 DMRs with these methylation patterns that correspond to 1,439 unique genes. Then, median bi-weight values of methylation from each gene in TR and CR tumors were compared to calculate the median difference (MD). Promoter regions with an MD upper than 1.4 times, corresponding to the 13 genes enlisted in Table 2, were chosen for further analysis. The promoter sequence of the selected genes was analyzed as mentioned in Section CpG Island Density Determination, and the three genes with highest CpG density were selected for validation. A CpG island is defined as a DNA region highest than 500 bp that contains 50% or more of CG dinucleotides (27). The promoter region of BRD9, CTU1, and DOCK8 genes showed 88%, 63%, and 90% of CpG density, respectively (Figure 3). Additionally, the MD for these promoter regions was 1.61, 2.25, and 1.73 for BRD9, CTU1, and DOCK8, respectively (Figure 3, boxplots). Interestingly, H3K4me3 mark (chromatin compaction mark) was found near to these promoter regions (Figure 3).

To validate the methylation levels of BRD9, DOCK8, and CTU1 promoter regions as therapy response biomarkers, an MSP assay was performed. Bisulfite-treated DNA from 30 TR and 40 CR tumors samples were used to analyze methylated status. The results showed that the BRD9 promoter region was methylated in all CR tumor samples (40 CR tumor samples, 100% of cases), whereas it was hemimethylated in 25 TR tumor samples and unmethylated in 5 TR tumor samples (83% and 17% of cases, respectively) (Figure 4A). Instead, the promoter region of CTU1 gene was detected to be unmethylated in all CR tumor samples (40 CR tumor samples, 100% of cases), while in 27 TR tumors samples, it was hemimethylated; in 2 TR tumor samples, it was unmethylated and only in 1 TR tumor sample that it was methylated (90%, 6.6%, and 3.4% of cases, respectively) (Figure 4B). On the other hand, the DOCK8 promoter region was unmethylated in 29 TR tumors samples (97% of cases), whereas in 31 CR tumors samples, it was hemimethylated, and in 9 CR tumors samples, it was unmethylated (77.5% and 22.5% of cases, respectively) (Figure 4C).

Additionally, a chi-square analysis was performed to compare the methylation status of BRD9, CTU1, and DOCK8 genes with demographic characteristics of LACC patients (Table 3). The methylation of the BRD9 gene was associated with tumor stages II and tumor size <5 cm. In contrast, unmethylation of the CTU1 promoter region gene was associate with stages II and with tumor size <5 cm. The unmethylation status of the DOCK8 promoter region showed an association with stages III–IV; however, no significant relationship was found between the methylation status of this promoter with tumor size.

Finally, we determined if the methylation status of BRD9, CTU1, and DOCK8 genes could be an OS and the PFS biomarker of LACC patients. The results showed a better OS (p < 0.0041) and PFS (2.28 months in the hemimethylated group, p < 0.0001) in patients with methylation of BRD9 promoter (Figures 5A, D). In contrast, worse OS (p < 0.025) and PFS (3.12 months in the unmethylated group p < 0.0001) was observed in patients with the methylation of the DOCK8 promoter (Figures 5B, E). Moreover, patients with a unmethylated CTU1 promoter showed a better OS and PFS (1.76 months in the hemimethylated group p < 0.0001) (Figures 5C, F). These data highlight that the methylation status of BDR9, CTU1, and DOCK8 have the potential as biomarkers of OS and PFS in LACC patients.