Global gene expression of histologically normal primary skin cells from BCNS subjects reveals “single-hit” effects that are influenced by rapamycin

Studies of dominantly heritable cancers enabled insights about tumor progression. BCNS is a dominantly inherited disorder that is characterized by developmental abnormalities and postnatal neoplasms, principally BCCs. We performed an exploratory gene expression profiling of primary cell cultures derived from clinically unaffected skin biopsies of BCNS gene-carriers (PTCH1 +/-) and normal individuals. PCA and HC of untreated keratinocytes or fibroblasts failed to clearly distinguish BCNS samples from controls. These results are presumably due to the common suppression of canonical HH signaling in vitro. We then used a relaxed threshold (p-value <0.05, no FDR cut-off; FC 1.3) that identified a total of 585 and 857 genes differentially expressed in BCNS keratinocytes and fibroblasts samples, respectively. A GSEA identified pancreatic β cell hallmark and mTOR signaling genes in BCNS keratinocytes, whereas analyses of BCNS fibroblasts identified gene signatures regulating pluripotency of stem cells, including WNT pathway. Significantly, rapamycin treatment (FDR<0.05), affected a total of 1411 and 4959 genes in BCNS keratinocytes and BCNS fibroblasts, respectively. In contrast, rapamycin treatment affected a total of 3214 and 4797 genes in normal keratinocytes and normal fibroblasts, respectively. The differential response of BCNS cells to rapamycin involved 599 and 1463 unique probe sets in keratinocytes and fibroblasts, respectively. An IPA of these genes in the presence of rapamycin pointed to hepatic fibrosis/stellate cell activation, and HIPPO signaling in BCNS keratinocytes, whereas mitochondrial dysfunction and AGRN expression were uniquely enriched in BCNS fibroblasts. The gene expression changes seen here are likely involved in the etiology of BCCs and they may represent biomarkers/targets for early intervention.


INTRODUCTION
Basal cell carcinomas (BCCs) of the skin represent nearly one half of all cancers diagnosed in the United States [1,2].Although these tumors rarely metastasize and are an infrequent cause of cancer mortality [3], they are associated with significant morbidity due to local invasion and tissue destruction.The basal cell nevus syndrome (BCNS), also known as Gorlin syndrome or the nevoid basal cell carcinoma syndrome, is an autosomal dominant disorder characterized by multiple BCCs, occasional childhood malignancies, and developmental defects including palmar pits, jaw cysts, coarse facial appearance, malformations of the ribs, spine and brain, macrocephaly, and generalized overgrowth [4][5][6].
The majority of BCNS cases can be attributed to heterozygous germline mutations in PTCH1; the human homolog of the Drosophila gene, patched [7].Patched (PTCH1) encodes the receptor for hedgehog (HH), a secreted morphogen involved in establishing the basic framework of developing embryos in Drosophila and in other organisms, including vertebrates [8,9].In the canonical pathway, binding of HH to PTCH releases inhibition on SMO, which then signals through a series of steps to zinc finger trans-activators that ultimately convey the signal to the nucleus.In humans, the transactivator function is performed by GLI1, GLI2 and GLI3.Genes directly regulated by GLI include PTCH1, HIP (hedgehog interacting protein), and GLI1 itself, in many tissues.A host of other potential target genes in normal cells have been identified by expression microarrays where it appears that HH signaling is influenced by the stage of development and tissue-specific factors [8,9], leading to different global expression patterns in different organs [10].
Appropriately, tumor suppressor genes and their effector pathways have been identified in several dominantly heritable cancers, providing insights about cancer initiation and potential intervention [11][12][13].Analogous to retinoblastoma [11,12] and other autosomal dominant forms of heritable cancer [13][14][15][16][17][18], BCCs in BCNS patients, arise through a two-hit mechanism in which "one hit" is an inherited, inactivating mutation in PTCH1, and the second hit is a somatically derived mutation in the remaining PTCH1 allele [7].The great majority of sporadic BCCs have two somatically derived mutations in PTCH1, or less often, an activating mutation in SMO that mimics loss of PTCH1 [7,[19][20][21][22][23].Although additional molecular changes are probably necessary for the development of BCCs, it is likely that loss of PTCH1 function and consequent activation of the hedgehog pathway is a necessary early step [24,25].
This "one-hit model" for BCNS is supported by many features of this disorder such as developmental defects that have a generalized or symmetric manifestation which would be difficult to reconcile with a two-hit mechanism.Generalized overgrowth, acromegalic bone structure, and relative macrocephaly, for example, suggest that heterozygous loss of function leads to cellular hyperproliferation in bone or cartilage [26].In this regard, expression analysis of BCNS fibroblasts and fibroblasts from unaffected individuals [27] implicated genes that were hypothesized to contribute to a growth advantage in BCNS fibroblasts.
BCNS patients can develop thousands of skin tumors over their lifetimes, and treating this multiplicity of tumors with conventional surgical or other ablative therapies is a challenge.Preventive or curative medical therapy with an effective targeted agent that is known to have limited adverse, off-target effects would represent a great advance.Although several compounds are known to inhibit hedgehog signaling [28][29][30], their overall preventative efficacy of BCC remains to be established, and side effects limit their long-term use.In this regard, it is important to note that cell identity switch allows residual BCC to survive hedgehog pathway inhibition, necessitating the use of an additional agent(s) [31].
Rapamycin, also known as sirolimus, is an antibiotic isolated from Streptomyces hygroscopicus [32] that blocks mitogenic signaling by growth factor receptors and is a promising anti-neoplastic agent for a variety of tumors [33][34][35][36][37].The best studied pathway affected by rapamycin is that conducted by the mTOR complex 1 (mTORC1), which consists of the mTOR protein (a member of the (PI(3)K) family) DEPTOR, PRAS40, RAPTOR, mLST8, and TTI1-TEL2 [33][34][35][36][37]. Importantly, rapamycin has also been shown to have an effect on the canonical HH pathway wherein crosstalk between mTOR and HH pathways has been indicated [38][39][40].Thus, while blocking mTOR may take away a signaling component that is necessary but not sufficient for transformation, the mTOR and HH pathways coincidentally may target an overlapping set of genes that are likely to inhibit BCC carcinogenesis.Incidentally, rapamycin has been shown to inhibit BCC development in kidney transplant patients [41,42] and has also been shown in a case report to inhibit BCCs, advocating its use when surgical intervention is counter-indicated [43].
The first aim of this study was to profile baseline global gene expression in cells obtained from histologically normal skin of BCNS subjects and compare these with skin cells from unaffected individuals to detect "single-hit" effects caused by PTCH1 (+/-).Since hedgehog activation has been proposed to function in carcinogenesis through paracrine activity on stroma as well as cell autonomous activity [44,45], both keratinocytes and fibroblasts were analyzed in this study.Our second aim was to determine if a particular set of genes was differentially influenced by rapamycin in tissues that might have been "sensitized" to its effect by loss of one copy of PTCH1.www.oncotarget.com

Exploratory analysis at baseline reveals no clear association of samples based on PTCH1 mutation status
The first goal of this study was to explore changes in gene expression that can be attributed to mutations in PTCH1 and can possibly contribute to molecular abnormalities of basal cell carcinomas.Here we profiled gene expression patterns of PTCH1 (+/-) keratinocytes and PTCH1 (+/-) fibroblasts derived from histologically normal tissues of BCNS subjects and compared them to cells derived from normal individuals at baseline and following treatment with rapamycin (Table 1).A principal components analysis (PCA) (Supplementary Figure 1) revealed that although the two treatments clustered separately, i.e., with and without rapamycin, the differences in gene expression due to PTCH1 mutation did not result in a significant clustering and separation of BCNS from normal samples, indicating that the PTCH1 mutation itself, whether truncating or missense, does not alter global gene expression of keratinocytes or fibroblasts enough to provide clear separation between BCNS and normal specimens.The lack of clear delineation of samples at baseline can be attributed to an attenuated canonical HH pathway in PTCH1 (+/-) cells grown in culture [46,47].However, rapamycin treatment caused significant changes in global gene expression pattern that provided a clear separation of treated and untreated samples derived from either BCNS patients or normal subjects (Figure 1).The PCA results were further supported by unsupervised hierarchical clustering of genes (coefficient of variation of more than 0.1 in the microarray data) at baseline and after rapamycin treatment (Figure 1).While control samples distinctly cluster apart from rapamycin treated samples for both the keratinocytes (Figure 1A) and fibroblasts (Figure 1B), the clustering pattern was independent of PTCH1 mutation status or rapamycin concentration.Together, these results reaffirm our observation that mono-allelic PTCH1 mutation did not affect the global gene expression pattern to distinguish BCNS from normal samples in vitro.Secondly, the analysis suggests that the effect of rapamycin at a more stringent level of statistical scrutiny (FDR<0.05) is dependent on the intrinsic gene expression pattern unique to the sample derived from each individual.

Differentially regulated genes in BCNS and normal samples at baseline under more relaxed statistical constrains
Following 5-way mixed-model ANOVA-REML (FC >1.3 and <-1.3; unadjusted p<0.05), we found several gene sets that were differentially expressed.Briefly, the baseline comparison of microarray data of keratinocytes and fibroblasts showed 585 and 857, respectively of differentially expressed genes between BCNS and normal subjects.Hence,

GSEA specifies baseline differences between BCNS and normal keratinocytes or fibroblasts
Next, we performed baseline gene set enrichment analyses (GSEA) of both keratinocytes and fibroblasts samples derived from BCNS patients and compared them to the control group (Figure 2).The results indicated a positive correlation of BCNS keratinocytes with predefined pancreatic β cell hallmark gene set that comprises 40 genes, specifically upregulated in pancreatic beta cells and known to be involved in glucose metabolism and dysglycemia (Figure 2A).Interestingly, BCNS keratinocytes, in contrast to normal keratinocytes, correlated negatively with the oncogenic gene signature mTOR_UP.N4.V1_DN in the MSigDb (Figure 2A and Supplementary File 2), comprising genes shown to be down-regulated upon treatment with rapamycin [48].In the case of BCNS fibroblasts, compared with normal fibroblasts, GSEA pointed to a negative correlation with three oncogenic signatures in the MsigDb that comprised genes down-regulated in response to activated CTNNB1 overexpression (BCAT_BILD_ET_ AL_DN) and RNAi-mediated JAK2 knockdown (JAK2_ DN.V1_DN), and genes up-regulated upon RNAi-mediated knockdown of PCGF2 (MEL18_DN.V1_UP) (Figure 2B).

Biological differences attributed to rapamycin effects in BCNS and normal keratinocytes
Here we used GSEA algorithm on comprehensive microarray dataset without preprocessing to determine an a priori defined gene sets that exhibit statistically significant and biologically relevant differences between untreated and rapamycin treated BCNS keratinocytes and fibroblasts.The GSEA results indicated several significant pathways that were altered by rapamycin treatment in both BCNS keratinocytes and BCNS fibroblasts (Figure 3, Supplementary Files 2-4).
In BCNS keratinocytes, enriched signatures in the presence of rapamycin signal the involvement of and regulation by MYC and G2/M checkpoint genes, cell cycle related targets of E2F transcription factor and mitotic spindle assembly [49].The oncogenic signatures enriched in BCNS keratinocyte treated with rapamycin showed sets of genes up-regulated in neuronal precursors after Shh stimulation, PIGF treatment, genes governed by overexpression of E2F1 and E2F3, and genes upregulated in primary keratinocytes from RB1 and RBL1 skin-specific knockout mice.Interestingly as well, the YAP conserved signature [49][50][51], was also shown to be concordant with rapamycin treated BCNS keratinocytes  3A, Supplementary Files 2-4).The oncogenic signatures comprising of down-regulated genes leading to CTNNB1 overexpression, SRC overexpression, oncogenic KRAS expression in context of TBK1 knockdown, RNAi mediated knockdown of genes such as EIF4G1, HOXA9 and RPS14, and genes down-regulated after VEGFA treatment, positively correlated with BCNS keratinocytes after rapamycin treatment (data not shown).

Biological differences attributed to rapamycin effects in BCNS and normal fibroblasts
The GSEA pointed to genes upregulated by activation of PI3K/AKT/mTOR pathway and mTORC1 complex, as well as genes involved in mitotic spindle assembly, DNA repair, G2/M checkpoint, and genes regulated by MYC which appear to be in agreement with (A) GSEA shows that the 'HALLMARK_PANCREAS_BETA_CELLS' gene-set comprising of genes specifically up-regulated in pancreatic beta cells is enriched in BCNS keratinocytes as compared to the normal keratinocyte samples.The GSEA analysis also indicates that the 'MTOR_UP.N4.V1_DN' gene set containing down-regulated genes post rapamycin treatment is enriched in the normal keratinocytes (A).(B) For the fibroblasts, GSEA shows that the 'BCAT_BILD_ET_AL_DN' gene set comprising of genes down-regulated due to activated CTNNB1 overexpression; the 'MEL18_DN.V1_UP' gene set comprising of up-regulated genes in medullablastoma cells after PCGF2 knockdown; and the 'JAK2_DN.V1_DN' gene set containing the genes down-regulated by JAK2 knock-down are all enriched in the normal fibroblast samples.For the enrichment plots, profile of running ES score and position of GeneSet members on the rank ordered list are denoted for HALLMARK_PANCREAS_BETA_CELLS (A), MTOR_UP.N4.V1_DN (A), BCAT_BILD_ET_AL_DN (B), MEL18_ DN.V1_UP (B), JAK2_DN.V1_DN (B), and the respective normalized enrichment scores (NES), nominal p-values (NOM p-val), and q-values (FDR q-val) are summarized in table below the plots.the transcriptomic profiles of BCNS fibroblasts post rapamycin treatment (Figure 3B, Supplementary Files 3-4).Additionally, genes upregulated by APC knockout, RNAi mediated ELK3 knockdown in response to hypoxia, MYC and IRF4 targets, along with genes down regulated in lesional skin biopsies from mycosis fungoides patients compared to normal skin samples positively correlated with the results on treated BCNS fibroblasts.Oncogenic signatures comprising of up-regulated genes in primary keratinocytes from RB1 and RBL2 skin specific knockout, genes upregulated upon EED or EZH2 knockdown, MYC overexpression, stimulation with Shh, and genes upregulated by everolimus are found to be correlated with treated BCNS fibroblasts.Also, significant immunologic signatures were enriched in the treated BCNS fibroblasts.Furthermore, GSEA indicated genes involved in glycogenesis and gluconeogenesis and genes upregulated by PI3K/AKT/mTOR pathway and mTORC1 complex activation, genes important for and involved in mitotic spindle assembly, DNA repair, G2/M checkpoint, genes responding to estrogen and genes regulated by MYC [49] were negatively correlated with the transcriptomic profiles of BCNS fibroblasts post rapamycin treatment.Additionally, genes expressed in multipotent progenitors, genes upregulated by APC or BCL3 knockout or RNAi mediated ELK3 knockdown in response to hypoxia, MYC and IRF4 targets, along with genes downregulated in lesional skin biopsies from mycosis fungoides patients compared to normal skin samples and genes upregulated in peripheral blood monocytes of Sezary syndrome, were negatively correlated with rapamycin-treated BCNS fibroblasts.Also, oncogenic signatures comprising of YAP [50,51] conserved signature, up-regulated genes in RB1 and RBL2 skin specific knockout primary keratinocytes, genes upregulated upon EED or EZH2 knockdown, MYC overexpression, stimulation with Shh, and genes upregulated by everolimus, were found to be negatively correlated with treated BCNS fibroblasts.Gene sets consisting of downregulated genes upon SRC overexpression, RPS14 knockdown, HOXA9 knockdown, and VEGFA treatment were also negatively correlated with rapamycin treated BCNS fibroblasts.On other hand, significant immunologic signatures were enriched in untreated BCNS fibroblasts.

Ingenuity Pathway Analysis (IPA) shows distinct signaling pathways enriched in BCNS samples treated with rapamycin
To study whether rapamycin treatment affected cells derived from BCNS patients with PTCH1 mutation differently than normal samples, we generated a gene list of the top differentially regulated and statistically

Gene Symbol
Gene  3  and 4 (the entire gene list is available in Supplementary File 4).Thus, a total of 4797 genes were differentially expressed in normal fibroblasts and expression of 3214 genes were significantly altered in normal keratinocytes after rapamycin treatment.In the BCNS group, a total of 4959 genes and 1411 genes were differentially expressed in fibroblasts and keratinocytes, respectively after treatment with rapamycin.To tease out genes that were usually altered in both normal and BCNS groups following rapamycin treatment, we used Venn diagrams (Figure 4) to generate gene lists for further analysis.

IPA core analysis indicates effect of rapamycin treatment on keratinocyte and fibroblasts
Here we explored the genes/signaling pathways differentially expressed in response to rapamycin.To this end, we performed core IPA analysis on comprehensive   5B, Table 5B).
In the case of the normal fibroblasts, an IPA core analysis of 4797 differentially regulated genes in response to rapamycin resulted in enrichment of role of BRCA1 in DNA damage response, role of CHK proteins in cell cycle checkpoint control, mitotic roles of Polo-like Kinases, cell cycle control of chromosomal replication and protein ubiquitination pathway as canonical signaling pathways (Figure 5C, Table 5C).On the other hand, a core analysis of 4959 genes enriched in rapamycin treated BCNS fibroblasts pointed to top canonical pathways that included those playing a role of BRCA1 in DNA damage response, protein ubiquitination pathway, hereditary breast cancer signaling, cell cycle control of chromosomal replication and mitotic roles of Polo-like Kinases (Figure 5D, Table 5D).Although it may appear that there was some overlap among canonical signaling pathways between normal and BCNS samples, the focus molecules in these pathways In order to understand how canonical pathways that are uniquely enriched in normal or BCNS specimens respond differently to rapamycin due to differentially regulated genes, we dissected the gene list guided by the Venn diagram (Figure 4) and performed core analysis on each group.The canonical pathways enriched in both normal and BCNS keratinocytes, based on 812 common genes, were p53 signaling, DNA damage induced 14-3-3 signaling, EGF signaling, IL-3 signaling and PTEN signaling.The IPA core analysis of the 812 differentially regulated genes common to both BCNS and normal keratinocytes after rapamycin treatment predicted CDKN1A, KDM5B, let-7 and PTEN to be activated   Hepatic fibrosis/hepatic stellate cell activation, HIPPO signaling, natural killer cell signaling, galactose degradation (Leloir Pathway) and ILK signaling were identified as the top canonical pathways unique to the 599 focus genes in the BCNS keratinocytes samples after rapamycin treatment.The top up-regulated molecules were BIRC3, ZBTB10, IFITM1, GART, CDK19, PLEKHM3, RSF1, FN1, ZCCHC7, TGFBR3.The top down-regulated molecules were KRT1, CLCA2, ARTN, TPM4, VPS53, S100A8, ACOT11, HDHD2, FKBP5, GNG8.The upstream regulators inhibited in the BCNS keratinocytes samples were microRNAs miR-1-3p (no activation score), miR-381-3p, miR-19b-3p, miR-19b-3p, miR-374b-5p, miR-340-5p.JNK signaling and miR-101-3p regulator network was predicted to be associated with the signaling response of BCNS keratinocytes to rapamycin.
Similar analysis of fibroblast samples that was based on 3496 differentially expressed genes common to both normal and BCNS fibroblasts resulted in enrichment of canonical pathways involving cell cycle control of chromosomal replication, role of BRCA1 in DNA damage response, mitotic roles of Polo-like kinase, hereditary breast cancer signaling and role of CHK proteins in Cell cycle checkpoint control.The top up-regulated molecules are ADH1B, ITGB8, DIO2, EFEMP1, GRIA1, SLC40A1, RGCC, FMO2, SVEP1, CCL2, GLDN, RGCC, STEAP4, DAPK1.The top down-regulated molecules were ANLN, PBK, DLGAP5, RRM2, HMMR, TOP2A, SHCBP1, CEP55, TTK, NDC80, KIF20A, CCNB1, BIRC5.In both normal and BCNS fibroblasts samples, upstream regulators CDKN1A and TP53 were activated whereas the E2F1 and RABL6 were predicted to be inhibited in response to rapamycin.It is of interest to note that RABL6 plays a role in cell growth and survival and its overexpression is associated with tumorigenesis.

Networks enriched in BCNS samples after rapamycin treatment
Here we performed an IPA network analysis that revealed several important such networks that are significantly altered post-rapamycin treatment (Figures 6 and 7).The most notably altered pathway in BCNS keratinocytes compared with normal keratinocytes, based on 35 focus molecules (score 50), was linked to digestive system development and function (Figure 6A and 7A).Several genes such as APP, SPON1 (cell adhesion promotes neurite growth and axonal guidance), CISD2 (associated with neurodegenerative Wolfram syndrome 2), KLK9 (associated with collagen biosynthesis and modification), CDCA3 (involved in cell cycle regulation mediated by APC) were included in this network and are related to development, neuronal growth and function.
Similar analyses of BCNS and normal fibroblasts, based on 35 focus molecules, showed enrichment of networks related to cell cycle, DNA replication, recombination and repair (Figure 6B and 7B).Notably, CASC5 (kinetochore-microtubule attachment and chromosome segregation, MAPK1 (involved in proliferation, differentiation, transcriptional regulation), and CROCC (centrosome cohesion before mitosis, associated with medulloblastoma) were included in this network.

DISCUSSION
The cellular mechanisms whereby hedgehog (HH) pathway activation leads to BCCs are yet to be established [29].Excess HH signaling may directly cause proliferation of BCC precursors [9,52].Alternatively, the HH pathway may be involved primarily in stem cell maintenance as has been shown in tumor stem cells for some tissues [8,53,54].As either a direct effector of proliferation or maintainer of stem cells, HH would be acting in a cell autonomous role in BCCs.An alternative model involves a paracrine function in which HH signaling acts upon stroma to generate a tumor-promoting microenvironment while having a modest role in the tumor cells themselves [44,45].This model is appealing in the case of BCCs because of evidence that conditioned stroma is necessary for maintenance of tumor viability and growth.There are data supporting all three of these biological roles in cancer, and it is not clear which may apply to epidermal cell transformation [29].Likewise, the fine details of the molecular mechanisms by which HH pathway activation exerts its effects are not entirely understood.Most investigations have focused on downstream targets of the GLI transactivators.A number of specific candidate genes have been studied in human and mouse BCCs including PDGF receptor, BCL2, cyclin D1 and D2.However, unbiased attempts by global gene expression to determine key HH target genes important in BCCs have yielded inconsistent results (reviewed in [22]).
Using exploratory analysis at baseline, we were unable to identify gene expression changes that molecularly distinguish affected individuals from controls.The lack of clear delineation of BCNS subjects from normal controls is depicted by PCA and HC analyses of fibroblasts or keratinocytes (Supplementary Figure 1 and Figure 1).Our results showing lack of BCNS sample clustering on the basis of mutation type or gender are in line with previous observations demonstrating lack of evidence for genotype-phenotype correlation in BCNS cases, including no association between mutation type, age of onset, or number of BCCs developed in BCNS patients [55].Conceivably, subtle gene expression changes caused by single copy loss of PTCH1 may enable the intact PTCH1 gene to exert negative feedback control, particularly under in vitro growth conditions.Along these lines, previous studies showed that the up-regulated HH-signaling in vivo was suppressed in BCC and medulloblastoma (MB) cells grown in vitro [46,47].Moreover, HH signaling of MB cells grown in vitro was not restored upon transplantation in nude mice, nor did they respond to treatment with HH antagonists [47].It is also interesting to consider whether apart from the classical tumor suppressor 'two-hit' model, the 'continuum model' [56] that accounts for subtle dosage effects, 'obligate haploinsufficiency', is associated with the effect of PTCH1 tumor suppressor gene in BCNS.
Hence, our results at low stringency (p-value <0.05, no FDR cut-off; FC 1.3) have enabled us to unmask candidate baseline differences between BCNS and normal control cells that are likely to play a role in BCC growth under conditions in which the HH pathway has been suppressed [47].Thus, a total of 585 and 857 genes differentially expressed in keratinocytes and fibroblasts samples from BCNS subjects, respectively were identified.For example, a Gene Set Enrichment Analysis (GSEA) of highly regulated candidate genes recognized pancreatic β cell hallmark genes and mTOR signaling genes in BCNS keratinocytes, whereas analyses of BCNS fibroblasts candidate genes showed gene signatures that affect pathways regulating pluripotency of stem cells, including the WNT pathway.Moreover, these differentially affected gene signatures included canonical and non-canonical HH pathways in BCNS keratinocytes and fibroblasts that are likely to play a role in BCC development.One such example that is likely relevant in this context is our finding that MKL1 is up-regulated in BCNS keratinocytes, since its activation within the non-canonical HH pathway has recently been shown to promote drug resistance of BCCs [57].
Importantly, the phenotypic features of BCNS patients shown in our study and elsewhere [4,[58][59][60][61], provide strong support for a two-hit mechanism wherein the first hit is associated with early molecular changes and a variety of developmental defects occurring concurrently.While some of these defects may in fact arise through a two-hit mechanism, it is difficult to imagine how such symmetric features as generalized overgrowth, macrocephaly, coarse facies, and hypertelorism could arise through any mechanism other than one invoking an effect of haploinsufficiency.Hence, gene expression patterns in tissues at some stages of development must almost certainly differ between BCNS subjects and controls.Accordingly, networks of differentially regulated genes at baseline shown here are of interest as they were found to be enriched in processes which might play important roles in the development of congenital phenotype (e.g., macrocephaly, skeletal abnormalities, bossing of forehead).For example, BEX1 (Brain Expressed, X-Linked 1) plays a role in cell cycle progression and neuronal differentiation, and is related to p75NTR/ NGFR-mediated signaling and it inhibits neuronal differentiation in response to nerve growth factor (NGF).SOX4 (SRY-related HMG-box) is a member of SOX family of transcription factors that regulate pathways in www.oncotarget.comembryonic development and cell fate determination and is related to ERK signaling.MEG3 (Maternally Expressed 3) may function as a lncRNA tumor suppressor as it inhibit tumor cell proliferation, interacts and regulates p53 gene expression and its down-regulation has been experimentally observed in cancer cell lines of various origins.SYN1 (Synapsin 1) is member of neuronal phosphoproteins that function to regulate axonogenesis and synaptogenesis.SIM2 (Drosophila single-minded gene homolog) encodes for transcription factor that regulates neurogenesis and may have pleiotropic effects in tissues during development.Thus, perturbations caused by loss of function PTCH1 mutations although subtle, seem sufficient to result in developmental anomalies and BCCs in BCNS patients.The occurrence of PTCH1 mutations that result in developmental anomalies concomitant with BCC expression serve as poignant reminder of the role of PTCH in normal tissue development and as a tumor suppressor gene (TSG), a basic tenet in cancer biology, including all TSGs identified to date [62].
Prior studies suggest that signaling through the hedgehog and mTOR pathways have separate, albeit synergistic effects on malignant transformation [33][34][35]37].Indeed, our present GSEA studies (vide supra) showed an activated mTOR pathway in BCNS PTCH1 (+/-) keratinocytes.Accordingly, mTOR/S6K1 and HH pathways coincidentally target an overlapping set of genes that are key modulators of HH-related carcinogenesis [38][39][40].In this regard, a recent study has demonstrated the inhibition of human rhabdomyosarcoma xenografts by rapamycin through its effect on HH effector genes such as Gli1 and Gli 2 as well as PTCH1 and PATCH2 [63].These effects were associated with over 80% reduction in cyclin D1, a downstream transcription target for both HH and mTOR signaling pathways [63].
To investigate the differential effect of rapamycin in BCNS, we compared the gene expression profiles of keratinocytes and fibroblasts derived from normal subjects and from unaffected skin cells of BCNS patients.Accordingly, rapamycin affected the expression of 1411 genes and 3214 genes which were altered in BCNS and normal keratinocytes, respectively (FC >1.3 and <-1.3;FDR <0.05), including 812 genes whose altered expression was shared by both group.A similar analysis of BCNS and normal fibroblasts resulted in 4959 genes and 4797 genes, respectively wherein expression in each group was altered following rapamycin treatment, including expression of 3496 genes that was common in both groups.The increased stringency in analyzing the rapamycin effect as compared with baseline's lesser stringency was meant to focus on the most salient alterations seen due to rapamycin.We focused on these gene lists via an IPA analysis in order to discern differences in molecular functions and canonical signaling.For example, our results show that rapamycin treatment of BCNS keratinocytes affected Wnt/β-catenin Signaling, ILK signaling, epithelial adherens junction signaling, Sonic Hedgehog signaling, IGF-1 signaling, caveolar-mediated endocytosis signaling, and amyloid processing (Supplementary Data).On the other hand, pathways that were associated with response to rapamycin in BCNS fibroblasts include paxillin signaling, DNA damageinduced 13-3-3σ, integrin signaling, Cdc42 signaling, telomere extension by telomerase, HGF and VEGF signaling.The overlap between hedgehog targets and rapamycin targets provides evidence for a model in which the anti-neoplastic effect of rapamycin relates to suppression of HH signaling.Rapamycin is the subject of approximately 1500 clinical trials for cancer treatment and prevention listed on the NIH Clinical Trials web site (clinicaltrials.gov).Among these trials, five ongoing or recently completed studies relate to the use of this compound in prevention of skin cancer in transplant recipients.
This study underscores the need to catalog baseline molecular differences in BCNS cohorts in order to anticipate the numerous manifestations characteristic of BCNS, particularly the locally invasive but seldom metastatic nature of BCCs.Perturbations caused by loss of function of PTCH1 mutations unmasked herein, including alterations in canonical and noncanonical HH pathways, although subtle, seem sufficient to result in developmental anomalies and BCCs in BCNS patients.This is the first study involving skin derived keratinocytes and fibroblasts from BCNS patients wherein global gene profiling may be used to generate rationale hypotheses and design of future experiments, including noncanonical HH pathways.

Ethics statement
The clinical protocol under which this study (https:// clinicaltrials.gov/ct2/show/record/NCT00433485)was conducted was approved by Yale University's Institutional Review Board, the Clinical Trials Protocol Review Committee of the Yale Comprehensive Cancer Center, and by the sponsoring agency, the Division of Cancer Prevention, National Cancer Institute.Informed consent was obtained to collect tissue from study participants meeting clinical diagnostic criteria for familial BCNS and normal individuals in accordance to the ethical standards and principles expressed in the Declaration of Helsinki.

Study participants
Four female and five male BCNS cases were ascertained through the clinical practices and previous research programs (AEB and DJL).Eight age-and sexmatched, unaffected controls included relatives of cases and normal volunteers at Yale University.The mean ages of cases and controls were virtually identical -42.3+/-14.0 and 41.6+/-12.5 -although female cases were older on average than female controls while male cases were younger on average than male controls (Table 1).
Cases were chosen on the basis of meeting clinical criteria for BCNS [4].In addition, PTCH gene sequencing was performed by the CLIA-certified DNA Diagnostics Laboratory in the Department of Genetics at Yale.Four affected subjects had a mutation predicted to truncate the PTCH protein.Three had missense alterations in codons conserved from humans to invertebrates and not seen in approximately 125,000 normal controls (gnomAD browser: http://gnomad.broadinstitute.org/gene/ENSG00000185920).One of these mutations was reported previously in BCNS [40].Subject number 18 had no clearly deleterious mutation.However, the diagnostics methods used in this study would not have detected exonsize or larger deletions and duplications or other largescale rearrangements, such as inversions.This subject was heterozygous for SNPs throughout the gene, ruling out "hemizygosity" (heterozygous deletion) of the whole gene.Subject 22 had an alteration in exon 1E predicted to create a highly efficient, abnormal splice acceptor site (http://www.fruitfly.org/seq_tools/splice.html).This variant was not seen in 125,000 normal controls.Among the controls, four were relatives of cases involved in this study, and all tested negative for the mutation found in their affected relative.Among the remaining four controls, none had any of the exclusion criteria listed below: • One or more basal cell carcinomas

Specimen acquisition and tissue culture
Four 5 mm punch biopsies were obtained from unaffected skin of the upper inner arm of each case and control.Keratinocyte and fibroblast cultures were established as previously described [41].Briefly, keratinocytes were harvested from epidermis by overnight exposure to Dispase at 4°C and seeded onto mitomycin C-treated feeders of 3T3 cells.Once colonies of keratinocytes were growing, they were expanded in low calcium MCDB medium.Fibroblasts were harvested from minced pieces of dermis obtained after dispase treatment of skin biopsies and expanded in DMEM plus 10% serum.Twenty-four vials of low-passage cells with approximately 5×10 5 cells/vial were frozen in liquid nitrogen from each subject.
For experiments examining the effects of rapamycin, primary keratinocyte and fibroblast cultures were plated at approximately 25% confluence in 60 mm plates.One day after seeding, they were treated with 0μM (vehicle only), 10μM or 50μM rapamycin in DMSO (Sigma, Inc.) for 48 hours.Cells were nearly confluent when harvested.

RNA extraction
Total RNA was extracted from the keratinocyte and fibroblast cultures using the standard TRIzol ® (Invitrogen) protocol and was further purified using the RNeasy cleanup procedure (Qiagen) according to manufacturer's instructions.The quality of total RNA was assessed by Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA) for visual absence of genomic DNA contamination and integrity of 28S and 18S bands.Only samples with an A260/A280 ratio of at least 1.9 were used for microarray analysis.Samples were stored in 5μg aliquots at -70° C until use.

Microarray processing and data analysis
Gene expression profiling was performed on total 102 samples on Affymetrix Human Genome U133 Plus 2.0 GeneChip Array ® , comprising of fibroblasts and keratinocytes treated with vehicle, low-and high-dose rapamycin.Microarray processing was performed at Yale.The cDNA and cRNA preparation, cRNA labeling and hybridization to Affymetrix Human Genome U133 Plus 2.0 GeneChip ® microarrays (Affymetrix, Santa Clara, CA) was performed as per standard protocols according to manufacturer's instructions.The arrays were scanned using Hewlett-Packard GeneArray Scanner and the scanned output files were visually inspected for hybridization artifacts.The data quality was checked for the presence of spiked control cRNAs, background values caused by array autofluorescence and non-specific binding, and Q value.The microarray gene expression data have been deposited in Gene Expression Omnibus (GSE#120242).
The microarray data in.celfiles was imported and RMA normalization performed using Partek ® Genomics Suite ® software, version 6.6 © ; 2014 (Partek Inc., St. Louis, MO, USA) followed by QC metrics and preliminary exploratory analysis.Attributes were assigned to the microarray files to the random effects namely subject and scan dates, and the fixed effects namely the cell type, gender, PTCH mutation type and rapamycin treatment followed by exploratory principal components analysis to observed clustering based on various sample attributes.
For unsupervised analysis, first probe level data was collapsed to gene level data by retaining for each gene, amongst multiple possible probes, the probe that had maximum coefficient of variation.All other probes were discarded.After collapsing the probe level data to gene level data, a secondary filter to exclude genes with coefficient of variation less than 0.1 (exclude genes with CV<0.1) was applied to normalized data resulting in 5384 and 3920 genes in keratinocytes and fibroblasts respectively that were clustered using Euclidean and average linkage method as similarity measure to identify patterns of clustering in the samples.Batch effect arising due to samples being run at different time points was studied using the statistical batch effect removal tool as per the manual (Partek GS's User Manual).To ensure that the batch effect is not the major source of variation, scan dates were included in the ANOVA model as well as sources of variation were studied.To identify differentially expressed genes, a mixed-model ANOVA with restricted maximum likelihood (REML) method to estimate variance components was used in this unbalanced, mixed model study that takes into account both 'fixed' experimental factors such as with disease status/mutation type, cell type, treatment and 'random effects' such as scan dates and subject.To generate gene lists for contrasts included in the ANOVA-REML design, significance of change cut-offs of FDR <0.05 and size of fold change >1.3 or <-1.3 were used.

Gene Set Enrichment Analysis (GSEA)
GSEA was performed on the comprehensive microarray datasets for control and treatment comparisons to determine differences and enriched gene sets in the normal and BCNS groups.GSEA was performed as per software instruction to estimate the association between predefined gene sets in the MSigDb and the phenotype defined by the gene expression profiles of BCNS and normal keratinocytes and fibroblasts before and after treatment with rapamycin [64,65].The comparison of the groups was performed using GSEA by permuting the phenotype for 1000 times, and selecting the weighted scoring scheme with signal to noise statistical metric to rank genes and complete the GSEA analysis.

Ingenuity Pathway Analysis (IPA ® )
Following ANOVA, gene lists were created by applying false discovery rate (FDR), calculated using Benjamini's method [66] cut-off of <0.05 before uploading to IPA.For pathway analysis, these datasets with p-values and fold changes for each logical comparison were analyzed through QIAGEN's Ingenuity ® Pathway Analysis ® (IPA ® , QIAGEN Redwood City, https://www.qiagen.com/ingenuity)using flexible file format, Affymetrix as identifier, and Human Genome U133 Plus 2.0 array as the platform used.

Figure 1 :
Figure 1: Unsupervised hierarchical clustering of microarray data of keratinocytes (A) and fibroblasts (B) for each subject with or without (control) rapamycin treatment.Coefficient of variation of more than 0.1.The subject numbers match the participant ID numbers in Table 1.S10 indicates sirolimus/rapamycin low dose and S50 indicates sirolimus/rapamycin high dose.

Figure 2 :
Figure 2: Gene set enrichment analyses (GSEA) on BCNS keratinocytes and fibroblasts compared to normal samples.

Figure 3 :
Figure 3: Gene set enrichment analyses (GSEA) on rapamycin-treated BCNS keratinocytes and fibroblasts compared to normal samples.(A) GSEA shows that the 'GCNP_SHH_UP_EARLY.V1_UP' gene set consisting of genes up-regulated in granule cell neuron precursors (GCNP) after Shh stimulation for 3h are positively correlated and enriched in BCNS keratinocytes treated with rapamycin.Interestingly, GSEA also indicates the enrichment of a YAP conserved signature in the rapamycin treated BCNS keratinocytes.(B) For the fibroblasts, GSEA indicates that normal fibroblasts treated with rapamycin correlate with the hallmark gene sets namely 'GLYCOLYSIS', 'PI3K/AKT/MTOR SIGNALING', and 'MTORC1 SIGNALING' as compared to rapamycin treated BCNS fibroblasts.For the enrichment plots, profile of running ES score and position of GeneSet members on the rank ordered list for are denoted for oncogenic signatures 'GCNP_SHH_UP_LATE.V1_UP' (A), CORDENONSI_YAP_CONSERVED_SIGNATURE' (A), 'GCNP_SHH_UP_EARLY.V1_UP' (A), and the hallmark gene sets namely 'GLYCOLYSIS' (B), PI3K/AKT/MTOR SIGNALING (B), MTORC1 SIGNALING (B)and the respective normalized enrichment scores (NES), nominal p-values (NOM p-val), and q-values (FDR q-val) are summarized in table below the plots.

Figure 4 :Figure 5 :
Figure 4: Venn diagrams of genes commonly altered in both normal and BCNS group, with or without rapamycin treatment.Normal treated (Rx) samples were compared to normal control samples, and these samples were then compared BCNS treated (Rx) samples to BCNS control samples.(A) Keratinocyte samples.(B) Fibroblast samples.

Figure 6 :
Figure 6: Networks enriched in normal cells after rapamycin treatment.(A) Keratinocytes and (B) fibroblasts.Gene expression variation by at least 2-fold is depicted by color (red, up-regulated; green, down-regulated; gray, no significant change).

Figure 7 :
Figure 7: Networks enriched in BCNS cells after rapamycin treatment.(A) Keratinocytes and (B) fibroblasts.Gene expression variation by at least 2-fold is depicted by color (red, up-regulated; green, down-regulated; gray, no significant change).