Differential Orthopedia Homeobox expression in pulmonary carcinoids is associated with changes in DNA methylation.

Limited number of tumor types have been examined for Orthopedia Homeobox (OTP) expression. In pulmonary carcinoids, loss of expression is a strong indicator of poor prognosis. Here, we investigated OTP expression in 37 different tumor types, and the association between OTP expression and DNA methylation levels in lung neuroendocrine neoplasms. We analyzed publicly available multi-omics data (whole-exome-, whole-genome-, RNA sequencing and Epic 850K-methylation array) of 58 typical carcinoids, 27 atypical carcinoids, 69 large cell neuroendocrine carcinoma and 51 small cell lung cancer patients and TCGA (The Cancer Genome Atlas) data of 33 tumor types. 850K-methylation analysis was cross-validated using targeted pyrosequencing on 35 carcinoids. We report bimodality of OTP expression in carcinoids (OTPhigh vs OTPlow group, likelihood-ratio test P = 1.5 × 10-2 ), with the OTPhigh group specific to pulmonary carcinoids while absent from all other cohorts analyzed. Significantly different DNA methylation levels were observed between OTPhigh and OTPlow carcinoids in 12/34 OTP infinium probes (FDR < 0.05 and β-value effect size > .2). OTPlow carcinoids harbor high DNA methylation levels as compared to OTPhigh carcinoids. OTPlow carcinoids showed a significantly worse overall survival (log-rank test P = .0052). Gene set enrichment analysis for somatically mutated genes associated with hallmarks of cancer showed robust enrichment of three hallmarks in the OTPlow group, that is, sustaining proliferative signaling, evading growth suppressor and genome instability and mutation. Together our data suggest that high OTP expression is a unique feature of pulmonary carcinoids with a favorable prognosis and that in poor prognostic patients, OTP expression is lost, most likely due to changes in DNA methylation levels.

atypical carcinoids

false discovery rate

formalin‐fixed paraffin‐embedded

fragments per kilobase of transcript per million mapped reads

glioblastoma

gene ontology

Gene‐Set Enrichment Analysis

immunohistochemistry

interquartile range

large‐cell neuroendocrine carcinoma

low‐grade glioma

lung neuroendocrine neoplasms

multi‐omics factor analysis

Orthopedia Homeobox

small‐cell lung carcinomas

typical carcinoids

The Cancer Genome Atlas

tumor node metastasis

whole‐exome sequencing

World Health Organization

INTRODUCTION

Lung neuroendocrine neoplasms (LNENs) comprise a heterogeneous group of malignancies, mainly arising from pulmonary neuroendocrine cells, which account for approximately 20% of all primary lung cancers. According to the World Health Organization (WHO) 2021 classification, LNENs are subdivided into four different entities namely, typical carcinoids (TC, 1.8%), atypical carcinoids (AC, 0.2%), large‐cell neuroendocrine carcinomas (LCNEC, 3%) and small‐cell lung carcinomas (SCLC, 15%). , Even though pulmonary carcinoids only encompass 1% to 2% of all invasive lung malignancies, their occurrence has increased significantly over the past decades. ,

LNENs share morphological and immunohistochemical features but they show a broad clinical‐pathological spectrum with a different biological behavior. In contrast to the high‐grade LNENs (eg, LCNEC and SCLC), carcinoids are characterized by a low mitotic frequency (≤10 mitosis/mm2), none or punctuated necrosis, and a relatively favorable prognosis. Although carcinoids are considered as low‐ or intermediate grade tumors, they may show distant disease recurrence (TC: 1%‐6% and AC: 14%‐29%) in patients who initially underwent curative surgical resection. ,

Known predictive factors for recurrence of disease after surgery for carcinoids are histopathological type (AC), tumor size and lymphatic involvement. Unfortunately, none of these factors enables the clinician to reliably predict patients at risk for recurrence up‐front. Therefore, extensive follow‐up is required for all patients with a carcinoid up to 15‐year after surgical resection resulting in frequent oversurveillance. , , A possible alternative is to subdivide carcinoids into prognostically relevant categories (at risk for recurrence vs no risk) using tumor‐specific molecular features. Previously, we identified Orthopedia Homeobox (OTP) as a prominent molecular marker to accurately identify patients at risk for disease recurrence and showed that OTP transcription levels and nuclear protein expression levels were strongly correlated (P = 5.9 × 10−11, n = 60). In addition, our research group and others have shown that OTP protein expression is absent in normal neuroendocrine cells of the lung, whereas OTP is highly expressed in carcinoid precursor lesions (ie, neuroendocrine cell hyperplasia and tumor lets), as well as in most pulmonary carcinoids. , , , , , Other NENs (ie, gastrointestinal/pancreatic) turned out to be OTP negative. Intriguingly, pulmonary carcinoids without expression of OTP at the time of diagnosis develop more frequently metastases during follow‐up. These data suggest that OTP is a prognostic marker in pulmonary carcinoids and may play a role in carcinoid tumorigenesis. Although OTP has been described as a key player in the development of the hypothalamic neuroendocrine system of vertebrates, its specific role in tumorigenesis and carcinoid disease recurrence remains to be elucidated.

Here, we performed a comparative analysis on OTP expression throughout 37 different tumor types, using publicly available The Cancer Genome Atlas (TCGA) data (n = 33/37 tumor types) combined with publicly available LNEN data (n = 4/37 tumor types), to identify groups with similar expression profiles. Subsequently, we used multi‐omics publicly available data on LNENs to study the molecular differences between OTPhigh vs OTPlow LNENs, with a special focus on carcinoids. In addition, we correlated transcriptomic and methylomic data to shed light on the regulatory mechanisms of OTP expression.

MATERIALS AND METHODS

Public LNEN and TCGA data

LNEN whole‐exome sequencing (WES), transcriptome sequencing (RNAseq), Epic 850K Illumina arrays methylation data, mutational data and corresponding clinical and histopathological annotations were available from the lungNENomics project within the Rare Cancers Genomics initiative http://rarecancersgenomics.com/lungnenomics/. Detailed information on data generation and quality controls can be found in Gabriel et al GigaScience 2020. This dataset includes 88 carcinoids of which 58 TC, 27 AC and 3 unclassified carcinoids, 69 LCNEC and 51 SCLC. , , , All samples were collected from surgically resected tumors.

Transcriptomic data (In units of Fragment Per Kilobase of transcript per Million mapped read, FPKM) and corresponding clinical data of 33 different tumor types were gathered from the publicly available The Cancer Genome Atlas (TCGA) platform. Data were downloaded using the Bioconductor R package TCGAbiolinks (version 2.9.5). Duplicated samples were removed as well as one lung squamous cell carcinoma sample (TCGA‐37‐4129) because this sample was previously reported as misclassified. TCGA OTP expression data was merged with LNEN OTP expression data to gather an overview of OTP expression patterns throughout 37 different tumor types.

Data processing of expression data

Data processing was performed as described previously by Alcala et al, providing quantification of expression at the gene level in two formats. (a) FPKM, one of the most popular formats for expression quantification, that facilitates comparisons across cohorts by mitigating technical batch effects through normalization based on gene length and sequencing library size; this is the format we used for comparing with TCGA cohorts, and for interpreting absolute expression levels. (b) Normalized read counts, obtained through the variance stabilization procedure described in Alcala et al using the DESeq2 R package, that facilitates downstream statistical analyses of the expression data within a cohort by reducing the relationship between mean and variance in expression level; this is the format we used for the statistical analysis of the expression (see determining an expression cut‐off point procedure, and gene set enrichment analysis below).

Clinical data

Clinical cohort data were retrieved from Alcala et al. Data included, among others, age (in years), sex (male or female), smoking status (never smoker, former smoker, passive smoker and current smoker), Union for International Cancer Control/American Joint Committee on Cancer stage (Stage I‐IV) and survival (calculated in months from surgery to last day of follow‐up or death).

Gene‐set enrichment analysis of somatic mutations

Gene‐set enrichment analysis of somatic mutations (GSEA) was performed as described previously. In short, GSEA for somatic mutations was evaluated independently for each set of genes linked to a hallmark of cancer, taking into account genes with multiple mutations, using the Fisher's exact test (Fisher.test R function, stats package version 4.0.4). P values for both the OTPhigh and OTPlow group were adjusted for multiple testing. Altered hallmarks, including the mutated genes, are presented in Table S1.

In addition, we performed robustness analyses to assess the validity of the results, especially to outlier samples/genes that would have a high influence on the statistical results, e.g., that would alone drive the significance of a particular hallmark. First, we assessed the influence of each individual sample using a jackknife procedure (ie, for each sample, we performed the GSEA test after removing this sample). Second, we assessed the influence of each gene using a jackknife procedure (ie, for each gene, we performed the GSE test without this gene). A threshold adjusted P‐value of .05 was fixed to select enriched pathway.

Methylome analyses

DNA methylation levels at each OTP Infinium probe (publicly available 850K array data) were analyzed for 51 LNEN samples (n = 19 TC, n = 15 AC, n = 17 LCNEC) of the total cohort, combined with 528 publicly available TCGA normal lung adenocarcinoma (n = 274) and squamous cell carcinoma (n = 254) tissues (both 450K array data). Differential DNA methylation was tested using the Wilcoxon rank‐sum test on β values, to allow comparisons between the nonnormally distributed distributions of β values. False discovery rate (FDR) was controlled using the Benjamini‐Hochberg procedure to correct for multiple hypothesis testing. In addition, an overview of the DNA methylation data at all OTP probes together with the genomic location of these probes was presented. The exon and transcript annotations were retrieved from Ensembl using the BioMart tool (Ensemble Genes 75, Homo sapiens genes GRCh38.p12), CpG island locations were downloaded from the UCSC genome browser, and the EPIC probe locations were retrieved from Zhou et al using the human genome version GRCh38.p12.

Validation of the methylome analyses using pyrosequencing

To cross‐validate the methylome analyses findings, we have performed pyrosequencing on n = 35/88 fresh frozen carcinoid samples of the sequencing cohort. DNA was isolated using the Gentra Puregene tissue kit 4g (Qiagen, Hilden, Germany), following the manufacturer's instructions. DNA samples were quantified by the fluorometric method (Quanti‐iT PicoGreen dsDNA Assay, Life Technologies, CA) and integrity was checked by gel electrophoresis using a 1.3% agarose gel. Next, sodium bisulfite‐modification, which converts unmethylated cytosine residues to uracil residues, was performed on 100 ng genomic DNA using the EpiTect Bisulfite Kit (Qiagen, Hilden, Germany), according to the manufacturer's instructions. Bisulfite modified DNA was amplified using methylation‐specific primers, which were designed using the Pyromark Assay design SW 2.0 (Qiagen, Hilden Germany; Table 1). PCR amplification was performed using the PyroMark PCR kit (Qiagen, Hilden Germany), according to the manufacturer's protocol, on a T100 thermal cycler (Bio‐Rad). Afterward, PCR products were loaded on a 2% agarose gel and visualized using Proxima AQ‐4 imaging software (Isogen, De Meern, The Netherlands). Next, the methylation status of cg26576712 (genomic location 77640156‐77640158) and cg02493167 (genomic location 77639893‐77639895), both located in the promoter region of OTP, were determined by pyrosequencing using PyroMark Q24 Advanced CpG reagents and a PyroMark Q24 Instrument upgraded with the PyroMark Q24 Advanced software (Qiagen, Hilden, Germany) according to the manufacturer's protocol. After pyrosequencing, the methylation percentage of each site was determined by analyzing the pyrograms using the PyroMark Q24 Advanced software (Qiagen, Hilden, Germany). Every PCR and pyrorun included a bisulfite converted and methylated DNA obtained from a RKO cell line, EpiTect unmethylated and bisulfite‐converted controls (Qiagen, Hilden, Germany) and a H2O control. In addition, an internal cytosine control was incorporated in the sequence to be analyzed. Pyrosequencing on normal lung samples (n = 3) was performed to determine the cut‐off methylation level on cg26576712 and cg02493167 (mean positivity threshold 32.2%).

TABLE 1

Overview of the primer combinations for both CpG sites within the promoter region of OTP and their specifications

Location	Primer	Sequence	Nt	T m (°C)	% CG	Product length (bp)
OTP cg02493167	Forward	GGGAGTAGTAAATATTAGTTTTTATTGTGA	30	58.8	26.7	160
	Reverse	ATTCTATACCATTTCTAATCTACTCCTAAA	30	57.3	26.7
	Pyrosequencing	ATGTTTTGTTATAAATATAATTG	23	39.2	13.0
OTP cg26576712	Forward	GTTTTTAGTTAGTATTTTTAATGTTTTGTTTAAGT	35	57.2	17.1	116
	Reverse	CCTTCCACAAAAAAATAACCCAATAA	26	58.4	30.8
	Pyrosequencing	ATGTTTTGTTTAAGTTAATTGG	22	44.6	22.7

Abbreviations: bp, base pairs; CG, cytosine‐guanine content; Nt, nucleotides; T m, melting temperature.

Immunohistochemistry

Immunohistochemistry (IHC) on formalin‐fixed paraffin‐embedded (FFPE) tissue sections was performed using the rabbit anti‐OTP primary polyclonal antibody (HPA039365; Atlas Antibodies). The primary antibody was diluted 1:3000 and incubated overnight at 4°C. Antibodies were detected by Bright Vision Poly‐HRP‐anti‐mouse/rabbit/rat immunoglobulin G (IgG; Immunologic) followed by peroxidase‐DAB (3,3′‐diaminobenzidine) visualization. Expression of OTP IHC was scored as described previously.

Statistical analysis

Statistical analyses were conducted using SPSS for Mac version 25 (SPSS Inc., Chicago, IL) and RStudio for Mac version 3.6.1. The chi‐squared test and Fisher's exact test were used to compare categorical data. One‐sample t‐test was performed to test the significance between measures in the Bland‐Altman plots. Survival analysis for LNENs was performed using Cox's proportional hazard model (R package survival version 3.2‐11), with statistical differences between groups performed using the Wald test and global fit of the model assessed using the log‐rank test. Two‐sided P values < .05 were considered significant. Bimodality in gene expression data was tested using the function normalmixEM from the R package mixtools. A likelihood ratio test was performed, using the chi‐squared test to choose between a null model (unimodal Gaussian distribution) and the bimodal model, rejecting the null when P‐value < .05. The optimal cut‐off was then defined as the lowest density point of the two Gaussian mixture distributions (in our case, 8.7).

Independent dataset of pulmonary carcinoids

To confirm our findings, we used the publicly available gene expression, mutation and Illumina 450K methylation array data from an independent pulmonary carcinoid cohort. The cohort comprised of 30 primary lung carcinoids including 17 TC and 13 AC. Methylome data were available in 17 patients (n = 10 TC and n = 7 AC). The data were downloaded from Supporting Information data files (https://www.omicsdi.org/dataset/geo/GSE118133) reported in Reference 23. Data processing and subsequent analysis were performed as previously described above.

RESULTS

Patterns of expression in pulmonary carcinoids

To gain insights into the expression patterns of OTP throughout 33 different TCGA tumor types and four LNEN subtypes, we performed comparative data analyses. Results showed that OTP was expressed in pulmonary carcinoids with higher levels in TC (median 126.4, interquartile range [IQR] 72.9‐193.4 FPKM) than in AC (median 0.16, IQR 0.06‐57.7 FPKM). However, both TC and AC groups include samples with lower and higher OTP expression, respectively (Figure 1A). All other tumor types, including pulmonary LCNEC (median 0.09, IQR 0.02‐0.17 FPKM) and SCLC (median 0.19, IQR 0.08‐0.55 FPKM), showed very low to no OTP expression. Both the glioblastoma (GBM; median 0.64, IQR 0.19‐1.7) and low‐grade glioma (LGG; median 0.007, IQR 0.001‐0.04) cohort show an overall low median expression of OTP, some samples within these groups display higher OTP expression levels (Figure 1A). Nevertheless, these data suggest that OTP is a highly specific marker for pulmonary carcinoids as compared to all other TCGA tumor types.

FIGURE 1

(A) RNA gene expression of OTP in 37 different tumor types highlighted using boxplots in fragment per kilobase million (FPKM). Center lines represent the median and box bounds represent the interquartile range (IQR). (See Table S4 for the abbreviation list of the tumor types.) (B) Histograms presenting the OTP expression pattern, in units of normalized read counts, in carcinoids (upper panel), LCNEC (middle panel) and SCLC (lower panel). The striped curve represents the distribution fit of the two Gaussian mixture distributions (component 1, in blue and component 2, in orange). The OTP cut‐off is determined as the lowest density point of the two Gaussian mixture distributions (upper panel, x = 8.7). (C) Kaplan‐Meier curves of overall survival probability for the OTPhigh and OTPlow group in pulmonary carcinoid patients [Color figure can be viewed at wileyonlinelibrary.com]

Considering that high‐ and low‐OTP expression samples were detected in both TC and AC, we assessed the OTP gene expression of 88 sequenced carcinoids to determine the optimal cut‐off point for OTPhigh‐ and OTPlow ‐classification. A mixture model of two Gaussian distributions was fitted to the distribution of OTP expression (normalized read counts) in LNENs (Figure 1B). Data showed a clear bimodal distribution of OTP in carcinoids (P = 1.5 × 10−2) suggesting that within the carcinoid cohort, two distinct groups exist, which can be separated based on OTP expression levels (Figure 1B, upper panel). Next, we confirmed the absence of OTP bi‐modality in both the LCNEC‐ and SCLC sample cohort (Figure 1B, middle and bottom panels, respectively). By applying this cut‐off point (lowest density point of the two Gaussian mixture distributions 8.7, in units of normalized read counts) to the high‐grade LNENs, only one LCNEC sample and one SCLC sample were considered OTPhigh, whereas all other high‐grade LNENs were classified as OTPlow. In fact, Alcala et al. recently showed that this LCNEC sample (S00602) clustered with carcinoid tumors according to Multi‐Omics Factor Analysis (MOFA).

Molecular and clinical characteristics of the OTP and OTP carcinoids

Baseline clinical characteristics are presented in Table 2. The OTPhigh cohort was enriched for females (69%), while the OTPlow cohort contained relatively more males (79%, Fisher's exact test P‐value = 7.9 × 10−5). Both age and smoking status were comparable within the two groups (Fisher's exact test both P‐value = .16). In addition, patients in the OTPhigh group were more frequently diagnosed as TC, whereas the OTPlow group was enriched for AC histopathology (Fisher's exact test P‐value = .00013). Compared to OTPlow, patients in the OTPhigh group were more frequently diagnosed in lower Tumor Node Metastasis (TNM) stages (Fisher's exact test of stage I‐II vs stages III‐IV P‐value = .01).

TABLE 2

Patient characteristics of the OTPhigh and OTPlow group

Variable	Groups		P‐value
	OTPhigh	OTPlow
Patients, n	64	24
Age, years
Mean ± SD	51.6 ± 18.3	58.3 ± 12.8
Median (IQR)	54 (16‐80)	58 (29‐80)	.16
Gender			7.9 × 10−5
Female	44 (68.7)	5 (20.8)
Male	20 (31.3)	19 (79.2)
Smoking status			.16
Current	14 (21.9)	3 (12.5)
Former	13 (20.3)	8 (33.3)
Never	23 (35.9)	6 (25.0)
Passive		1 (4.2)
Histopathological classification			1.30 × 10−4
Typical	50 (78.1)	8 (33.3)
Atypical	12 (18.8)	15 (62.5)
Unclassified	2 (3.1)	1 (4.2)
TNM stage			.01
I‐II	60 (94)	18 (75.0)
III‐IV	3 (5)	6 (25.0)
Unknown	1 (1)
Survival censor			.03
Alive	50 (78.1)	14 (58.3)
Death	6 (9.4)	8 (33.3)
Unknown	8 (12.5)	2 (8.3)
Median survival in months	79.3	59

Data are presented as n (%) unless otherwise stated. Abbreviations: IQR, interquartile range; TNM, tumor node metastasis.

Kaplan‐Meier survival analysis was performed for the pulmonary carcinoids clustered in OTPhigh and OTPlow. Results showed that carcinoids with low OTP expression harbor an unfavorable survival as compared to carcinoids with high OTP expression (log‐rank test P‐value = .0052, Figure 1C).

To investigate whether the survival difference between OTPhigh and OTPlow carcinoids might be associated with somatic gene mutations, we performed a GSEA for mutated genes related to hallmarks of cancer. The analysis showed that mutated genes in carcinoids clustered in the OTPlow group were strongly enriched in three hallmarks, that is, sustaining proliferative signaling (P = 1.05 × 10−3), evading growth suppressor (P = 1.11 × 10−8) and genome instability and mutation (P = 6.53 × 10−5, Figure S1, Table S1). Jackknife analyses of samples showed that OTPlow enrichments for somatic alterations were not influenced by a single sample (Table S1). In addition, jackknife analyses of mutated genes underlined that MEN1 is necessary for the enrichment of genome instability and mutation and sustaining proliferative signaling hallmarks in the OTPlow group (Table S1). In the OTPhigh cluster, three hallmark enrichments nearly reached significance including activating invasion and metastasis (P = .03), evading growth suppressor (P = .02) and sustaining proliferative signaling (P = .04, Table S1). However, it must be considered that jackknife analyses for both samples and mutated genes in the OTPhigh group revealed poor stability of enrichments, indicating that the enrichments are driven by specific samples and/or mutated genes. Interestingly, MEN1 mutations were found in 6 out of the 14 carcinoids clustered in the OTPlow group while only one case with a MEN1 mutation was found in the OTPhigh group of 40 patients. This is in agreement with several studies proving that MEN1 mutations are associated with poor prognosis. , , These data suggest that the clear and robust enrichment of hallmarks of cancers in the OTPlow group may explain tumor aggressiveness and thereby the difference in survival.

Mechanisms of OTP inactivation

Genomic OTP inactivation

To examine which underlying regulatory mechanism may cause the expression differences between OTPhigh and OTPlow carcinoids, we investigated publicly available gene mutation data. Results, to date, have shown no gene‐inactivating somatic mutations, alterations by chimeric transcripts or genomic rearrangements in the OTP gene (Supplementary data, Alcala et al ).

To further investigate the regulatory mechanism underlying OTP expression, we evaluated a possible association between DNA methylation and OTP expression using both transcriptome and methylome data (Illumina Infinium 850K for tumor samples and 450K for normal samples). Transcriptomic and methylomic data were available in 51 samples including 24 OTPhigh, 10 OTPlow and 17 LCNEC samples. Results showed 12/34 OTP Infinium probes harboring a significantly different methylation level (FDR < 0.05 and delta > 0.2) between OTPhigh and OTPlow carcinoids, of which three probes were located in the promotor region (cg02493167, cg26576712 and cg01763890; Figure 2A, Table S2). At these 12 loci, OTPhigh carcinoids harbor a lower methylation level (based on β values) as compared to OTPlow carcinoids. In addition, normal TCGA lung samples showed a baseline methylation level, which was greater than OTPhigh and lower than OTPlow carcinoids (Figure 2A, specified in green). This is interesting since OTP expression have not yet been reported in normal lung tissues. Evaluation of the LCNEC samples showed a median RNA expression level of 0.09 FPKM and high methylation levels (Figure S2). Of note, the methylation level of LCNEC samples was higher as compared to OTPlow carcinoids. Hence, we hypothesized that epigenetically hypermethylation of OTP could lead to gene‐silencing and subsequent loss of protein expression.

FIGURE 2

(A) Plot showing the DNA methylation levels at each OTP infinium probe (850K) of all carcinoid samples (OTPhigh and OTPlow) combined with TCGA normal lung adenocarcinoma and squamous cell carcinoma tissues (450K). The y‐axis on the right shows the β values; a horizontal bar was drawn at the median β‐value for each probe. Differential DNA methylation between OTPhigh and OTPlow carcinoids was calculated using the Wilcoxon rank‐sum test (significant different cg‐sites are presented in yellow). (B) Representative images illustrating OTP IHC of a pulmonary carcinoid patient showing nuclear OTP positivity and a low methylation level (left panel) and a pulmonary carcinoid patient showing absence of OTP protein expression and a high methylation percentage (right panel). (C) Heatmap of the methylation level (in β values) for the OTPhigh and OTPlow group (x‐axis) for each cg‐site (y‐axis). The cg‐sites which harbor a significantly different methylation level between the groups are presented in yellow. The upper green legend bar represents the OTP level measurement, the middle bar represents the histopathological diagnosis of each sample and the lower bar indicates the OTP group. Somatic mutations are represented in the lower rectangle [Color figure can be viewed at wileyonlinelibrary.com]

To cross‐validate the OTP methylation within the groups of pulmonary carcinoids, we performed pyrosequencing targeting two single promoter probes present in the Infinium 850K assay (cg02493167 and cg26576712). We have targeted these two CpGs since it has been frequently described that methylation within the promoter region of genes is associated with transcriptional silencing. Results revealed that 26% (n = 9/35) of patients harbored an increased methylation status while 74% (n = 26/35) showed no or low methylation percentages on cg02493167, using the methylation percentage of normal lung samples as the methylation threshold (average methylation percentage of 32.3%). Targeting cg26576712 showed an increased methylation status in 8/23 (35%) of the pulmonary carcinoids while 15/23 (65%) showed a low methylation percentage. Twelve cases failed due to technical issues. Additional IHC analysis (n = 14) confirmed loss of OTP expression in OTP‐methylated carcinoids while expression was present in unmethylated carcinoids (Figure 2B). Bland‐Altman plots showed no proportional difference between pyrosequencing and 850K arrays (Figure S3; P = .562 for cg26576712 and P = .069 for cg02493167).

The identified association between DNA methylation and OTP expression raised the question whether one or multiple CpGs are associated with expression. For this purpose, we have generated a heatmap of the β values of all cg‐sites for both carcinoid clusters (ie, OTPhigh and OTPlow) and LCNEC (Figure 2C). Complete‐linkage clustering using the Euclidean distance metric revealed a cluster with cgs that strongly correlates when applying a tree height cut‐off of three (Figure 2C, cg21907107, cg16703762, cg02086801, cg26365545, cg17374364, specified in green). These cgs might possibly explain together the regulation of OTP expression. Noteworthy is cg17374364, which shows a β‐value above .5 in all OTPlow samples (Figure 2C). While OTPhigh samples show overall a low methylation level, data reveal that some samples (eg, LNEN013, LNEN014, S01502 and S01539) tend to show a higher methylation level (β‐value > .5). Somatic mutation analysis revealed that these samples harbored mutations in genes associated with Gene Ontology (GO) terms related to histone demethylase activity (HDM that is, KDM4A, KDM5C and FBXL19) and histone methylation (HMT that is, DOT1L; Figure 2C). In addition, all cases tested with OTP IHC expression consistently showed that the OTP expression groups matched methylation groups (Figure S4).

Independent confirmation of our findings

We verified our findings using published pulmonary carcinoid data from Laddha et al, which included gene expression, mutation and Illumina 450K methylation array data of 30 samples (17 TCs and 13 ACs). Bimodality testing confirmed our previous findings showing a clear bimodal distribution of OTP expression within the group of carcinoids (Figure S5A, P = 6.9 × 10−9). After the OTP cut‐off (lowest density point of the two Gaussian mixture distributions 10.2, in units of normalized read counts), 7 patients were allocated to the OTPlow group and 23 patients to the OTPhigh group. To investigate whether OTP expression is also associated with methylation in this cohort, we have generated a heatmap using the Illumina 450K methylation array data. Results showed that the highest methylation levels, in β values, were observed in carcinoids clustered in the OTPlow group (Figure S5B). Results showed that, in line with our findings described above, all OTPlow samples harbor a β‐value above .5 on cg17374364, while OTPhigh samples show overall a low methylation level (i.e., a β‐value below .5). Furthermore, we observed a variable methylation level throughout all samples for cg05832593 (Figure S5B, specified in red), also consistent with our data (Figure 2C, specified in red).

DISCUSSION

Pulmonary carcinoids are rare lung tumors with a relatively indolent course of disease although a subgroup shows a more aggressive disease course. Previously, we have shown that low OTP gene and protein expression is associated with a poor prognosis and others have shown that OTP expression may be utilized to identify patients at risk for disease recurrence. However, thus far, the mechanisms underlying the regulation of OTP expression have not been clarified. Here, we evaluated OTP expression and methylation levels within 208 LNEN samples, 33 other tumor subtype cohorts and normal lung tissue using publicly available transcriptomics and methylomics data and identified a unique and bimodal expression of OTP in lung carcinoids. To date, no mutations or other genomic modifications (ie, chimeric transcripts and/or genomic rearrangements) have been reported in the OTP gene. Therefore, we comprehensively analyzed epigenomic data, revealing, for the first time, that differential OTP expression patterns could be explained by epigenetic modifications. Our findings were verified in 30 additional pulmonary carcinoids samples from Laddha et al.

Previously, Alcala et al correlated gene expression and promoter methylation in pulmonary carcinoids to identify genes, which expression can be explained by methylation of CpG islands. Results highlighted several top candidate genes including two homeobox genes (HNF1A and HNF4A). However, OTP, although methylated, was not among these top candidate genes. This might be the result of the fact that DNA methylation does not occur exclusively at CpG islands. Most of the tissue‐specific DNA methylation seems to occur at CpG island shores (region of lower CpG density that lies in close proximity [~2 kb] of CpG island) , , , as analyzed in‐depth here.

One of the cg sites, cg17374364, harbored a β‐value above .5 in all OTPlow samples. Nevertheless, also five of the 39 OTPhigh samples showed a β‐value above .5 (LNEN013, LNEN014, S01502, S02330 and LNEN021). Some of these samples carried mutations in lysine demethylases of histone methyltransferases, which might have impacted the methylation status of the above‐mentioned cg site.

In relation to the spectrum of LNENs (low‐grade to LCNEC), we and others have shown that downregulation of OTP expression is correlated with poor prognosis. Both a subgroup of carcinoids as well as LCNEC are OTP negative. It is interesting to investigate whether this observation matches with other observations in these tumors suggesting a temporal transition from low‐grade to high‐grade NE carcinomas. , Alcala et al and Laddha et al identified through different analyses and in a different dataset, three equivalent molecular groups of carcinoids (i.e., A1, A2 and B vs LC1, LC2 and LC3, with different clinical features). , , Importantly, in both the A1‐A2 and LC1‐LC3 groups OTP was generally highly expressed, while OTP was downregulated in the B and LC2 group. Moreover, they showed that within cluster A1 a subcohort existed, also referred to as “supra”‐carcinoids, with molecular and clinical characteristics most similar to LCNEC; these supra‐carcinoids also showed low OTP expression. Alcala et al showed that these supra‐carcinoids are, albeit their shared low OTP expression, vastly different from the carcinoid B cluster on the molecular level. Despite these distinct genomic features based on both genome‐wide expression and methylation profiles, we here observe that these supra‐carcinoids showed all high methylation levels similar to the OTPlow carcinoids. Most remarkable is LNEN021, because this supra‐carcinoid clustered as OTPhigh based on IHC analysis, but harbored high methylation levels on cg17374364 in contrast to the other OTPlow samples. This suggests that despite the high genomic difference between OTPlow and supra‐carcinoids, they share this common feature of high methylation levels on a specific CpG in the OTP gene. It remains to be investigated whether these clusters of carcinoids evolve or whether they are distinct entities from the beginning and progress as a result of punctuated tumor evolution.

The high expression levels in pulmonary carcinoids as compared to all other cancers, the correlation between loss of expression and poor prognosis and the simplicity of IHC to detect expression mark OTP as a highly suitable diagnostic and prognostic marker for daily clinical practice. However, different studies have reported additional genes which expression correlated with overall patient survival. , , , , It remains to be studied whether a panel of markers can further improve the prognostication of pulmonary carcinoids.

Our study implicates DNA methylation as possible regulatory mechanism of OTP; nevertheless, our study has several limitations. First, it is a retrospective study with a rather small cohort. Albeit pulmonary carcinoid is an orphan disease, future studies should validate our findings using larger cohorts containing extensive clinical data (ie, relapse‐free survival). Second, our study did not contain data on the methylation level of normal neuroendocrine cells of the lung. It would be very informative to investigate the methylation levels in these specific cells and in other neuroendocrine tumor types in future studies. Furthermore, we should investigate whether methylation may occur as the result of temporal tumor evolution by analyzing the methylation levels in tissue of primary carcinoid tumors at diagnosis that later developed disease recurrence as well as the corresponding metastatic lesion. Third, even though we did not observe significant differential expression levels in the TET (TET1, TET2 and TET3) and DNMT family (DNMT1, DNMT3A and DNTM3B), two enzyme families that play a major role in DNA methylation, between OTPhigh and OTPlow carcinoids, future studies should investigate the members of the methylation machinery in more detail (Table S3). Finally, future studies might investigate the role of miRNAs as additional regulators of OTP expression.

CONCLUSION

We show that within the group of pulmonary carcinoid patients, two distinct groups exist which can be separated based on OTP expression. To the best of our knowledge, we are the first to prove that the differential OTP expression within pulmonary carcinoids is associated with changes in DNA methylation. These findings arouse curiosity about whether epigenetic therapies might be useful for pulmonary carcinoid patients in the future.

CONFLICT OF INTEREST

Dr Speel reports personal fees from Amgen, Eli Lilly and Novartis, outside the submitted work. Nonfinancial support from Biocartis and Abbvie and grants from AstraZeneca, Pfizer, Novartis and Bayer, outside the submitted work. Dr Dingemans reports financial conflicts (ad board, paid to institute) of Eli Lilly, Roche, Amgen, Jansen, Pfizer, Astra‐Zeneca, Bayer, Boehringer, Ingelheim, Sanofi, outside the submitted work. Dr van Engeland is co‐founder and shareholder of Epify BV and MLA Diagnostics BV, companies which develop cancer DNA methylation markers. All other authors declare no potential conflicts of interest. Where authors are identified as personnel of the International Agency for Research on Cancer/World Health Organization, the authors alone are responsible for the views expressed in this article and they do not necessarily represent the decisions, policy or views of the International Agency for Research on Cancer/World Health Organization.

ETHICS STATEMENT

All new specimens were collected from surgically resected tumors, applying local regulations and rules at the collecting site, and including patient consent for molecular analyses as well as collection of de‐identified data, with approval of the IARC Ethics Committee.

Appendix S1 Supporting Information

Click here for additional data file.

Table S1 Hallmarks of cancer gene set enrichment analysis

Click here for additional data file.

Table S2 Methylation overview for each OTP cg site

Click here for additional data file.

Table S3 Differential gene expression analysis between OTPhigh and OTPlow carcinoids for genes of the TET and DNMT family

Click here for additional data file.

Table S4 TCGA abbreviations

Click here for additional data file.