CC tumor biopsies were obtained from 30 untreated women with histologically proven cervical SCC (n = 25) or ADC (n = 5) referred to the Gynecological Service of the Instituto do Cancer do Estado de São Paulo (ICESP, São Paulo, Brazil) (Supplementary Table 1).

All patients signed an informed consent form before biopsies were obtained. Histological analysis was carried out using hematoxylin and eosin-stained sections at ICESP.

According to the manufacturer's instructions, DNA and RNA were isolated from each cervical tissue sample using the DNA/RNA All Prep kit (Qiagen, Germany).

Isolated DNA concentration and purity were evaluated using a NanoDrop 2000 spectrophotometer (Thermo Scientific, USA), and samples were diluted to a final concentration of 50 ng/μl. HPV detection and genotyping were performed using 100 ng of each DNA sample and the Inno-LiPA HPV Genotyping kit (FujiRebio, Belgium), according to manufacturer instructions.

RNA isolation was followed by DNase treatment. RNA integrity and quality were evaluated using the 2200 Tape Station System (Agilent Technologies Inc, USA), and RNA integrity number ≥ 8 was considered the quality threshold. In total, seven cervical SCCs and three ADCs with < 20% necrotic tissue were submitted to RNA sequencing (Supplementary Table 1).

All qualified CC samples were submitted to library preparation using the TruSeq Stranded messenger RNA (mRNA) LT Sample Prep kit (Illumina Inc., USA) and paired-end (2 × 100 bp) sequencing using the HiSeq SBS Kit v4 and the HiSeq Flow Cell v4 with HiSeq 2500 system (Illumina Inc., USA) at the Centro de Genômica Funcional Esalq—University of São Paulo. On average, 85 million paired-end reads were sequenced per tumor sample (Supplementary Table 2).

Transcriptome data of HPV-negative non-tumoral cervical samples (named non-CC) were downloaded from the Gene Expression Omnibus portal (accession PRJNA454568) (10).

Because gene-gene interaction analysis must be performed with a sample size bigger than the one used in this study, we used The Cancer Genome Atlas (TCGA) biolinks to extract data from the cervical tumors (TCGA-CESC). We downloaded 280 CC transcriptome data from the public TCGA repository under the tag TCGA-CESC, including 249 cervical SCC and 31 cervical ADC samples. We further compared gene expression profiles from our study samples with the TCGA dataset (Supplementary Table 3).

RNA sequence data output (FASTQ files) was trimmed with Trim Galore (v. 0.0.4.0) (11) and aligned onto the GRC37/hg19 version of the human genome and the HPV-16 genome (NC_001526.4, GenBank) using Hisat2 (v. 2.1.0) (12) with the following parameters: “--dta --fr -q --no-mixed --no-discordant”. The latest human gene annotation dataset (hsapiens_gene_ensembl, available at http://grch37.ensembl.org) was obtained using the biomaRt package (v. 2.34.2) (13) from R statistical environment. SAM files were converted, merged, sorted, and indexed using Samtools utilities (14). The aligned reads were visualized using the IGV software (v. 2.3.82) (15). Reads counts were accessed using HTseq (v. 0.5.4p3) (16) applying the following parameters: “-m intersection-noempty -i id_gene -s reverse.”

Differential gene expression analysis was performed to compare between the normalized and filtered read count among three groups: cervical SCC, ADC, and non-CC (10) using the Bioconductor package edgeR (v. 3.20.9) (17), setting the log₂ fold change module (|log₂FC|) ≥ 2 and false discovery rate as ≤0.01. To overcome the batch effect due to combining data from independent cohorts, we used edgeR package and BatchEffect function from Limma Package, available at the Bioconductor project. Non-expressed and weakly expressed genes were defined as having ≤1 count (read) per million in n of samples in the smallest group for each comparison and were excluded from the differential expression analysis.

The regulatory network for both cervical SSC and ADC was inferred by computing mutual information between annotated TFs and all potential target genes using the RTN package (18). In the network architecture inferred, each TF was assigned to a list of candidate targets, which could be linked to multiple TFs. As regulation can occur from both direct (TF–target) and indirect interactions (TF–TF–target), the Algorithm for the Reconstruction of Accurate Cellular Networks was used as an additional step to enrich the regulatory networks with direct TF–target interactions (19). The resulting regulatory networks were plotted with Cytoscape (v. 3.7.2) (20).

Enrichment of Gene Ontology biological processes of cervical SCC and ADC gene expression profiles were investigated using Enrichr (21). The GSEA software (version 4.0.3.4) (22) was used to identify Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways showing an overrepresentation of up- or downregulated genes in cervical ADC or SCC. Briefly, an enrichment score was calculated for each gene set (i.e., KEGG pathway) by ranking each gene according to their expression difference using Kolmogorov–Smirnov statistic, computing a cumulative sum of each gene ranked in each gene set and recording the maximum deviation from zero as the enrichment score.

To validate genes with different expression levels revealed by whole RNA sequencing, we performed quantitative polymerase chain reaction (qPCR) of six selected differentially expressed lncRNA using RNA obtained from 20 cervical SCC and 5 cervical ADC samples (Supplementary Tables 1, 4). We selected for validation by qPCR lncRNAs with the highest or lowest log₂FC values. Complementary DNA synthesis from 500 ng of total RNA isolated from CC samples was performed using the SuperScript VILO kit (Life Technologies, USA). qPCR was carried out using the SYBR Green PCR Master Mix on an ABI 7500 real-time PCR system (Applied Biosystems, USA). LncRNAs transcript levels were normalized to the mRNA levels of the ribosomal protein lateral stalk subunit P0 (RPLP0) constitutive gene. The lncRNA gene annotation from LNCipedia (23), Ensembl transcript identification, qPCR pair primers pair sequences, and expected amplicon length are available in Supplementary Table 4. The Mann–Whitney statistical test was applied to compare cervical SCC and ADC distribution, considering α = 0.01 as the significance level. Statistical analyses were performed using GraphPad Prism version 7 for Windows (24) and ggplot2 package (25) in the R statistical environment.

The following databases of miRNA targeting mRNA predictions have been investigated to detect potential mRNA transcripts regulated by the miR-28-3p or miR-96-5p: mirDIP with the minimum score set to medium (26); mirMap with the default parameters (27); TargetScan with a parameter to search for conserved sites for miRNA families conserved only among mammals (28); miRDB with the default parameters (29); and DIANA-microT with a threshold set to 0.5 (30). If miR-28-3p or miR-96-5p returned any positive result for the TF genes in the MODY pathway, a bidimensional matrix was populated (Supplementary Table 5).

We developed a classifier based on TCGA data, considering the ADC histological type (possible responses: yes or no, when it was an SCC case) as our outcome of interest. According to differential expression analysis, the set of predictors was composed of those genes identified as enriched in cervical ADC gene profile followed by TRN analysis. We applied logistic regression to adjust the classifier and leave-one-out cross-validation (CV) technique to evaluate its performance in data not used for its adjustment. Briefly, on each CV iteration, the observations are randomly divided into a training and a test set, then we adjust a logistic model in the training data and estimate the risk (probability) of being classified as ADC on the test data. The predictive performance was evaluated by calculating the area under the receiver operating characteristic (AUROC) curve with 95% confidence intervals. Sensitivity (S), specificity (E), positive predictive value (PPV), and a negative predictive value (NPV) were measured for different estimated probability (p) cutoff points: one that maximizes model's sensitivity and specificity, and other aiming to improve model's specificity and, accordingly, to reduce the false-positive results (i.e., an SCC case classified as ADC). The latter approach was based on a cutoff point that distinguishes the 10% of the observations with the highest predicted risk for ADC. The adjusted model was also applied to our dataset, and S, E, PPV, and NPV were estimated according to the cutoff points described earlier. Analyses were performed using R software. Model training was carried out using the caret (Classification and Regression Training) package and model evaluation through pROC and ROCR packages.

The upregulation of genes involved in epithelium maintenance among cervical SCC samples was further endorsed by the TRN analysis, which identified regulons whose lncRNA genes are highly expressed in SCC in contrast to ADC or non-CC samples. Epithelial differentiation and homeostasis are governed by several TFs (37–39). Among our findings, POU Class 3 Homeobox 1 (POU3F1) (40), Basonuclin 1 (BNC1) (41), Human Forkhead-box N1 gene (FOXN1) (42), hes family bHLH transcription factor 2 (HES2) (43), and zinc finger protein 185 with LIM domain (ZNF185) (44) were previously described to be involved in keratinocyte differentiation, or skin development, homeostasis, and wound healing.

One important regulator of the expression of several genes involved in keratinocyte proliferation and differentiation is p63 (45), which is encoded by the tumor protein p63 (TP63) gene. We found TP63 upregulated in cervical SCC as compared with ADC and non-CC (Supplementary Table 6, Figure 2B). Indeed, p63 was previously detected by immunohistochemistry in 97% of cervical SCC, whereas it was absent in ADC (46). In this study, TP63 integrates the SCC TRN as a gene potentially regulated by BNC1 and FOXN1. Expression of BNC1 and ZNF185 have been previously described to be positively regulated by p63 (47, 48). Higher expression levels of BNC1 and SRY-Box Transcription Factor 15 (SOX15) genes have also been previously detected in SCC samples from other anatomical sites compared with their related non-tumoral tissues (49–51). In contrast, ZNF185 displayed a discordant expression pattern in cervical SCC compared with our findings (48). However, because the expression of the protein encoded by this gene was shown to vary widely among skin layers (44), we suggest that our findings must be investigated in-depth in the future. Two of the six TFs within the cervical SCC TRN were also previously reported to be positively regulated by p63. We thus conclude that the cervical SCC TRN proposed in the present study encompasses a set of TFs with functions associated with epithelial differentiation and homeostasis, supporting its role in cervical SCC.

In addition to the putative role of different TFs upon cervical SCC TRN, several lncRNAs genes regulated by these proteins might also affect integrity maintenance and differentiation of the epithelium. Using qPCR, we validated the higher expression of TINCR ubiquitin domain-containing (TINCR), CALML3 antisense RNA 1 (CALML3-AS1), DSG1 antisense RNA 1 gene (DSG1-AS1), and IGF-like family member 2 antisense RNA 1 (IGFL2-AS1) lncRNA genes in cervical SCC in comparison with cervical ADC samples. Nevertheless, it needs to be highlighted that, to date, the functions of most of the cataloged human lncRNA are not clarified due to the lack of specific literature.

Increased TINCR lncRNA levels have been found within differentiated superficial epidermis layers and were shown to positively regulate the expression of the MAF:MAFB TF dimer (52). Corroborating our findings, higher expression of TINCR in cervical SCC was previously observed (53), although these data contrast with data obtained from TCGA (54). All other validated lncRNA genes in our study were classified as antisense transcripts (55). We observed higher expression levels of the IGFL2-AS1 in cervical SCC compared with cervical ADC and non-CC. Higher levels of this antisense gene have been described in renal cell carcinoma (56) and gastric cancer (57); in contrast, a lower expression of IGFL2-AS1 was observed in breast ADCs (58). The CALML3-AS1 lncRNA was shown to increase proliferation, migration, and invasion of CC cells and to decrease the apoptosis of these cells. Indeed, high CALML3-AS1 levels were previously observed in CC samples of TCGA (59). In our study, lower CALML3-AS1 levels were detected in cervical ADC and non-CC compared with cervical SCC, indicating this molecule may be further explored as a molecular target for specific diagnosis of this histological tumor type.

In conclusion, although the function of most of the lncRNAs remains unclear, a selected set of such transcripts were successfully validated as highly expressed in cervical SCC as compared with ADC. Among them are TINCR and CALML3-AS1, whose role in CC is under investigation. Taken together, our data reinforce the molecular differences between cervical SCC and ADC, whose significance should be explored in the prognostic and therapeutic management of these tumors.

ADC arises from the cervical glandular tissue and leads to abnormal epithelium structure and production of secreted molecules. In this study, we describe a set of 10 TF genes upregulated in ADC as compared with SCC or non-CC. All TFs within the ADC TRN have been previously associated with the genesis of ADCs, in particular pancreatic ductal adenocarcinoma (PDA), with the exception of the insulinoma-associated protein 1 (INSM1). Although no literature regarding INSM1 in PDA, it was described to have a key role in pancreatic development (60), similar to the GATA binding protein 4 (GATA4) (61).

GATA4 gene expression levels were shown to vary in ADCs arising from distinct tissues (62, 63). In PDA, GATA4 protein levels were detected with higher intensity in women than in men (64). In addition to GATA4 gene expression, the peroxisome proliferator-activated receptor gamma gene (PPARG) (65), the hepatocyte nuclear factor-4-alpha (HNF4A) (66), the hypoxia-inducible factor 3 subunit alpha (HIF3A) (67), the MYB proto-oncogene (MYB) (68), and the single-minded 2 (SIM2) (69) genes were also previously described to be highly expressed in PDA. To our knowledge, this is the first report to identify the surprising similarity of TFs overexpressed in cervical ADC and PDA. Such similarity encourages additional investigations. Interestingly, although highly expressed in PDA, high levels of SIM2 in cervical SCC have been associated with better overall survival compared with tumors expressing low levels of this gene (70). The Meis homeobox 2 (MEIS2) gene encodes for a TF that was identified as a key signal transduction pathway regulator in breast ADC (71). Recently, MEIS2 was identified as highly expressed in cervical SCC and showed prognostic value (54). Therefore, we believe that additional studies must be carried out to clarify the possible roles of SIM2 and MEIS2 in the distinct histological types of CC.

The pathway enrichment analysis in cervical ADC unveiled the MODY pathway. MODY TFs have also been extensively investigated in the context of pancreas and liver development. In addition to HNF4A expression in PDA, HNF4A and HNF4G genes were reported to be positively regulated by HNF1 homeobox A (HNF1A) protein in the pancreatic cell (72). Corroborating our findings, HNF1A was shown to have higher expression levels in CC compared with the normal cervix (73), despite the fact that the histological type was not indicated in their study. Nevertheless, we may speculate that HNF1A was identified as differentially expressed in their study due to an overrepresentation of cervical ADC among samples. Another MODY pathway member, the gene encoding the TF pancreatic and duodenal homeobox 1 (PDX1), had the literature reviewed in the context of pancreas embryogenesis (74), including PAD (75). In an extensive analysis regarding pancreatic cancer, PDX1, MNX1, HNF4G, HNF4A, HNF1B, HNF1A, FOXA2, FOXA3, and HES1 were proposed as the key TFs associated with the pancreatic progenitor subtype (76). Except for FOXA2 and HES1, all of these TFs were identified to be overexpressed in cervical ADC (log₂FC ≥ 2 and false discovery rate < 0.05; Supplementary Table 8), providing stronger support to the hypothesis that cervical ADC and PDA share molecular similarities. It is important to stress, however, that such relationships regarding cervical ADC regulation are predictions and must be further investigated and validated.

The construction of a TRN including all MODY TFs highly expressed in cervical ADC revealed PAX6, HNF4A, HNF4G, and FOXA2 genes with the potential to positively regulate other protein-coding or lncRNA genes (Figure 4B). At the moment, the role of these TFs within the context of cervical ADC is unknown. However, the MODY pathway was previously associated with other ADCs, such as EAC and PAD. Our TRN analysis showed that PAX6 and HNF4G TF genes might positively regulate the expression of the lncRNA LINC02381 gene (also known as LOC400043). Although there is limited literature regarding LINC02381's biological role, it was described to have a suppressive effect upon the tumorigenesis of human colorectal ADC (77). These data indicate that MODY TF genes may be associated with other molecular processes in addition to those previously described in onset diabetes. Additionally, we were able to validate that LINC02381 and LINC01833 had higher expression levels in cervical ADC than in cervical SCC by qPCR.

Although the role of lncRNAs in CC development is not completely elucidated, it is well-established that they can compete for binding to mRNAs with other regulatory molecules, such as miRNAs (78, 79). In the context of HPV-associated tumors, several studies indicate that E6/E7 viral oncoproteins target lncRNAs, subverting cellular processes seminal to carcinogenesis (80). Previous in silico miRNA to mRNA targeting predictions analyses, followed by experimental validation, identified that miR-96-5p and miR-28-3p could directly bind to LINC02381 (79).

The inspection of the entire MODY pathway (KEGG ID hsa04950) regarding the highly expressed TFs in cervical ADC allowed for the identification of the following set of MODY TF genes regulated in sequence: PAX6, PDX1, HNF4A, HNF1A, HNF4G, and FOXA3 (Supplementary Table 10, Supplementary Figure 6). Surprisingly, according to at least one of five prediction tools, this set of TF genes within the MODY pathway is a potential target of either miRNAs miR-96-5p or miR-28-3p (Figure 6A, Supplementary Table 5). Interestingly, miR-96-5p has been associated with CC tumorigenesis progression by regulating the PTEN pathway (81).

Integration of MODY TRN and MODY TFs regulated in sequence in cervical ADC revealed that HNF4G and PAX6 might lead to increased expression of LINC02381. We speculate that the enhanced expression of LINC02381 in cervical ADC could dysregulate the repression role of miR-96-5p and miR-28-3p, acting as a sponge of such miRNAs (79). Therefore, the repressed MODY pathway previously silenced in the cervix would become active in cervical ADC (Figure 6B).

Taken together, our TRN analysis revealed a set of MODY TFs upregulated in cervical ADC. Among these, two TFs could positively regulate the lncRNA LINC02381. According to the literature, LINC02381 can directly bind to two miRNAs acting as a sponge, reducing the availability of miR-96-5p and miR-28-3p. Computational predictions showed that miR-96-5p and miR-28-3p might negatively regulate this set of TFs. Therefore, we speculate LINC02381 may act by decreasing the levels of miR-96-5p and miR-28-3p and promoting the upregulation of the MODY pathway in cervical ADC. Such findings, when considered as predictors in a classifier to differentiate ADC and SCC cases, resulted in encouraging results. Nonetheless, it is important to note that external validation, based on a larger dataset, is required to corroborate our results because, due to the small sample size (especially for the positive ADC type), we chose to use leave-one-out CV process and, therefore, effectively all data were used to estimate the predictive performance of the classifier, as this process implies a rotation of roles (training and testing) for each sample at each iteration performed.