Dynamic Expression of Genes Involved in Proteoglycan/Glycosaminoglycan Metabolism during Skin Development

Glycosaminoglycans are important for cell signaling and therefore for proper embryonic development and adult homeostasis. Expressions of genes involved in proteoglycan/glycosaminoglycan (GAG) metabolism and of genes coding for growth factors known to bind GAGs were analyzed during skin development by microarray analysis and real time quantitative PCR. GAG related genes were organized in six categories based on their role in GAG homeostasis, viz. (1) production of precursor molecules, (2) production of core proteins, (3) synthesis of the linkage region, (4) polymerization, (5) modification, and (6) degradation of the GAG chain. In all categories highly dynamic up- and downregulations were observed during skin development, including differential expression of GAG modifying isoenzymes, core proteins, and growth factors. In two mice models, one overexpressing heparanase and one lacking C5 epimerase, differential expression of only few genes was observed. Data show that during skin development a highly dynamic and complex expression of GAG-associated genes occurs. This likely reflects quantitative and qualitative changes in GAGs/proteoglycans, including structural fine tuning, which may be correlated with growth factor handling.


Introduction
During various cell signaling processes, glycosaminoglycans (GAGs), such as heparan sulfate (HS), chondroitin sulfate (CS), and dermatan sulfate (DS), play a role in binding, guiding, and modulating signaling molecules, e.g., growth factors and morphogens [1][2][3]. In skin this role can be illustrated by the importance of GAGs in adult wound healing [2,4] and in the extracellular matrix architecture formed during dermal development [5,6]. A further example to illustrate the importance of GAGs comes from mice overexpressing heparanase, an enzyme involved in the degradation of HS, showing accelerated hair growth [7], indicating its involvement in hair follicle morphogenesis and homeostasis. Other observations show that HS is involved in hair follicle cycling, sebaceous gland morphogenesis, and homeostasis [8]. Finally, HS and heparanase influence wound healing in adult mice by enhancing keratinocyte migration and stimulating blood vessel maturation [9]. Taken together, GAGs play an important role in skin healing and development and this prompted us to evaluate the expression of GAG related genes during (embryonic) development in skin.
Inhibition of the expression of genes coding for enzymes involved in GAG modification reactions clearly indicates the importance of GAGs during organogenesis [10], especially with respect to growth factor handling. For example, mice deficient in Ndst1 (N-deacetylase sulfotransferase isoenzyme 1) die neonatally due to several defects in which defective sonic hedgehog (Shh) signaling is implicated [11,12]; mice deficient in Hs2st (heparan sulfate 2-O sulfotransferase) or Glce (glucuronic acid epimerase) display renal agenesis [13,14], whereas mice deficient in Hs6st1 (heparan sulfate 6-O sulfotransferase isozyme 1) show aberrant signaling of VEGF (vascular endothelial growth factor) and impaired lung development [15]. A skin phenotype of the above mouse models, however, has not been reported. In general, it is thought that specific modifications of the GAG chain are involved in the binding and modulation of signaling molecules resulting in cell-type and/or tissue specific reactions [2,3]. GAG mimetics like the RGTAs (regenerating agents) have been used to treat skin disorders and improve skin healing [16,17]. To obtain insight in GAG metabolism during skin development we studied the expression of GAG related genes covering six functional classes ranging from the synthesis of precursor molecules to the synthesis and degradation of GAGs. In addition, we probed the expression of a number of (GAG binding) signaling molecules.

Materials and Methods
An overview of the experimental setup on the gene expression during murine skin development is given in Figure 1.

Animals for the Study on Skin
Development. NIH guidelines for the care and use of laboratory animals (NIH publication 85-23 Rev. 1985) were followed. The study was approved by the Ethics Committee of the Radboud university medical center (DEC 2005-111, project: 81027). C57BL6/j mice were obtained from Elevage Janvier (Le Genest Saint Isle, France). Mice aged 90 days (90 days post birth [P90]) were used for timed mating and dorsal skin was collected at 14 days (E14) and 16 days after conception (E16). At E14 hair follicle development is initiated, and at E16 this process is almost completed in combination with a stratified epidermis and organized dermis [18,19]. For the RNA samples of E14, dorsal skin of seven embryos from one female was pooled and used for RNA isolation. Skin was isolated at E14 by snap freezing the whole embryo in liquid nitrogen followed by scraping the skin layer in a cryomicrotome with a scalpel to minimize contamination with other embryonic tissues (skin is very thin at this time point). Samples were stored at -80 ∘ C. RNA samples for E16 were taken from two females, collecting dorsal skin form 7 embryos each. In addition, skin from 1-day old pups (P1) and adult mice (P90) was collected. At P1 skin is more organized and has been exposed to air [18,19]. For the two dorsal skin samples for P1, three pups from two females were taken per sample. Two adult three-month old mice were used for the two dorsal skin samples at P90. Samples for RNA isolation for E16, P1, and P90 were collected by removing dorsal skin and snap freezing it in liquid nitrogen and storage at -80 ∘ C.  [7,20]. For RNA isolation two wildtype and two mutant mice were used of both mouse models.  [18,21]. Gene level expression data were calculated for the CORE transcripts (probe sets supported by RefSeq mRNAs) using Affymetrix Expression Console software with quantile normalization (all arrays are considered to have an equal intensity distribution), GC-content background correction (probes with high GC-content hybridize better, corrected for with built-in probes with different known GC-contents) and summarization with the RMA algorithm [22]. Data were imported into GeneSpring GX 7.3 (Agilent Technologies), duplicates were averaged, and the expression of each transcript was normalized to the median per array.

RNA Isolation
Real Time-Quantitative PCR (qPCR) was performed using custom designed Taqman Low Density arrays (TLDA) (Applied Biosystems, Nieuwerkerk aan de IJssel, the Netherlands) containing probes against genes involved in GAG metabolism and GAG binding proteins (Supplementary data  Table 1). Glce and Hpse samples were analyzed using qPCR using custom designed TLDA with an adapted design containing additional GAG related genes (Supplementary data  Table 2).
For the TLDA cards, 100 ng cDNA in Taqman Universal PCR Master Mix (Applied Biosystems) was loaded on the TLDA card per slot and run on a 7900HT Fast Real Time PCR System (Applied Biosystems). Expression was analyzed based on the threshold cycle (Ct) which was obtained using the SDS 2.3 software and RQ Manager 1.2 of Applied Biosystems using the combined expression data of the tested TLDA cards. In Microsoft Excel the reference genes for ΔCt calculation were checked for stability of expression by analyzing the results of the reference genes across all used TLDA cards and selecting the reference genes with the smallest deviation across the cards tested. Subsequently ΔCt values were calculated using a reference gene with the smallest difference between the average Ct found for the gene of interest and for the reference gene. The obtained ΔCt values were further processed using the 2 −ΔΔCt method using P90 as a calibration point in case of the developmental study and the wild type background (C57BL6 mice) data in case of the two mouse models [23].

Statistical Methods.
Statistical significance of the exon array data was analyzed using ANOVA and Benjamini-Hoghberg multiple testing correction [24]. Statistical significance of the TLDA card data was tested with an unpaired T-test (2-tailed) using Microsoft Excel. Data with a statistical threshold of p<0.10 and a fold threshold of >2.0 were considered statistically significant.

Results
Genes involved in glycosaminoglycan (GAG) synthesis, modification, and degradation were studied during skin development at 14 and 16 days after conception (E14 and E16, respectively) and at one day after birth (P1) and compared to mature skin of a 3 month old mouse (P90). In addition, two mouse models, a Glce knockout mouse (E18.5) and an Hpse overexpression mouse model (P70), were analyzed. Taqman Low Density Array (TLDA) cards were designed to contain genes involved in GAG metabolism (Supplementary  data Tables 1 and 2). The expression data obtained using TLDA cards and exon arrays were screened for genes with 2fold differential expression at a statistical threshold of p<0.10 (Tables 1, 2, 3, and 4, Supplementary data Tables 4, 5, and 6).
In Tables 1-3 and Supplementary data Table 3, an overview is given of the differentially expressed genes applying TLDA cards and exon arrays. In all categories of genes involved in GAG metabolism, i.e., production of precursor molecules, core proteins, synthesis of linkage region, polymerization, modification and degradation of the GAG chain, and differences in expression were found (Tables 1-3). This indicates a highly dynamic expression pattern during skin development. Some isoenzymes were upregulated, whereas other isoenzymes were downregulated, further stressing metabolic complexity. This is, for instance, the case with GFPT1 and 2, both rates limiting enzymes involved in the production of hexosamines, and the isoenzymes HS 3-O sulfotransferase 6 and 3b1.
The upregulated expression of genes involved in the synthesis of the linkage region may signal increased GAG synthesis during development since after the formation of the linkage region the GAG chain is formed. For HS polymerization differential expression was found for, e.g., Extl1 and Extl2. Extl1 showed downregulation at E14 while Extl2 was upregulated, and both enzymes are involved in the initiation and elongation of the HS chain [2].
During and after synthesis of the glycosaminoglycan chain, disaccharide units within the chain are specifically modified. These modifications determine which effector molecules can bind to the chain and thus play a role in cell signaling [1,3]. Upregulated expression was found for three of the four N-deacetylase/N-sulfotransferases (Ndst; TLDA cards, Table 2), especially isoenzyme Ndst3. Upregulation was also found for two out of seven genes coding for 3-O-sulfotransferases (Hs3st1 and Hs3st3b1), involved in 3-O sulfation of GlcNS and GlnNAc residues, whereas one was downregulated (Hs3st6).
The GlcNS and GlcNAc residues can also be 6-O sulfated by 6-O-sulfotransferases (Hs6st) [2] and selectively desulfated extracellularly by two sulfatases (Sulf1 and Sulf2) aided by two cofactors (Sumf1 and Sumf2) [2,27]. Hs6st2 was upregulated at all time points (Table 3). Sulf1 was upregulated during embryonic development, whereas Sumf2 was upregulated at E14 (Tables 2 and 3). These results indicate that specific expression of GAG modifying enzymes may play a role in specific cellular signaling during skin development.
Within the class of genes encoding for GAG chain degradation enzymes, two genes were differentially expressed. Heparanase expression was downregulated at E14 and P1 (Tables 2 and 3), whereas N-sulfoglucosamine sulfohydrolase    Numbers in italic are significant (p<0.10); numbers in bold are >2-fold differentially expressed. Gene symbols indicated with a †-symbol are normalized using GAPDH as a reference gene. Gene symbols indicated with a ‡-symbol are normalized using TBP as a reference gene. Genes, for which a signal was not or only partly detected at a given time point or multiple time points and therefore a fold change and/or p value could not be calculated based on the available data, are given as "not detected." Gene symbols for which all time points were classified as "not detected" do not show a symbol for the used reference gene due to lack of data for a calculation.
In addition to genes involved in GAG metabolism, the TLDA card contained 37 genes encoding growth factors, which were also present in the microarray (Tables 2 and 3). Differential expression was found by both TLDA card and microarray analysis for 10, 9 and 4 growth factors at E14, E16 and P1 respectively. Examples are insulin-like growth factor 2 (Igf2), wingless-related integration site 6 (Wnt6), and Wnt7b. Igf2 was dramatically upregulated at all time points, as expected based on previous research [18]. Wnt6 was also upregulated at all time points, while Wnt7b was upregulated only during embryonic development.
Next to their expression during development, gene expression of GAG-associated genes was studied in a Glce (glucuronyl epimerase) knockout mouse model and a heparanase overexpression mouse model using TLDA cards. In the Glce knockout mice six genes were differentially expressed ( Table 4). Three of them are involved in CS and DS proteoglycans and were downregulated, i.e., aggrecan (Acan), asporin (Aspn), and chondroitin sulfate N-acetylgalactosaminyltransferase 2 (Csgalnact2). Up/downregulation was not found for HS related genes, except for Glce, which was downregulated as expected. For the heparanase overexpression mouse model, in which a human heparanase was overexpressed [7], the results showed only one gene to be differentially expressed, i.e., aggrecan (Acan) which was 2.5-fold upregulated. The complete results of both the Glce knockout mouse and the Hpse overexpression mouse are given in Supplementary data Table 6.

Discussion
GAGs play a regulating role during embryonic development of various organs [1][2][3]. Therefore, we examined the expression of genes involved in GAG metabolism during skin development using custom designed Taqman Low Density Arrays (TLDA card) and exon arrays. To structure the data we studied gene expression in six functional classes, viz. the production of precursor molecules, the synthesis    Numbers in Italic are significant (p<0.10); numbers in bold are >2-fold differentially expressed. Genes indicated as "not measured" represent genes for which probes were not available on the used exon array version. of core proteins and the linkage region, and the synthesis, modification, and degradation of the GAG chain proper. In addition we studied a number of growth factors, since GAGs are involved in their regulation including growth factor diffusion and signaling [3,28].
With respect to core proteins, the heparan sulfate proteoglycans syndecan and glypican showed notable differential expression (Tables 2 and 3). Glypicans play an important role in development and cell signaling [12,26,29], and we found upregulation of 3 out of 6 glypican core proteins. Gpc3 was upregulated during embryonic development and one day postbirth, suggesting that this glypican has a role during skin development. A possible function of Gpc3 in skin has been suggested for the Gpc3-null mouse, which showed pigmentation defects [30]. Humans deficient in Gpc3 suffer from the Simpson-Golabi-Behmel syndrome (SBGS). Based on the symptoms of SGBS and the phenotype found for the Gpc3-null mice, it has been suggested that Gpc3 is involved in the regulation of hedgehog signaling [31], a signaling pathway involved in hair follicle development [32]. Surprisingly, the Gpc3-null mice did not show a defect in appendage formation [30], indicating a functional but not essential role. Further Numbers in Italic are significant (p<0.10); numbers in bold are >2 fold differentially expressed. All genes were normalized using 18S RNA as a reference gene.
research is needed to elucidate the role of Gpc3 and the two other differentially expressed glypicans, i.e., Gpc2 and Gpc6.
Syndecans are described to take part in adult wound healing [33]. We found downregulation of the core proteins of two syndecans during embryonic development, which could indicate that these proteoglycans do not play a major general role during skin development. Specific roles, such as the involvement of Sdc1 in hair follicle development, as described on basis of immunohistochemical data [34], can, however, not be excluded.
In the class of GAG chain polymerization, we found differential expression of genes encoding for the initiation of HS or CS/DS synthesis. HS chain polymerization is initiated by the addition of GlcNAc by Extl2 [35] or Extl3 [36], while CS/DS chain polymerization is initiated by the addition of GalNAc by Csgalnact1 [37]. Extl2 was upregulated during early skin development (Table 2), while Csgalnact1 was downregulated (Supplementary data Table 5), which suggests that during early skin development HS production is stimulated in comparison to CS/DS production.
Enzyme mediated chemical modifications of the GAG chains result in the creation of specific binding sites for effector molecules [38]. Enzymes forming the class of Ndeacetylase/sulfotransferases (Ndst's) are initiating elements in this respect. Especially Ndst3 was upregulated, being one of four enzymes responsible for the removal of the acetyl group from the N-acetylated glycosamine and for the addition of a sulfate group. The additional expression of Ndst3 in combination with Ndst1 and Ndst2 points to the fine tuning of HS chains for specific recognition of ligands. Ndst3 has a higher deacetylation activity in comparison to the N-sulfotransferase activity, while Ndst1 and Ndst2 have a slightly higher N-sulfotransferase activity [39]. In addition, the data on the expression of heparan sulfate 3-O sulfotransferases (Hs3sts) [40] and heparan sulfate 6-O-sulfotransferases (Hs6sts) [41] suggest dynamic and specific modification of HS chains.
Three genes encoding for enzymes involved in HS and CS/DS degradation were differentially expressed, one of them being Hpse (heparanase). Hpse is downregulated at E14 and at P1, but not at E16 at which time point hair follicle development is taking place. Hpse has been reported to be involved in this process [42,43].
Glycosaminoglycans are involved in growth factor regulation during developmental processes [1,2]. We therefore studied 37 growth factors implied in skin development. A number of genes encoding growth factors were differentially expressed during development and the data are in line with earlier results for, e.g., Igf2 [18], Wnt6, and Wnt7b [44]. Although speculative, the dynamics in GAG structure may be correlated with the dynamics of growth factors.
Next to skin development we also studied gene expression in skin of a Glce (glucuronyl epimerase) knockout mouse and an Hpse (heparanase) overexpression mouse [7,13]. In the Glce knockout mice relatively few genes were differentially expressed, suggesting that skin is relatively unaffected by the lack of Glce in line with the observation that skin in these mice is phenotypically normal [20]. The skin phenotype of the Hpse overexpression mouse shows accelerated hair growth [7]. Gene expression analysis of this model showed only one differentially expressed gene (aggrecan). These results may touch upon the regulation of translation of mRNAs coding for GAG related enzymes. Enzymes involved in the synthesis and modification of GAGs as well genes coding for (some) growth factors share a common alternative translation mechanism via IRES sites [45,46]. In general mRNAs are translated by the ribosomal scanning mechanism which scans for short leader sequences of 50 to 70 nucleotides [46,47]. The leader sequences of the HS modifying enzymes and growth factors are characterized by long but structured sequences, which are not recognized by the ribosomal scanning mechanism [46,47]. Within these sequences internal ribosomal entry sites (IRES) allow alternative translation, e.g., under stress conditions [47]. This indicates that in addition to mRNA levels an additional control mechanism on the translational level may be present. In addition, other types of regulatory levels are known including the interaction of biosynthetic enzymes with each other and the (possible) presence of large biosynthetic complexes (GAGosomes) [48]. This makes the regulation of GAG biosynthesis very complex, gene expression being only a part of it.
Taken together, it is concluded that a highly dynamic expression of genes involved in GAG metabolism and in GAG binding growth factors is associated with skin development. This indicates the importance of fine tuning of GAG structures during developmental processes. Further studies should focus on the biochemical analysis of the GAGs chains themselves.

Data Availability
The EXON array data used to support the findings of this study are included within the article and are provided via [18]. The Taqman low density array data used to support the findings of this study are included within the article. The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.