Functional characterization of DLK1/MEG3 locus on chromosome 14q32.2 reveals the differentiation of pituitary neuroendocrine tumors

Pituitary neuroendocrine tumors (PitNETs) represent the neoplastic proliferation of the anterior pituitary gland. Transcription factors play a key role in the differentiation of PitNETs. However, for a substantial proportion of PitNETs, the etiology is poorly understood. According to the transcription data of 172 patients, we found the imprinting disorders of the 14q32.2 region and DLK1/MEG3 locus associated with the differentiation of PitNETs. DLK1/MEG3 locus promoted somatotroph differentiation and inhibited tumor proliferation in PIT1(+) patients, furthermore, the level of DLK1 played a critical role in the trend of somatotroph or lactotroph differentiation. Anti-DLK1 monoclonal antibody blockade or siMEG3 both indicated that the DLK1/MEG3 significantly promoted the synthesis and secretion of GH/IGF-1 and inhibited cell proliferation. In addition, loss of DLK1 activated the mTOR signaling pathway in high DLK1-expressing and PIT1(+) GH3 cell lines, a mild effect in the low DLK1-expressing and PIT1(+) MMQ cell lines and no effect in the PIT1(-) ATT20 cell line. These findings emphasize that expression at the DLK1/MEG3 locus plays a key role in the differentiation of PitNETs, especially somatotroph adenomas, and provide potential molecular target data for patient stratification and treatment in the future.

Genomic imprinting, a unique mechanism of epigenetic regulation that results in parent-of-origin-specific gene expression, is essential for normal mammalian development and growth [12]. In addition to microRNAs (miRNAs), long noncoding RNAs (lncRNAs) (>200 nt in length) have attracted attention due to their potential regulatory roles in biological processes [13]. Human chromosome 14q32.2 plays crucial roles in cell differentiation and tissue development, which includes an imprinted region containing paternally expressed genes (PEGs) encode delta like non-canonical Notch ligand 1 (DLK1), retrotransposon gag like 1 (RTL1), and iodothyronine deiodinase 3 (DIO3), maternally expressed genes (MEGs) encode MEG3/8/9, and several large clusters of miRNAs [14,15]. The DLK1/MEG3 locus is mostly silenced in human nonfunctional PitNETs but not in functional PitNETs [16]. DLK1 is expressed throughout the developing pituitary gland and becomes restricted mostly to somatotroph cells within the adult gland; loss of DLK1 leads to a significant reduction in GH content throughout the lives of DLK1null mice [17]. DLK1 is involved in PIT1-mediated transcription, and GH is a DLK1-regulated target gene in the GH3 cell line [18].
Increasing evidence has suggested alterations in gene expression as major contributors to the identification of disease-specific patterns in PitNETs, not infrequent recurrent somatic mutations [19,20]. Through diversity analysis of the transcripts per million (TPM) values of the transcriptomic data and clinicopathological characteristics of 172 PitNETs, the imprinting disorders of the 14q32.2 region and DLK1/MEG3 locus associated with the differentiation of PitNETs were filtered. The DLK1/MEG3 locus functions in PitNETs, and its underlying mechanisms were largely assessed through anti-DLK1 monoclonal antibody blockade or siMEG3 in PitNET cell lines. Here, we showed the critical role of the DLK1/MEG3 locus in the differentiation of PitNETs in patients depending on the level of PIT1 and suggested the DLK1/MEG3 locus as a potential therapeutic target for somatotroph adenoma patients.

Identification of differentially expressed genes in PitNETs
In this study, we divided 172 patients into four groups according to the 2017 WHO classification of PitNETs for transcriptomic experiments; these groups included patients with somatotroph adenomas, lactotroph adenomas, corticotroph adenomas and gonadotroph adenomas. An unsupervised hierarchical clustering of the top 352 most variable genes revealed distinct gene expression profiles corresponding to the discretion of morpho-functional classification categories ( Figure 1A). The volcano maps of the characteristic difference profiles of somatotroph adenomas are shown based on pairwise comparisons ( Figure 1B-1D). A total of 110 molecules overlapped based on DEG filter parameters (|log 2 FC|>2, adjusted P value<0.001), and DLK1 (somatotroph vs. lactotroph: log 2 FC=7.397, somatotroph vs. gonadotroph: log 2 FC=9.611, somatotroph vs. corticotroph: log 2 FC=9.832) were the most significantly different in somatotroph adenomas compared with other subtypes ( Figure 1E). Using the Metascape database, DEG pathway enrichment showed that the top three pathways were Growth hormone receptor signaling (R-HSA-982772), RA biosynthesis pathway (R-HSA-5365859) and Neuroactive ligand-receptor interaction (hsa04080), and the enriched GO terms focused on negative regulation of synapse assembly, forebrain development, and regulation of catecholamine secretion ( Figure 1F, 1G). GSEA of 172 patients showed that the most positively correlated pathway was related to protein export and that the most negatively correlated pathway was related to methylated histone arginines in somatotroph adenomas compared with other subtypes ( Figure 1H).

Clinical relevance of DLK1 in somatotroph adenoma
All patients were confirmed by MR image T1WI+C sequence and transmission electron microscopy ( Figure  4A, 4B). To confirm the possibility of DLK1 can stratify PitNETs, we measured the levels of DLK1, PIT1 and Somatostatin receptor 2 (SSTR2) in these samples through immunohistochemistry staining ( Figure 4C).  Figure 4D, 4E). We divided 31 patients AGING into two groups based on the median DLK1 level: the high DLK1 group and the low DLK1 group. Patients in the high group had higher preoperative serum GH levels (P<0.001) and lower average tumor volumes ( Figure  4F). Although there was no significant difference in tumor volume between the two groups, considering that individual extreme values do exist, we believe that they still have a negative correlation trend. Furthermore, we analyzed the potent risk factors related to the serum level of GH. Based on the area under the receiver operating characteristic (ROC) curve, we identified HRAS (0.825), WNT5B (0.819), Notch2 (0.799), DLK1 (0.795) and ITPR2 (0.782) in the PI3K/AKT/mTOR pathway, Notch pathway and growth hormone synthesis, secretion and action ( Figure 4G).

The effect of anti-DLK1 antibody on PitNET cell lines
In this study, we measured the levels of PIT1 and DLK1 in the GH3, MMQ and ATT20 cell lines. Western blot experiments showed that the GH3 and MMQ cell lines were PIT1(+) cell lines and that the ATT20 cell line was a PIT1(-) cell line ( Figure 5A). Based on these results, according to the protein level, we defined GH3 cells as a high DLK1-expressing and PIT1(+) cell line and MMQ cells as a low DLK1-expressing and PIT1(+) cell line. We found that the anti-DLK-1 antibody could promote cell proliferation in a dose-and time-dependent manner in the GH3 cell line (P<0.05), and only a high concentration of the anti-DLK1 antibody (20 μg/ml) could promote cell proliferation in the MMQ cell line (P<0.05); there was no effect in the ATT20 cell line ( Figure 5B). ELISA showed that the anti-DLK1 antibody caused the strongest inhibition of GH secretion (>70%), followed by IGF-1 secretion (almost 40%), and a slight inhibition of PRL secretion ( Figure 5C). Clone forming experiments also suggested that the anti-DLK1 antibody promoted cell proliferation in the GH3 cell line ( Figure 5D). Confocal images showed a decline in PIT1(green) accompanied by the blockade of DLK1 (red) in the GH3 cell line ( Figure 5E). In addition, we found that the mTOR pathway was activated by the anti-DLK1 antibody in a dose-dependent manner in the GH3 cell line ( Figure 5F).

The effect of siMEG3 on PitNET cell lines
To measure the effect of MEG3 on the bioactivity of PitNET cell lines, we synthesized fragments of small interfering RNA. The RT-qPCR assays proved that the RNAi efficiency of the siMEG3-1 and siMEG3-2 fragments in GH3 cells was more than 75% compared with the control or vector fragments. The RT-qPCR and  AGING western blot assays both showed inhibition of the expression levels of DLK1 and PIT1 after siMEG3 treatment ( Figure 6A, 6B). Furthermore, the western blot assays showed that siMEG3 could inhibit the synthesis of GH but not of PRL in the GH3 cell line. siMEG3-1 and siMEG3-2 both promoted cell proliferation in the GH3 cell line after 48 h or 72 h of treatment (P<0.05), and only siMEG3-2 mildly promoted cell proliferation in the MMQ cell line after 72 h of treatment; however, no effect on the proliferation of the ATT20 cell line was observed, regardless of whether siMEG3-1 or siMEG3-2 was used ( Figure 6C). ELISA also showed that siMEG3-1 and siMEG3-2 both inhibited the GH/IGF-1 levels in cell culture and had no effect on the secretion of PRL or ACTH by these cell lines ( Figure 6D). Confocal assays showed a decrease in fluorescence intensity for DLK1 (red) and PIT1(green) after siMEG3-1 and siMEG3-2 treatment, especially siMEG3-2 ( Figure 6E).

DISCUSSION
Pituitary adenomas, which were recently renamed PitNETs, mostly represent benign neoplastic proliferation of the anterior pituitary gland. The 2017 WHO classification classifies PitNETs into seven morphofunctional types and three lineages. The French fivetiered classification emphasizes tumor proliferation and invasion, which provide prognostic value for the clinical treatment of PitNETs [5]. Mario et al. found gonadotroph signatures in some corticotroph and somatotroph PitNETs by integrated pangenomic analyses [4]. The low mutation burdens of PitNETs based on Whole genome sequencing (WGS) suggesting a limited impact of single nucleotide variations (SNVs) on the PitNET classifications [20][21][22]. In this study, we found that expression at the DLK1/MEG3 locus played a key role in the differentiation of PitNETs in patients depending on the level of PIT1, and we suggested that DLK1 may be used in addition to the current array of PitNET classifications.
Based on the DEGs of 172 PitNETs, DLK1 had the most obvious effect on distinguishing somatotroph adenomas from other subtype PitNETs with a high expression level over 2 7 times. AC126177.8, AC355974.2, LINC02475, MEG3, MEG9 and MiR7-3HG were the most significantly different lncRNAs in somatotroph adenomas. We confirmed imprinting disorders in the 14q32.2 region by analyzing the landscape of chromosome 14. Combined with the clinical phenotype, DLK1 and MEG3 were identified as characteristic molecules in somatotroph adenomas. The transcriptomic and immunochemistry data also suggested that the expression level of DLK1 played a key role in the differentiation of PIT1(+) PitNETs. DLK1/MEG3 is an imprinted locus consisting of multiple maternally expressed lncRNA genes and PEGs [23,24]. The expression dosage of DLK1 is functionally important, and modulation of its expression in a stem cell niche can profoundly alter the self-renewal properties of cells [25][26][27]. DLK1 plays important roles in the physiological adaptation associated with early life and is implicated in metabolic disease resistance; overexpression of DLK1 reduces fat stores, IGF-1 resistance, and defects in feedback regulation of GH in mice [28]. Ansell et al. reported GH as a target gene that is regulated by DLK1 through a PIT1-dependent mechanism [18]. These data supported our opinion that the level of DLK1 is the key factor for determining the differentiation of PIT1(+) PitNETs into somatotroph or lactotroph adenoma.
Traditionally, the DLK1/MEG3 locus plays a tumor suppressor role in human gonadotroph adenomas [16,29]. Hypermethylation of IG-DMR at the DLK1/MEG3 locus is the reason for the loss of MEG3 that occurs only in gonadotroph adenomas, and no gene in the DLK1/MEG3 locus was significantly downregulated in somatotroph adenomas compared to other PitNETs [16,30]. Overexpression of MEG3 could inhibit tumorigenesis and induce apoptosis in a pituitary-derived folliculostellate cell line [31]. It has been difficult to understand why somatotroph patients with higher serum GH levels have smaller tumor sizes and less invasive potential, while the DLK1 had the same trend of serum GH levels [32]. We speculated that the DLK1/MEG3 locus promoted somatotroph differentiation and inhibited tumor proliferation of PitNETs.

AGING
In this study, we measured the effect of DLK1/MEG3 on the bioactivity of the GH3 cell line (somatotroph), MMQ cell line (lactotroph) and ATT20 cell line (corticotroph). The western blot assays proved that the GH3 cell line was PIT1(+) with high DLK1 levels, the MMQ cell line was PIT1(+) with low DLK1 levels, and the ATT20 cell line was PIT1(-) with low DLK1 levels. Accompanied by cell proliferation, the blockade of the DLK1/MEG3 locus strongly inhibited the synthesis and secretion of GH/IGF in the GH3 cell line and mildly inhibited the synthesis and secretion of PRL, regardless of whether anti-DLK1 antibody or siMEG3 treatment was used. Blockade of the DLK1/MEG3 locus mildly increased cell proliferation but did not affect PRL secretion in the MMQ cell line. The change in the DLK1/MEG3 locus did not affect the bioactivity of the ATT20 cell line. Based on these results, we confirmed that DLK1/MEG3 promoted somatotroph differentiation and inhibited the tumorigenesis of PitNETs in PIT1mediated transcription.
The DLK1/MEG3 locus suppresses the entire PI3K/ AKT/mTOR pathway in hematopoietic stem cells by regulating mitochondrial biogenesis and metabolic activity [33]. We also noticed that DLK1 was negatively correlated with the citrate cycle in PitNET patients. This observation may explain the difference in the metabolic profile of somatotroph adenomas compared to other subtypes of patients. In our in vitro experiment, the anti-DLK1 antibody activated the mTOR pathway in GH3 cells.
In fact, the current 2017 WHO classification divides PitNETs into three transcription categories mainly based on the biological relevance of pituitary lineage factors, namely, PIT1, TPIT and SF1 [3]. However, the reliability of SF1 has been widely questioned due to variability across different research institutes, for example, in cytoplasmic and granular staining [34,35]. There is no unique molecule that specifically distinguishes somatotroph adenoma from lactotroph adenoma, even for ER-ɑ, as mentioned in the 2017 WHO classification criteria [2,36]. Our lab reported that 27/50 somatotroph adenomas patients (54%) had high ER-ɑ levels (definition: more than 50% positive cells), and only 5/42 lactotroph adenoma patients (11.9%) had high ER-ɑ levels [37]. ER-ɑ appeared to be more methylated in functioning corticotroph tumors than in functioning somatotroph adenomas (15 patients vs. 40 patients; P=0.025), which also provided evidence of high ER-ɑ levels in somatotroph adenomas [38]. The evidence above provide that SF1 and ER-ɑ are widespread in various subtypes of PitNETs. Based on the transcriptomic analysis of 172 PitNETs, we suggest a solution as follows: 1) somatotroph adenoma: high PIT1 and high DLK1; 2) lactotroph adenoma: high PIT1 and low DLK1; 3) gonadotroph adenoma: high FSHB and low PIT1/DLK1; and 4) corticotroph adenoma: high TBX19 and low PIT1/DLK1 (Figure 7). Certainly, more evidence should be provided for the exclusion of SF1 or ER-ɑ in the morphological and functional classification of PitNETs.
This study presented an adjusted classification of PitNETs based on transcription factors. The DLK1/ MEG3 imprint locus plays a critical role in somatotroph differentiation. The identification of biological mechanisms and potential clinical applications should lead to important improvements in anti-PitNET treatments. At least, combining analogs of somatostatin and blockade of the DLK1/MEG3 locus could be an effective treatment strategy for somatotroph adenomas, especially for patients with overactivation of the DLK1/MEG3 locus and resistance to standard treatment. The GH3 and MMQ cell lines (ATCC, Manassas, VA, USA) cultured in a humidified incubator at 37° C and 5% CO 2 in F-12K medium (ATCC) supplemented with 2.5% fetal bovine serum and 10% horse serum. The ATT20 cell line (ATCC) was cultured in a humidified incubator at 37° C and 5% CO 2 in DMEM (ATCC) supplemented with 10% fetal bovine serum.

RNA extractions, sequencing, transcriptomic data processing and analysis
For the RNA extractions, patient samples were processed with an AllPrep® DNA/RNA Mini kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions.
The quantity and quality of the RNA was evaluated by an RNA Nano6000 assay kit (Aligent Technologies, CA, AGING USA) (RIN>6.8). RNA (3 μg/sample) was used for the RNA preparations, and the ribosomal RNA was removed using an Epicentre Ribo-zero™ rRNA Removal Kit (Epicentre, Madison, WI, USA). The sequencing library was generated using the NEBNext® Ultra™ Directional RNA Library Prep Kit (NEB, Ispawich, USA). The library fragments (150~200 bp) were purified by the AMPure XP system (Beckman Coulter, Beverly, USA) and then assessed by the Agilent Bioanalyzer 2100 system. The libraries were sequenced on an Illumina HiSeq X platform, and then, 150-bp paired-end reads were generated. Reads containing adapters, reads containing poly-N and low-quality reads were removed. The paired-end clean reads were aligned to the human reference genome (hg19) using Hisat2 (v2.0.5) [39]. HTSeq (v0.11.2) was used to count the read numbers mapped to each gene [40]. The R package limma was used to analyze the quantitative differentiation between two identified groups [41]. KEGG pathway enrichment and GO term results were exported from the Metascape website with filtered differential gene data input [42]. The R package ClusterProfiler was used to process the GSEA analysis [43].

Tissue microarray construction and immunochemistry staining
The formalin-fixed, paraffin-embedded tissue blocks were sectioned. Three core biopsies (2.0 mm in diameter) were selected from the paraffin-embedded tissue. The cores were transferred to tissue microarrays using a semiautomated system (Aphelys MiniCore, Mitogen, UK). The microarrays were cut into 4-μm sections and incubated with anti-DLK1 (rabbit monoclonal, 1:600, ab210471, Abcam), anti-PIT1(mouse monoclonal, 1:500, sc393943, Santa Cruz) and anti-SSTR2 (rabbit monoclonal, 1:400, ab134152, Abcam) primary antibodies. BondTM Ploymer Refine Detection (Leica Biosystems, DS9800) was used for the detection of the primary antibodies. The slides were scanned into digital pictures, and expression was examined using Aperio AT2 (Leica Biosystems). The H-score was obtained by multiplying the staining intensity by a constant to adjust the mean to the strongest intensity [H-score = 3×(percentage of strong staining)] (1.0%, weak; 2.0%, moderate; 3.0%, strong) to yield a score ranging from 0 to 300.

Cell proliferation and colony formation assay
The number of cells in suspension was adjusted to 1×10 5 /ml, and 100 µl cell suspension was plated into each well of a 96-well plate. After overnight incubation, anti-DLK1 monoclonal antibody blockade (1 μg/ml, 5 μg/ml, and 20 μg/ml) or siMEG3 fragments (Supplementary Table 1) were added.
The cells were then cultured for 24 h, 48 h and 72 h; 20 µl 3-(4,5-diethylthiazol-2-yl)-5-(3-carboxymethoxyphenyl)-2-(4-sulfophenyl)-2H-tetrazolium, inner salt (MTS) solution was then added into each well and incubated for an additional 4 h. The absorbance at 490 nm was measured with an ELISA plate reader (Thermo, USA). For the colony formation assay, a total of 200 cells were plated into each well of a 24-well plate and maintained for 7 days. The colonies were then fixed with 4% formaldehyde for 30 min and stained with 0.2% crystal violet for 10 min. Triplicate wells were measured in each treatment group.

Immunofluorescence and confocal microscopy
For the immunofluorescence experiments, GH3 cells were seeded on poly-L-lysine-coated confocal dishes, and an anti-DLK1 antibody (5 μg/ml) was added. After 24 h in culture, the cells were fixed with 4% paraformaldehyde for 30 min, washed with PBS, and permeabilized with 0.2% Triton X-100 in PBS for 15 min at room temperature. The fixed cells were incubated with an anti-DLK1 (1:600) antibody and anti-PIT1 antibody (1:500) overnight at 4° C followed by incubation with the fluorescently (Alexa Fluor 488 and 594) labeled secondary antibody (1:300) for 1 h at room temperature. The cell nuclei were visualized by DAPI counterstaining (Invitrogen).

ELISA assay
The adrenocorticotropic hormone (ACTH), GH, IGF-1, and prolactin (PRL) levels in 10 μl cell culture supernatant were detected by ELISA kits (APPLYGEN) according to the protocol. The absorbance of each well at 450 nm was measured using an ELISA plate reader (Thermo Fisher).

Reverse transcription and quantitative PCR
Microarray hybridization and RT-qPCR were performed as previously described [44]. The total RNA of 30 samples was extracted and purified using the Rneasy®Mini Kit (QiaGen, Hilden, Germany) following the manufacturer's instructions. RT-qPCR was performed on a QuantStudio5 (Applied Biosystems, Singapore) with primers of each genes (Supplementary Table 2). The fold-change in the differential expression of each gene was calculated using the comparative CT method (2−∆∆CT method) in the R package PCR, with GAPDH as the reference gene [45].

Statistical analysis
All the statistical analyses were conducted using SPSS Statistics Version 22 (IBM Corporation, Armonk, New York, USA). Unpaired Student's test and chi-square (Fisher's exact) test were used to compare quantitative and qualitative data. The P value of less than 0.05 was considered significant. The R package pROC was used to compute ROC curves, AUC, and P values to evaluate the predictive accuracy of selected genes [46].

Ethics approval and consent to participate
The study protocols were approved by the Internal Review Board of Beijing Tiantan Hospital, which was affiliated to Capital Medical University, and conformed to the ethical guidelines of the Declaration of Helsinki (No. KY2016-035-01).

Availability of data and materials
All the data generated or analyzed in this study are included in this published article and its Additional files.