Circadian transcriptome analysis in human fibroblasts from Hunter syndrome and impact of iduronate-2-sulfatase treatment

Background Hunter syndrome (HS) is a lysosomal storage disease caused by iduronate-2-sulfatase (IDS) deficiency and loss of ability to break down and recycle the glycosaminoglycans, heparan and dermatan sulfate, leading to impairment of cellular processes and cell death. Cell activities and functioning of intracellular organelles are controlled by the clock genes (CGs), driving the rhythmic expression of clock controlled genes (CCGs). We aimed to evaluate the expression of CGs and downstream CCGs in HS, before and after enzyme replacement treatment with IDS. Methods The expression levels of CGs and CCGs were evaluated by a whole transcriptome analysis through Next Generation Sequencing in normal primary human fibroblasts and fibroblasts of patients affected by HS before and 24 h/144 h after IDS treatment. The time related expression of CGs after synchronization by serum shock was also evaluated by qRT-PCR before and after 24 hours of IDS treatment. Results In HS fibroblasts we found altered expression of several CGs and CCGs, with dynamic changes 24 h and 144 h after IDS treatment. A semantic hypergraph-based analysis highlighted five gene clusters significantly associated to important biological processes or pathways, and five genes, AHR, HIF1A, CRY1, ITGA5 and EIF2B3, proven to be central players in these pathways. After synchronization by serum shock and 24 h treatment with IDS the expression of ARNTL2 at 10 h (p = 0.036), PER1 at 4 h (p = 0.019), PER2 at 10 h (p = 0.041) and 16 h (p = 0.043) changed in HS fibroblasts. Conclusion CG and CCG expression is altered in HS fibroblasts and IDS treatment determines dynamic modifications, suggesting a direct involvement of the CG machinery in the physiopathology of cellular derangements that characterize HS.

with idursulfase, a recombinant form of human IDS, represents the only affordable therapeutic approach presently.
The cellular processes in every life form from bacteria to humans show diurnal variations driven by an internal timing system, and the oscillation frequency has a near-24-hour period, so that the rhythmicity is defined circadian (from the latin circa, about, and dies, a day) [3][4][5]. The mammalian circadian timing system is composed by a central pacemaker and master oscillator in the suprachiasmatic nuclei (SCN) of the brain and self-sustained oscillators in the peripheral tissues [6,7]. The SCN are entrained to the environmental light/dark cycle by cues perceived through photon capture, signalling by melanopsin-containing retinal ganglion cells, and inputs conveyed by the retino-hypothalamic tract [8]. They synchronize peripheral oscillators and coordinate bodily functions by means of output pathways that may be neural projections (sympathetic and parasympathetic nerve fibers) or diffusible factors (for example hormones as cortisol and melatonin) [9,10].
The molecular mechanisms underlying the biological clock functioning consist of a transcriptional-translational feedback loop operated by a set of genes, called core clock genes, coding for proteins that in turn block gene expression, completing a cycle in approximately 24 hours [11]. The positive limb of the loop is represented by the CLOCK (or its paralog NPAS2) and ARNTL/BMAL1 (or its paralog ARNTL2/BMAL2) genes that code for the transcription factors CLOCK/NPAS2 and ARNTL/2. A heterodimer of the CLOCK/NPAS2 and ARNTL/2 proteins rhythmically activates the transcription of clock genes PER1, PER2, PER3, CRY1, CRY2, coding for PER and CRY proteins that represent the negative limb of the loop. PER and CRY proteins along with Casein Kinase Iδ/ε (CKIδ/ε) and Casein Kinase II (CKII A1, A2, B), responsible of multilevel posttranslational regulation of various clock components, cooperate to form a repression complex that translocates back into the nucleus, interacts directly with CLOCK/NPAS2 and ARNTL/2, and impede their transcriptional activity [12]. The clock gene machinery is connected to a supporting feedback loop operated by the nuclear receptors REV-ERBα/β (encoded by the genes NR1D1 and NR1D2) and retinoic acid-related orphan receptor (ROR) α/γ, where RORs positively regulate ARNTL/2 and REV-ERB expression and, in turn, REV-ERBs antagonize RORs [13,14]. Besides, SIRT1, a NAD + −dependent protein and histone deacetylase, is required for high-magnitude circadian transcription of several core clock genes, including ARNTL, PER2, and CRY1. SIRT1 counterbalances the histone acetyltransferase activity of CLOCK, binds CLOCK-ARNTL heterodimer in a circadian manner and promotes the deacetylation and degradation of PER2 [15].
TIMELESS is a core circadian clock gene in Drosophila melanogaster and is maintained in mammals, but its role in mammalian circadian clock function is not clear. TIMELESS and its partner TIMELESS interacting protein (TIPIN) interact with components of the DNA replication system to regulate DNA replication processes under both normal and stress conditions and are essential for ataxia telangiectasia and Rad3-related (ATR)-checkpoint kinase (Chk)1 and ataxia telangiectasia mutated (ATM)-checkpoint kinase (Chk)2-mediated signaling and S-phase arrest [16][17][18][19].
The clock gene oscillation drives the rhythmic expression of other genes defined clock controlled genes, for example DBP, TEF, HLF, HSF1, NFIL3 (also called E4BP4), DEC1-2 (also called BHLHE40-41), which steer downstream tissue specific genes regulating key cellular functions, such as cell cycle progression, proliferation, DNA damage response, autophagy, apoptosis, metabolism, redox equilibrium [20][21][22][23], and are involved in physiological processes such as inflammation/immune response [24,25], and paraphysiological phenomena, such as aging [26]. About 5% to 15% of genome-wide mRNA expression exhibits a circadian rhythm that is driven by the clock genes and there is a tissue specificity of cycling genes and a tissue specificity of clock controlled gene expression timing and level (output genes i.e., genes involved in the functional output of an organ) [27][28][29][30]. Accordingly, multicellular organisms are characterized by tissue-specific circadian regulation of transcription, with peripheral oscillators controlling the definite biochemical cascades relevant to their tissue or organ function and generating rhythms in various pathways, including those involved in intracellular activities and depending on the correct functioning of intracellular organelles [27][28][29][30].
The importance of the role played by the clock gene machinery in the regulation of the processes in the intracellular organelles is confirmed by the physiopathology of Niemann-Pick types A and B disease, related to deficit of sphingomyelin synthase 2 (SGMS2), the enzyme that catalyzes the transfer of phosphocholine from phosphatidylcholine onto ceramide to produce sphingomyelin, a major component of cell and Golgi membranes. Furthermore, Niemann-Pick type C disease is related to deficit of Niemann-Pick C1 (NPC1), which is crucial for the intracellular trafficking of cholesterol from the late endosome to the trans-Golgi network [31]. SGMS2 and NPC1protein levels oscillate with circadian rhythmicity driven by the clock controlled genes SGMS2 and NPC1 respectively [27,28]. Patients carrying a mutation in these genes develop a condition characterized by accumulation of sphingomyelin in spleen, liver, lungs, bone marrow and brain, causing irreversible neurological damage, and high cholesterol levels in the endosomal-lysosomal system, respectively [32].
The aim of our study was to assess the expression of clock genes and clock controlled genes in Mucopolysaccharidosis type II, and to evaluate the circadian pattern of variation and the effects of the treatment with idursulfase on the expression of clock genes and clock controlled genes at different time points.
We addressed these issues taking advantage of an in vitro model of HS, represented by human fibroblasts with the mutational features of this mucopolysaccharidosis, through evaluation of the mRNA expression levels of the core clock genes and of a panel of clock controlled genes selected by means of literature-mining [27][28][29][30]33] (for a complete list of circadian transcripts refer to: CircaDB at http://circadb. org). These evaluations are part of a whole transcriptome analysis performed by Next Generation Sequencing (NGS) in fibroblasts from healthy subjects and fibroblasts from HS patients before, 24 hours and 144 hours after idursulfase treatment. We also evaluated by qRT-PCR the expression levels of core clock genes upon serum-shock induced synchronization in normal human fibroblasts and fibroblasts of patients affected by HS before and after 24 hours of treatment with idursulfase.

Cells
Human fibroblasts from skin biopsy of five HS paediatric patients carrying different mutation in IDS gene were obtained from "Cell Line and DNA Bank from Patients Affected by Genetic Diseases", Gaslini Institute (Genoa, Italy). As healthy controls human fibroblasts from four children's circumcision were used; they were obtained from the Histology Unit of the Department of Histology, Microbiology and Medical Biotechnology (University of Padua, Padua, Italy). Written informed consents were obtained from patients at the time of biopsy and the study was approved by the Ethics Committee of the University of Padua, Padua, Italy. All cells were anonymously obtained.

Cell culture and treatment with idursulfase of HS fibroblasts
Primary fibroblasts were cultured at 37°C in 5% CO2 atmosphere in RPMI medium supplemented with 15% fetal bovine serum (FBS), 100 U/ml penicillin and 100 ng/ml streptomycin (Invitrogen Life Technologies, Milan, Italy). Fibroblasts of patients affected by HS were treated with idursulfase (Elaprase W , Shire Human Genetic Therapies, Inc, Lexington, MA, USA) at a concentration of 62.5 nM for 24 hours and the cells were harvested 24 h and 144 h after idursulfase treatment.

Whole transcriptome analysis performed by NGS technology
Total RNA was extracted using the TRIzol W Reagent (Sigma-Aldrich W ) according to the manufacturer's protocols. Isolated RNA was quantified by NanoDrop ND-1000 Spectrophotometer (Thermo Scientific, Barnstead, NH, USA) and further assessed for quality using the Agilent 2100 Bioanalyzer (Santa Clara, CA) prior to library construction. Total RNAs extracted from different cell lines were equally pooled into two specimens (Hunter and healthy) using 12.5 μg of RNA from each cell line. From total RNA of each pool mRNA purification was performed trough poliA(+) enrichment by Dynabeads W mRNA Purification Kit (Life Technologies™ Carlsbad, CA, USA). Total RNA was used for standard fragment-library preparation using the SOLiD Total RNA-Seq Kit (Life Technologies). Finally, emulsion PCR reactions were carried out by mixing appropriated amount of libraries with SOLiD beads. After PCR amplification, emulsions were broken using butanol, and the beads were washed, enriched, and terminal transferased before quantification and deposition onto a slide for sequencing. Templated beads were deposited onto one full slide, one sample for quarter (quad). Sequencing was carried out to 35 bases using SOLiD™ 3 (Sequencing by Oligo Ligation and Detection) System and following the manufacturer's instructions.

Alignment
Reads obtained from sequencing were aligned on human genome by CRIBI Biotechnology Centre (University of Padova) using the PASS software (http://pass.cribi.unipd.it). The February 2009 human genome assembly (GRCh37), publicly available at UCSC Genome Bioinformatics Site (http://genome.ucsc.edu/) was taken as reference sequence and the options best-hit, gap = 0 and maximun mismach = 3 have been selected. Only unique reads (aligning on only one gene) and only genes with coverage higher than 50% of gene length were considered.

Serum-shock induced synchronization and treatment with idursulfase of HS fibroblasts
The serum shock induced synchronization was performed as follows: approximately 4×10 5 cells/6 wells were plated the day before the experiments. Fibroblasts of patients affected by HS were treated with idursulfase at a concentration of 62.5 nM for 24 hours (Elaprase, Shire Human Genetic Therapies, Inc, Lexington, MA, USA). The day of the experiments, culture medium was exchanged with serum-rich medium (RPMI containing 50% FBS) and after 2 hours this medium was replaced with serum free RPMI [34]. The cells were harvested 1 h, 4 h, 10 h, 16 h, 22 h and 28 h after serum shock.
Quantitative real time polymerase chain reaction (qRT-PCR) Total RNA was extracted from normal primary human fibroblasts and fibroblasts of patients affected by HS before and after 24 hours of idursulfase treatment at the indicated time points using the RNeasy W Mini Kit (Qiagen S.p.a. Milan, Italy) and subsequently digested by DNase I. Quantitative Real Time PCR was performed starting from 100 ng of purified RNA using the one step quantifast SYBR Green RT PCR KIT (Qiagen). For real-time PCR, we used the following SYBR Green QuantiTect Primers purchased from Qiagen: ARNTL (QT00068250), ARNTL2 (QT00011844), CLOCK (QT00054481), CRY1 (QT00025067) PER1 (QT00069265), PER2 (QT00011207), PER3 (QT00097713). Expression levels of target gene were normalized using the housekeeping control gene TATA binding protein (TBP, QT00000721). Values of mRNA expression levels of clock genes were calculated using the formula 2 −ΔΔCt .

Statistical analysis
Results are expressed as means ± SE of at least three different experiments. Comparisons were made using Student's t-test as appropriate. The limit of statistical significance was set at to p < 0.001 for comparisons of levels of mRNA expression determined by NGS, and p < 0.05 for comparisons of levels of mRNA expression determined by qRT-PCR. Each time series of clock gene expression levels was analyzed for circadian rhythm characteristics by the single cosinor procedure involving the fit of a 24-h cosine curve to the data by least squares linear regression, in order to accurately describe waveforms and rhythm characteristics. An R2 value and a p value for the rejection of the zeroamplitude assumption were determined for each component in the cosine model separately and overall, with rhythm detection considered statistically significant if p ≤ 0.05 and borderline significant if p > 0.05 and p < 0.10 for any period tested. Circadian (24-h) characteristics were summarized from the 24-h cosine. Rhythm characteristics determined from the best-fitting cosine model include the MESOR (the middle of the cosine representing an adjusted average if unequal sampling), the amplitude (half the distance from the peak and trough of the best fitting curve), and the acrophase of the cosine model (the peak of a single component cosine). All analyses were performed using the MATLAB 6.5 statistical package (MathWorks, Natick, MA, USA).

Evaluation by qRT-PCR after Serum-Shock Induced Synchronization of changes of clock gene expression levels in healthy fibroblasts, untreated HS fibroblasts and HS fibroblasts after 24 hours of treatment with idursulfase
After serum shock synchronization and qRT-PCR analysis, statistically significant differences were evidenced in the mRNA expression levels of ARNTL2 at 10 h (p = 0.036), PER1 at 4 h (p = 0.019), PER2 at 10 h (p = 0.041) and 16 h (p = 0.043), between idursulfase treated HS fibroblasts and untreated HS fibroblasts ( Figure 1). Fitting of cosine curves with a 24-hour period to raw data and plotting as polarograms clock gene expression levels in control fibroblasts and fibroblasts of patients affected by HS before and after 24 hours of idursulfase treatment (Figures 2, 3 and 4) evidenced statistically significant and borderline significant rhythms of clock gene expression, and advance of phase os- Figure 1 x-y plots representing the time related profiles of expression of ARNTL, ARNTL2, CLOCK, CRY1 PER1, PER2, and PER3 after synchronization with serum shock in control fibroblasts and fibroblasts of patients affected by Hunter Syndrome before and after treatment with iduronate-2-sulfatase. Original units standardized to T0 (1 h after serum shock) and combined for analyses. * p < 0.05 for ANOVA followed by all pairwise multiple comparison procedures with Student-Newman-Keuls Method. cillation after idursulfase treatment. MESOR, amplitude, acrophase and p values are reported in Table 1.

Semantic hypergraph-based analysis of circadian gene expression
The aim of this analysis was to determine the functional groups of genes that play any key-role in the development of HS. The analysis workflow is here preliminary sketched: (i) RNA-Sequencing: Detection of core clock gene and clock controlled gene expression profiles; (ii) Descriptive analysis and feature selection: Selection of the most varying genes by the ensemble of principal component analysis and linear discriminant analysis; (iii) Qualitative hypergraph building: Inference of weighted and signed hypergraphs; (iv) Quantitative undirected graph construction: Mathematical agglomeration of weights over the edges; (v) Weighted topology indices calculation: Computation of several weighted topology indices for each node and normalization of individual fold changes by the most discriminant index; (vi) Determination of the strongest connected components: Hierarchical clustering and matching with critical biological functions [35,36].

Descriptive analysis and features selection
We have preliminarily assessed the expression variation of clock genes and clock-controlled genes, through five different states: HS (H) vs. Control (C), Treatment at time 1 (T1) vs. C, Treatment at time 2 (T2) vs. C, and T1 vs. H, T2 vs. H. We discarded genes that did not exhibit any relevant expression change between these states. We combine a principal component analysis to determine those genes that maximally vary through the first three components, and a linear discriminant analysis to project the genes and their expressions onto the most varying components' axes ( Figures 5 and 6).

Qualitative hypergraphs building
We have assembled a hypergraph connecting genes by querying and merging a number of heterogeneous data sources: (i) Interpro and PFAM (protein domains), (ii) Gene Expression Omnibus (Co-localization and co-expression), (iii) BIOGRID and IREF (genetic interactions), (iv) PathwayCommons, IMID, NCI_NATURE, REACTOME, KEGG and BIOCARTA (pathways), (v) BIOGRID, BIND, HPRD, INTACT, MINT, MPPI and OPHID (physical interactions) and (VI) curated literature (predicted interactions). Two genes were connected by an edge, whenever at least an interaction evidence of any of the above mentioned interaction categories was found. Several pair of genes resulted to be connected by more than one edge. We built a weighted and signed hypergraph with weights over the edges (carrying the reliability of the corresponding interactions) and weights over the nodes (carrying the fold change expression of the corresponding genes) [37][38][39].

Quantitative undirected graph
We deterministically inferred an undirectedm graph from our hypergraph, by applying an injective function to the sets of weights over the edges connecting any two pairs of genes. This simple agglomerative weighting function takes the weights of the edges connecting any two nodes in input, and gives a unique value in output. Constitutively, it gives more and more importance to the pairs that are connected by multiple edges.
where n is the number of edges connecting any two nodes A and B, and i refers to the ith edge. Thus, we obtain a new graph that contains the same set of nodes, and individual edges connecting any two nodes, with the newly calculated weights over them (Figure 7).

Weighted topology indices calculation
We calculated the closeness topology index for each gene. We recall that the un-weighted formulation of closeness between any two nodes of a graph measures the length of their shortest paths. The concept of closeness originates from the inverse definition of farness, namely as the inverse of the sum of the distance from a node to all other nodes. Thus, the more important a node is, the lower is its total distance to all other nodes. Closeness can be regarded as a measure of how fast it will take to spread information from a node to all other nodes sequentially. Contrarily, the weighted formulation of closeness takes into account paths whose lengths are not the sum of the minimum number of hops from a node to all other nodes, but the sum of the inverse weights over the edges of the shortest paths from a node to all other nodes. This formulation gives importance to the reliability of a connection, when determining the essentiality of each node, thereby discriminating between reliable and not-reliable interaction paths.
We have subjected the top 26 closest genes to other well-known topological indices (i.e., degree, betweenness, clustering coefficient and topological importance) and found that 5 (i.e., AHR, HIF1A, CRY1, ITGA5 and EIF2B3, see the plot below) were highly ranked by 4 out of 5 indices (Table 2 and Figure 8).

Strongly connected components determination
After having scaled the individual fold changes with the closeness values of each gene, we have finally determined the collaborative actions of such genes by extracting the most cohesive clusters out of their interaction network. We performed a hierarchical clustering on the expression profiles of clock and clock controlled genes, by considering the Squared Euclidean distance as metric. We have subjected the found clusters to functional enrichment analysis against the GO FAT category, and determined five clusters significantly associated to at least a biological process or pathway ( Figure 9).

Discussion
The molecular oscillator ticking in every cell regulates the customary sequence of intracellular activities, allowing the coordination of key pathways and the compartmentalization in the temporal dimension of poorly-compatible biochemical processes [40]. Altered functioning of the clock gene machinery may determine loss of time-of-day specific transcription of clock genes and clock controlled genes, causing severe deregulation of cellular homeostasis, cell dysfunction and biochemical and structural derangements that may lead to cell death and tissue dysfunction [41].
The aim of our study was to evaluate the functioning of the clock gene machinery in HS, and to address this issue we analyzed and compared mRNA expression levels of the core clock genes CLOCK, NPAS2, ARNTL1, ARNTL2, PER1, PER2, PER3, CRY1, CRY2, CSNK1D, CSNK1E, CSNK2A1, CSNK2A2,CSNK2B, NR1D1, NR1D2, RORA, SIRT1, TIMELESS, TIPIN, and a set of clock controlled genes in normal human fibroblasts and fibroblasts derived from subjects affected by Mucopolysaccharidosis type II. We also evaluated in vitro the effects of the treatment with idursulfase on circadian gene expression at different time points.
Analyzing data obtained from NGS we observed altered mRNA levels of some of the clock genes and clock controlled genes and a variable response to idursulfase treatment.
In HS fibroblasts CLOCK and NPAS2 showed lower expression levels. CLOCK and its paralog NPAS2 encode  transcription factors belonging to the basic helix-loophelix/Per-Arnt-Simpleminded (bHLH-PAS) domain family, and represent key elements in the positive limb of the transcriptional-translational feed-back loop that drives circadian rhythmicity of the oscillatory functions underlying cell homeostasis. In particular, their encoded proteins have histone acetyltransferase activity and are involved in chromatin remodeling with specificity for histones H3/H4 [42]. ARNTL, the heterodimerization partner of CLOCK and NPAS2, enhances their enzymatic function, and the decreased expression of these central components of the molecular clock could hinder the correct clock gene expression in the fibroblasts of Mucopolysaccharidosis type II patients.
A lower expression level was found in HS fibroblasts for CRY1, whereas PER1 was up-regulated. The mammalian Period and Cryptochrome genes are both E-box-regulated genes, but in peripheral tissues CRY1 mRNA expression peak is delayed by several hours with respect to that of PER1. CRY1 usually shows evening-time expression and serving as a strong repressor at morningtime elements ensures a delay in feedback repression in the molecular clock. This phase delay in CRY1 transcription is required for customary mammalian clock function [43]. The E-box (CACGTG) is crucial for daytime transcriptional activity and the delay originates from interactions between the proximal E-box and ROR response elements (RORE) present in the CRY1 promoter [44].
A higher expression level was found in HS fibroblasts for CSNK1D and CSNK1E, encoding the serine/threonine protein kinase casein kinase (CK) I δ and ε respectively, which are key regulators of metazoan circadian rhythmicity. CKI binds to and phosphorylates the Period proteins, wich in turn interact with a variety of circadian regulators, suggesting the possibility that CKI may interact with and phosphorylate additional clock components as well [45]. Cryptochrome proteins are phosphorylated by CKI only when both proteins are bound to mammalian Period proteins [12]. ARNTL is also a substrate for CKI in vitro, and CKI activity positively regulates ARNTL-dependent transcription from circadian promoters in reporter assays [12]. CKI phosphorylates multiple circadian substrates and may exert its effects on circadian rhythm in part by a direct effect on ARNTLdependent transcription [12,45]. The up-regulation of these protein kinases found in HS fibroblasts could play a crucial role in the derangement of the clock gene machinery and the downstream controlled pathways, but on the other hand CKs could represent important molecular targets for new therapeutic approaches.
In HS fibroblasts we found higher expression of NR1D1, and lower expression of NR1D2. NR1D1 and NR1D2, activated by CLOCK/NPAS2:ARNTL/ARNTL2 heterodimers, code respectively for the orfan nuclear receptors REV ERBα and REV ERBβ, which realize a supplementary loop amplifying and stabilizing the oscillation of the molecular clockwork. REV-ERBα and REV-ERBβ act as negative transcriptional regulators by binding RORE in gene promoters, preventing the binding of the positive transcriptional regulator RORα, and negatively regulating the expression of ARNTL, CLOCK, and CRY1 [46,47]. In addition, REV-ERBα is directly involved in lipid metabolism antagonizing the opposite role of RORα and inhibiting the expression of genes coding for apolipoprotein C-III, a constituent of very-low-density lipoproteins [48].
A lower expression level was found in HS fibroblasts for SIRT1, which regulates metabolic and stress responses acting in concert with the circadian rhythm machinery [49]. SIRT1 is involved in a number of cellular processes, including gene silencing at telomere and mating loci, and modulates cell survival by inhibiting apoptosis or cellular senescence induced by external challenges, including DNA damage and oxidative stress [50]. SIRT1 impinges on Ku70, peroxisome proliferator activated receptors, p53 and the forkhead box O (FOXO) family of transcription factors, and is modulated by active regulator of SIRT1 (AROS), hypermethylated in cancer 1 (HIC-1), deleted in breast cancer 1 (DBC1) and E2F transcription factor 1 (E2F1) [49][50][51]. Interestingly, E2F1 was greatly over-expressed in HS fibroblasts, but normalized after 144 hours of IDS treatment, whereas FOXO1 was severely down regulated, and responded faintly to IDS treatment. Besides, in HS fibroblasts there was severe down-regulation with no response to IDS treatment of NAMPT, encoding the enzyme controlling the synthesis of nicotinamide adenine dinucleotide (NAD + ), a cofactor of SIRT1, linking nutrient sensing and circadian regulation [52].
After 24 hours of treatment with idursulfase in fibroblasts of patients with Mucopolysaccharidosis type II the    the expression of a huge number of output genes involved in the control of important cellular and tissue processes. An example is represented by SERPINE1, which showed higher expression levels in HS fibroblasts when compared to healthy fibroblasts, increased after 24 hours of idursulfase treatment, and decreased after 144 hours. SERPINE1 codes for the plasminogen activator inhibitor (PAI)-1, the major component of inhibitors of fibrinolysis, whose activity shows a clear circadian oscillation peaking in the early morning. CLOCK:ARNTL and CLOCK:ARNTL2 heterodimers activate the SERPINE1 promoter, driving the circadian variation in circulating PAI-1, and CLOCK: ARNTL2 heterodimer, binding two E-box enhancers in the promoter, shows double capability to activate PAI-1 expression [53].
The timing of processes that underlay cell function and the customary sequence of activation of key pathways depend on the correct ticking of the biological clock. Accordingly, the altered expression of clock controlled genes found in HS fibroblasts could hinder normal output of the clock gene machinery (BHLHE40, BHLHE41, NFIL3), and cellular activities, such as DNA transcription (CTNNB1,  FOS, FOXL2, FOXO1, HIF1A, HSF1, ID2, JUN, KEAP1 The semantic hypergraph-based analysis of circadian gene expression brought out five gene clusters significantly associated to at least one biological process or pathway (cell cycle, DNA damage response, inflammation, liver development, cell motility, glucose homeostasis, response to cold and to decreased oxygen levels, the P53 signaling pathway, the positive regulation of epithelial cell proliferation, tissue morphogenesis, response to lipid, and the adipogenesis pathway). Besides, five genes, AHR, HIF1A, CRY1, ITGA5, and EIF2B3, were highlighted as top ranked by all the five considered topological indices and characterized by multifaceted centrality nature in the network analysis of the circadian transcriptome. The results suggest an interesting connection between circadian and hypoxic pathways and the toxicological signal mediator, the Aryl Hydrocarbon Receptor (AHR). Upon xenobiotic binding, AHR transcriptionally activates xenobiotic metabolizing enzymes and regulates molecules involved in the signaling of nuclear factor-erythroid 2-related factor-2 (Nrf2), p53 (TRP53), retinoblastoma (RB1), and NF-κB, influencing cellular responses to oxidative stress and inflammation. Through activation of these signaling pathways, regulation of cell cycle, and interactions with hypoxia-inducible factors (HIFs), AHR takes part in endogenous developmental functions during hematopoietic stem-cell maintenance and differentiation [55]. Intriguingly, the deficit of CRY proteins causes constitutive NF-κB and protein kinase A signaling activation, leading to unrelenting increase of proinflammatory cytokines, and ITGA5 encodes a receptor for fibronectin and fibrinogen, playing a role in chemotaxis, leukocyte migration, cell-cell adhesion, angiogenesis, and blood coagulation [56,57].
On the other hand, a cross-talk exists between circadian and AHR signaling pathways at genetic, epigenetic and proteomic level, and this interaction might play a role in the regulatory influences maneuvered by circadian rhythmicity on cell physiology [58]. Regarding to the interaction between hypoxic pathways and AHR complex, it has been proven that in endothelial cells AHR has a physiological function in the absence of exogenous ligands, and upon activation mediates toxicity of halogenated aromatic hydrocarbons. In addition, evidence suggests that hypoxia induces HIF dependent gene expression, significantly reducing AHR expression level and AHR-mediated expression of cytochrome P-450 phase I xenobiotic metabolizing enzyme CYP1A1 [59]. A cross-talk exists between hypoxic and circadian pathways operated by the PAS protein family members PER and CLOCK and HIF-1α. Interestingly, the components of AHR, circadian and hypoxic pathways are characterized by a PAS domain that serves as an interface for protein-protein interactions [60]. Moreover, cell proliferation and differentiation both in vitro and in vivo are highly influenced by oxygen concentration, and neuronal stem cell proliferation and neuronal and oligodendroglial differentiation are enhanced by a mild level of hypoxia [61]. Intriguingly, EIF2B3 encodes one of the subunits of initiation factor EIF2B, which catalyzes the exchange of eukaryotic initiation factor 2-bound guanosine-diphosphate (GDP) for guanosine-5'-triphosphate (GTP), representing an essential factor for protein synthesis, and mutations in this gene have been associated with neurodegenerative and white matter diseases [62][63][64].
Taking into account some limitations of the study, represented by the limited number of patients, suffering from a disease that is clinically heterogeneous in terms of onset, severity and progression, the dynamic variation of expression levels observed for core clock genes and clock controlled genes in the fibroblasts of patients affected by HS at different time points of treatment with idursulfase may be an indirect evidence of the key role played by the molecular clock in the regulation of the complex array of cellular functions in this mucopolysaccharidosis. A paradigm is represented by the altered expression of the clock controlled genes NPC1 and SGMS2, which are down-regulated in HS fibroblasts and increase temporarily after idursulfase treatment.

Conclusion
Biological processes show time related variations driven by genetically encoded oscillators operated by transcriptional/translational feedback loops hardwired by the clock genes and their coded circadian proteins. Data obtained in our study show that circadian gene expression is variably altered in the fibroblasts of subjects affected by HS, the treatment with idursulfase determines some modifications in the expression levels of clock genes and clock controlled genes, but these changes are temporary and fade in the course of treatment. Altered functioning of the clock gene machinery in HS may determine severe deregulation of cellular homeostasis, cell dysfunction and biochemical and structural derangements that may lead to cell death. These alterations are related to the key role played by the molecular clockwork in the control of downstream gene expression regulating a complex array of cellular functions, such as molecule biosynthesis, post-translational modification, processing, transport, conjugation, internalization and degradation, and cell processes such as cell cycle, autophagy, apoptosis and DNA damage response.The alteration of clock gene expression levels and the response to the deficient enzyme suggest a direct involvement of the molecular clock machinery in the physiopathology of cellular derangements that characterize Mucopolysaccharidosis type II, opening the way to possible new therapeutic strategies.

Additional file
Additional file 1: Table S1. Expression levels of core clock genes and clock controlled genes evaluated by whole transcriptome analysis performed through Next Generation Sequencing technology in normal human fibroblasts (Control) and fibroblasts of mucopolysaccharidosis Type II patients before (Hunter), 24 hours (T1) and 144 hours (T2) after idursulfase treatment.