extfieldsofstudy: Biology Medicine

text:
Lentiviral Vector Mediated Thymidine Kinase Expression in Pluripotent Stem Cells Enables Removal of Tumorigenic Cells

Embryonic stem cells (ES) and induced pluripotent stem (iPS) cells represent promising tools for cell-based therapies and regenerative medicine. Nevertheless, implantation of ES cell derived differentiated cells holds the risk of teratoma formation due to residual undifferentiated cells. In order to tackle this problem, we used pluripotent stem cells consisting of ES and iPS cells of mouse genetically modified by lentiviral vectors (LVs) carrying herpes simplex virus thymidine kinase (HSV-TK) under the control of different promoters of pluripotency genes. Cells expressing TK in turn are eliminated upon administration of the prodrug ganciclovir (GCV). Our aim was to study the conditions required for a safe mechanism to clear residual undifferentiated cells but using low MOIs of lentiviruses to reduce the risk of insertional mutagenesis. Our in vitro data demonstrated that TK expression in pluripotent stem cells upon treatment with GCV led to elimination of undifferentiated cells. However, introduction of hygromycin resistance in the LV transduced ES cells followed by pre-selection with hygromycin and GCV treatment was required to abolish undifferentiated cells. Most importantly, transplantation of pre-selected ES cells that had been transduced with low MOI LV in mice resulted in no teratoma development after GCV treatment in vivo. Taken together, our data show that pre-selection of ES cells prior to in vivo application is necessary if vector integration events are minimized. The study presented here gives rise to safer use of pluripotent stem cells as promising cell sources in regenerative medicine in the future.

Introduction
Embryonic stem (ES) cells are derived from the primitive ectoderm of the inner cell mass of blastocysts [1,2]. They are characterized by their self-renewal capability and their pluripotency, i.e. they can develop into the three primary germ layers (ectoderm, endoderm, mesoderm) [3]. Because of their capacity to differentiate into all cell types of the adult body, ES cells became a promising source for cell-based therapies for regenerative medicine over the past years. However, the application of differentiated pluri- [4] or multipotent stem cells [5] for these approaches carries a potential risk of tumor (teratoma) formation in vivo due to residual undifferentiated cells in the transplanted cell population. Hence, removal of residual undifferentiated stem cells from the differentiated cell population has been considered as an essential requirement for use of stem cell-based therapies. In the light of ethical controversies around the usage of human ES cells, a number of groups demonstrated successful generation of induced pluripotent stem (iPS) cells from adult somatic cells [6][7][8]. Thus, iPS cells might also be used as an alternative source for stem cells in regenerative medicine or cell replacement therapies [8]. Also for these cells safety concerns about their tumorigenic potential have to be addressed.
Previous reports have proposed elimination of the undifferentiated cells using suicide genes [9][10][11]. The transfer of a suicide gene has even been successfully used in clinical trials for tumor elimination [12]. One of the most thoroughly studied and widely used approach is based on the herpes simplex virus thymidine kinase (HSV-TK) that converts the prodrug ganciclovir (GCV) to a toxic metabolite [12]. Various routes to deliver the transgene, including transfection or viral transduction, have been studied [10,11]. Moreover, approaches using cytotoxic antibodies against undifferentiated ES cells [13,14] or an antibody against a surface antigen of ES cells combining flow cytometry-based separation were used to remove undifferentiated pluripotent cells [15] before cell transplantations.
Lentiviruses are members of the Retroviridae family, which can stably integrate their genetic information into the host genome of dividing as well as non-dividing cells [16,17]. HIV-1 is the best studied lentivirus and most of the currently used lentiviral vectors (LVs) are based on its sequence [16,[18][19][20]. Previous studies demonstrated that LVs allow for an efficient gene transfer in ES cells [21,22]. In addition, LVs have already been applied in first clinical gene therapy trials (e.g. [22][23][24]).
In the present study, we utilized LVs for the genetic modification of ES and iPS cells of mouse. To enable TK expression in undifferentiated pluripotent stem cells only, different promoters of pluripotency genes were used including Oct-3/4 [25,26], Nanog [11,27,28], EOS-C3 [29] or EOS-S4 [29]. Cells expressing TK are sensitive to GCV treatment. Using this approach, we successfully eliminated undifferentiated cells in vitro. However, to obtain a pure population of TK expressing ES cells, a hygromycin resistance cassette had to be incorporated in the LV. This led to elimination of undifferentiated cells after hygromycin preselection. Most importantly, in vivo transplantation of these LV transduced pre-selected ES cells led to loss of teratoma formation upon GCV administration to the mice.
For analysis of ES cell survival under undifferentiation conditions, cells were transduced with LVs (see below) and treated with or without GCV. Surviving undifferentiated cells were manually counted using three different fields of view that were counted twice.
For differentiation of ES cells into embryoid bodies (EBs) the mass culture protocol was used [30]. Briefly, lentiviral transduced or untransduced ES cells were split to single cells in differentiation medium and cell suspension was incubated at 37°C with ~80 rpm shaking for 14 days. The differentiation medium contained nonessential amino acids (0.1 mM), Lglutamine (2 mM), penicillin 100 (units/ml), streptomycin (100 µg/ml), β-mercaptoethanol (0.1 mM) and fetal calf serum (20% (v/v)) in IMDM (Gibco). GCV treatment was started after day 8. For analysis of ES cell survival under differentiation conditions, EBs were dissociated and stained as mentioned below (Immunocytochemistry). Oct3/4 positive cells were manually counted using one to three different fields of view that were counted five times, respectively.
For in vivo ES cell injection, undifferentiated ES cells were trypsinized and incubated on 10 cm tissue plates coated with 0.1% gelatin for 30 minutes. This enables a separation of feeder cells (reattachment on the plate) and ES cells (remaining in suspension). Afterwards, 1x10 6 ES cells/ 100 µl suspension were prepared in DMEM without LIF and kept on ice before ES cell injection.

Production of LVs
The production of LVs was performed as previously described [36]. Briefly, LV plasmids as well as the structure and packaging plasmids pMDLg/pRRE [37], RSV-rev [37] and pMD2.G [18] were co-transfected into HEK 293T cells (ATCC) seeded on Poly-L-lysine coated 150-mm 2 dishes. Cells were incubated at 37°C and 3% CO 2 overnight. Medium was changed 12 hours later. 48 hours and 72 hours after transfection, the released viruses in the supernatant were collected and cell debris was removed using a SFCA bottle-top filter. Supernatant was centrifuged by an ultracentrifuge with SW32Ti rotor at 61,700g at 17°C for 2 hours. Subsequently, supernatant was discarded and virus pellet was resuspended in Hank's Balanced Salt Solution (HBSS), respectively. Combined virus suspensions were concentrated by centrifugation over a 20% (w/v) sucrose cushion in a SW55Ti rotor at 53,500g at 17 °C for 2 hours. Finally, viruses were dissolved in HBSS and stored at -80 °C. The titers of the virus preparations were determined by measuring the amount of active viral reverse transcriptase (RT) using a colorimetric RT-ELISA (Roche) with a calibration curve.

LV transduction of ES cells
1x 10 4 ES cells were seeded in a 24-well plate with 1.2x10 5 feeder cells/well (see also section cell lines and cell culture). After ES cells attached on the feeder cells, medium was removed and replaced with 300 µl of medium containing LVs overnight. Next day, medium was replaced by fresh medium and ES cells were further cultured for experiments. Before complete confluence, the ES cells were split according to experimental procedures.

Determination of provirus integration per genome using quantitative Realtime PCR (qPCR)
Genomic DNA (gDNA) of ES cells was isolated by classic phenol-chloroform extraction. After isolation, DNA was digested with BamHI and EcoRI to obtain smaller DNA fragments that are more suitable as template in qPCR. The qPCR analysis for quantification of provirus integration was performed using the IQ5 Real-time PCR System (Bio-Rad) and the TaqMan method (iQ Multiplex Powermix, Bio-Rad) in accordance to the manufacturer's instruction. To detect the amount of integrated provirus DNA, Late-RT primers 5'-TGTGTGCCCGTCTGTTGTGT -3' (forward), 5'-GAGTCCTGCGTCGAGAGAGC -3'(reverse) and the 5'-FAMlabeled probe 5'-FAM-CAGTGGCGCCCGAACAGGGA-BHQ1-3' were used [38]. Primers and the Texas Red labeled probe of the house keeping gene Burkitt lymphoma receptor 1 (BLR1) were used to detect total gDNA as described elsewhere [39]. Since BLR1 is a single-copy gene, it can be used to determine absolute copy number per genome by comparison of the cycle threshold (Ct) values from virus specific probe versus the BLR1 probe. Copy numbers were evaluated by following formulas: Copy number/genome = 2 ΔCt with ΔCt = Ct (BLR1) -Ct (LateRT).

Generation and analysis of ES cell clones
1x10 4 ES cells were seeded on feeder cells in a 24 well plate and transduced with lentiviral vectors (NT or OT). 72h after transduction ES cells were trypsinized and 500 ES cells were seeded on 100-mm tissue dishes coated with feeder cells (4x10 6 feeder cells per dish). 7 days later, ES cell colonies were scraped carefully with a pipette tip to dislodge the colonies from the feeder cell layer. Subsequently, ES cell colonies were transferred into a well of a 24-well plate containing 120µl 0.05% trypsin/EDTA. After incubation at 37°C for 5 minutes, the cells were resuspended gently to dispense to single cells and split on three new 24 wells and further cultivated until they reached confluency. When reaching confluency, one well was used for freezing, one well for further cultivation and subsequent experiments and one well for gDNA isolation by classic phenol-chloroform extraction. After isolation, gDNA was digested with BamHI. Digested gDNAs of different clones were loaded on an agarose gel and further tested by using Southern Blot procedure with lentiviral integration specific probe (WPRE probe, amplified from a lentiviral construct using the following primers: 5'-GGACGAACATGGCGTGAAGC-3' (forward) and 5'-GGACGTCCCGCGCAGAAT-3' (reverse)). Experiments presented in the manuscript were performed with four clones (NT #8, NT #11, OT #4, OT #11) with only one integrant (one band on southern blot).

LDH assay
The effect of GCV treatment of TK transduced cells on cytotoxicity was analyzed using the LDH assay (Roche). The LDH in cell culture supernatants released from dead cells was either directly analyzed using LDH assay or LDH was measured after complete lysis of cells of each set of samples: 3x10 3 ES cells transduced with NT, OT, CT and ST or untransduced cells were seeded on a 24-well plate coated with a layer of feeder cells. One day after ES cell seeding, the medium was replaced by 500 µl fresh medium with 20 µM GCV. 36 hours later, further 500 µl fresh medium with 20 µM GCV was added on the cells. Further 36 hours later (=72 hours with GCV treatment), 50 µl of the supernatants from ES cells with GCV treatment and ES cells with GCV treatment after adding of lysis buffer for 15 minutes were each analyzed using the LDH assay in accordance to the manufacturer's instructions. To obtain "background" signals feeder cells without ES cells were seeded as well and treated with GCV. 50 µl of the supernatants from feeder cells alone and feeder cells alone after adding of lysis buffer for 15 minutes were analyzed. These values with or without lysis were each subtracted from the extinction values obtained from ES cells. The differences between values with or without lysis represent the number of ES cells that survived GCV treatment. Extinction was each measured using ELISA reader (TECAN).
The assay was applied for the quantification of ES cell survival under non-differentiation conditions.

Injection of ES cells in murine hind limbs
100µl cell suspension with 1x10 6 ES cells (see also section cell lines and cell culture) were injected s.c. into the hind limbs of 4 to 8 weeks old Fox, Chase SCID® Beige mice (Charles River). Three hours or seven days after ES cell injection, mice were administered i.p. with either saline solution (0.9% (w/v) NaCl) or GCV (20 mg/kg/day) for 12 days. After treatment, mice were sacrificed for the analysis of potential teratoma formation. The teratomas were analyzed by measuring their weight and size. For further analysis, teratomas were fixed in 4% (w/v) paraformaldehyde (PFA) in phosphate-buffered saline (PBS), embedded in O.C.T. (Tissue-Tec) and stored at -80 °C. Sections were prepared in 10 µm slices for H&E staining.
All experiments involving animals were done in accordance with the German legislation on protection of animals and the NIH Guidelines for the Care and Use of Laboratory Animals. The study was approved by the Landesamt für Natur, Umwelt und Verbraucherschutz Nordrhein-Westfalen, Germany. Cell injection was performed under isofluran inhalation anaesthesia, and all efforts were made to minimize suffering.

Hematoxylin and eosin stain
Teratoma sections were incubated in PBS for 5 minutes two times. Afterwards, slides were treated with hematoxylin (Mayers hemalaun, Merck) for 1 minute and then washed in fluent water for 5 minutes. After that, slides were stained with eosin (Eosin G, Merck) for 30 seconds and washed in fluent water for 5 minutes. Finally, slides were mounted with Roti®-Histokitt (Carl Roth). Images were taken using an AxioStart (Carl Zeiss) and AxioCam (Carl Zeiss).

Immunocytochemistry
EBs (see also section cell lines and cell culture) were collected and washed twice with PBS. After that, EBs were incubated in Collagenase B (1 mg/ml) (Roche) for 20 minutes in a shaker with 600 rpm at 37°C. 3x10 5 single cells from dissociated EBs were seeded in a 24-well plate with cover slips which were coated with 0.1% (w/v) gelatin. After reattaching of the cells on the cover slips (24 hours later), cells were fixed with 4% (w/v) PFA. Cells fixed on the cover slips were incubated with primary antibodies in solution containing 0.05% (v/v) donkey serum (Jackson ImmunoResearch), 0.1% (v/v) triton, anti Oct-3/4 (1:100, Santa Cruz Biotechnology) and anti α-actinin (1:400, Sigma-Aldrich) in PBS for 2 hours. After removing the primary antibodies and washing three times with PBS, cells were further incubated with secondary antibodies containing anti Dylight 549 and anti Dylight 649 (both 1:400, Jackson ImmunoResearch) in Hoechst stain solution (Sigma-Aldrich) for 1 hour. Images were taken using Axio Observer. Z1 with ApoTome system (Carl Zeiss) and AxioCam MRC5 (Carl Zeiss).

Statistical analysis
All values mentioned in the text are given as mean±SEM. Statistical significance was determined by one-way analysis of variance (ANOVA) or by using Student's t-test. We considered p<0.05 as statistically significant.

I: LV-mediated TK expression in pluripotent stem cells and GCV treatment
In the present study, we used a murine ES cell line that carries the cDNAs for puromycin resistance and green fluorescent protein (eGFP) connected via the IRES (internal ribosomal entry site) element. Transgene expression is driven by the cardiac specific α-myosin heavy chain (α-MHC) promoter. Thereby, successful differentiation of ES cells can be assessed by fluorescence microscopy after applying culture conditions that favor cardiac differentiation [40].
In order to generate ES cells stably expressing TK only in undifferentiated cells [25][26][27][28], we first applied LVs carrying HSV-TK cDNA under control of promoters of the pluripotency genes Nanog or Oct-3/4 (NT and OT) ( Figure 1A). Indeed, treatment of the transduced cells with GCV for 72 hours had a toxic effect on undifferentiated ES cells ( Figure 1B). Using different concentrations of GCV, we found that treatment with 20 µM GCV was sufficient to eliminate undifferentiated cells but had no adverse effect on wildtype ES cells ( Figure 1B). Therefore, 20 µM GCV was used in all further in vitro experiments.
As iPS cells also carry the potential risk of tumor formation after transplantation in in vivo applications, removal of iPS cells that did not undergo differentiation is an important task. Thus, we tested the NT-and OT-LVs in two different mouse iPS cell lines: iPS and iPS-Oct-GFP, the latter containing an eGFP expression cassette under control of promoter of the pluripotency gene Oct-3/4. Upon GCV administration, death of NT-and OT-transduced iPS-cells could be observed ( Figure  1C). Analysis of Oct-3/4 promoter driven eGFP expression in iPS-Oct-GFP cells also revealed efficient killing in the NT-or OT-transduced cells after GCV treatment ( Figure 1D). Taken together, besides murine ES cells, the TK/GCV system also led to an efficient elimination of undifferentiated iPS cells in vitro.

II: LV mediated TK expression in mixed ES cell populations
For clinical applications, it is desired to use a low copy number of LVs to minimize the potential risk of insertional mutagenesis by vector integration [41]. Therefore, we investigated ES cell populations with a low average copy number (e.g. one to two provirus integrations) per genome. A recent report showed that quantitative real time PCR (qPCR) can be used as a reliable method to determine the copy number of LVs [39]. Hence, we transduced ES cells with different LV concentrations and analyzed their average integration number using qPCR. These LV transduced ES cell populations without additional selection are further referred to as mixed ES cell populations. Mixed ES cell populations with an average copy number of 1.5 provirus integrations per genome (1.60±0.09 (NT) and 1.38±0.24 (OT)) were chosen for further analysis (Figure 2A). GCV treatment of these undifferentiated populations inhibited cell survival as shown by bright field images ( Figure 2B). Next, we investigated the effect of GCV in mixed ES cell populations under differentiation conditions. We differentiated untransduced as well as NT-or OT-transduced mixed ES cell populations using the mass culture protocol [30] for 14 days. The resulting embryoid bodies (EBs) were dissociated and characterized by immunostaining with antibodies against Oct-3/4 (pluripotency marker) and Hoechst (nuclear stain) ( Figure 2C Figure 2D). In the untransduced ES cells, we observed 13.1%±2.2 Oct-3/4 positive cells after GCV treatment indicating that GCV alone had no significant effect on cellular vitality ( Figure 2D). In contrast, only 3.6%±1.6 and 3.0%±0.9 of the NT-or OT-transduced mixed ES cell populations were Oct-3/4-positive when treated with GCV ( Figure 2D). Thus, most of the undifferentiated ES cells are eliminated during differentiation in the presence of NT-or OT-driven TK expression.
Differentiation of ES cell-derived cardiomyocytes was not influenced neither by the LV transduction nor the GCV treatment alone as shown by cardiac-specific α-MHC promoter driven eGFP expression (Supplementary Figure S1).
Taken together, ES cells transduced with NT-and OT-LVs at low copy numbers exhibited high killing efficiency upon GCV treatment after differentiation. Nevertheless, the mixed populations still contained around 3 to 4% Oct-3/4 positive cells on day 14 post differentiation ( Figure 2D) suggesting a potential risk of teratoma formation in vivo. To address whether higher levels of TK expression lead to more efficient elimination, we tried to increase TK abundance by using other promoters.

III: LV mediated TK expression in mixed ES cell populations using EOS promoters
In order to achieve higher levels of TK expression, we applied the promoters of the pluripotency genes EOS-C3 and EOS-S4, shown to drive stronger expression than Nanog or Oct-3/4 promoters in human or murine ES and iPS cells [29]. The promoters were cloned into the LV backbone leading to LV-constructs EOS-C3-TK (CT) and EOS-S4-TK (ST) ( Figure  3A). The four LVs (NT, OT, CT, ST) were used in parallel for transduction of murine ES cells and LV integration was determined by qPCR, respectively. Mixed ES cell populations with an average copy number of approximately 1.5 per genome were further investigated ( Figure 3B). TK expression on mRNA level was analyzed using qPCR ( Figure 3C). These experiments revealed that both, EOS-C3 and EOS-S4, promoters gave rise to higher TK expression levels compared to Nanog and Oct-3/4 in murine ES cells ( Figure 3C). However, the highest TK abundance was detected using EOS-S4 promoter ( Figure 3C). We observed a clear reduction of surviving ES cells in all four cell populations after GCV treatment ( Figure 3D). Next, we analyzed the TK/GCV induced killing efficiency by detecting Lactate Dehydrogenase (LDH) that is released by dead cells into the cell culture supernatant. Interestingly, a higher TK expression level did only slightly increase the effect of GCV on ES cells ( Figure 3E) demonstrating that the portion of surviving cells seems to be the consequence of incomplete transduction of the ES cell population.

IV: LV mediated TK expression in single integrant ES cell clones and GCV treatment
The mixed ES cell populations analyzed so far still resulted in incomplete elimination of undifferentiated cells after GCV treatment. However, a few pluripotent cells could lead to formation of teratomas in vivo [42]. To study whether a complete abolishment of undifferentiated cells could be reached in cells carrying one integrant, we next analyzed clones derived from individual cells from single integrant NT or OT expressing ES cell lines. Single ES cell clones that carry only one integrant per genome were identified by Southern blot analysis ( Figure S2A, B). Two different single-integrant clones of each NT-and OT-transduced ES cells were used for further investigations (NT #8, NT #11, OT #4, OT #11). Under undifferentiation conditions, cells of all four clones were eliminated upon GCV treatment ( Figure 4A). The differentiation of the ES cells of all four clones into ES cell-derived cardiomyocytes was not influenced by the LV transduction or the GCV treatment ( Figure S2C). Importantly, under differentiation conditions, only 1.5%±0.7 (NT #8), 1.5%±0.6 (NT #11), 0.5%±0.2 (OT #4) and 0.3%±0.3 (OT #11) of surviving Oct-3/4 positive cells were observed on day 14 post differentiation after GCV treatment for 6 days ( Figure 4B,C). These data demonstrated a much higher efficiency of the TK/GCV system in single integrant ES cell clones compared to the mixed cell populations where 3 to 4% of Oct-3/4 positive cells were observed under differentiation conditions (see also Figure 2D).

V: LV mediated TK expression and GCV treatment in mixed ES cell populations using Hygromycin preselection
The use of single integrant NT-and OT-transduced ES cell clones resulted in a more efficient elimination of undifferentiated cells after GCV treatment (see also Figure 4C) as compared to mixed ES cell populations (see also Figure  2D). Nevertheless, as a complete abolishment of undifferentiated cells is required, these approaches are not efficient enough. For excluding the untransduced ES cells from mixed populations an antibiotic pre-selection scheme was chosen: We incorporated a hygromycin resistance gene expression cassette in the ST construct leading to construct EOS-S4-TK-PGK-Hygromycin (STPH) ( Figure 5A). The hygromycin resistance gene was under the control of the ubiquitous PGK promoter enabling hygromycin selection of successfully transduced ES cells of each state.
The STPH vector was then used for ES cell transduction and mixed ES cell populations with an average copy number of 1.5 per genome were treated with GCV after hygromycin pre-selection for 9 days. The GCV treatment of the unselected, undifferentiated STPH-transduced ES cells resulted in a significant reduction of ES cells (47.3%±0.4 ES cell survival) whereas untransduced cells showed no difference with or without GCV as shown by brightfield images ( Figure 5B) and by cell counting (Figure 5C, D). Importantly, with hygromycin pre-selection, only 0.4%±0.3 undifferentiated STPH-transduced ES cells were observed in the presence of GCV indicating pre-   treatment is required. In addition to hygromycin selection for 9 days ( Figure S3B see also Figure 5E), we tested pre-selection for 6 days followed by 3 days without selection and subsequent GCV treatment ( Figure S3C) as well as 9 days pre-selection followed by 3 days without hygromycin before administration of GCV ( Figure S3D).
Importantly, treatment for 9 days with hygromycin followed directly by GCV administration resulted in an identical ES cell survival as the hygromycin treatment for 9 days followed by 3 day incubation without selection and subsequent GCV administration (0.15%±0.08 vs. 0.15%±0.15, respectively). The results indicate that after hygromycin selection for 9 days hygromycin can be withdrawn before GCV selection. On the other hand, if the duration of hygromycin pre-selection is reduced to 6 days and cells are subsequently cultivated for 3 days without antibiotics, 3.66% ±0.63 of ES cells survive after GCV application. Consequently, shorter pre-selection seems not to be sufficient to ensure survival of only transduced cells. Taken together, after 9 days of pre-selection a subsequent removal of selection pressure does not alter the efficacy of GCV-mediated elimination of undifferentiated cells. Alternatively, we speculated that an increase of lentiviral copy number might result in a similar efficiency: We transduced ES cells with higher LV concentrations of STPH leading to an average copy number of 3.8. As high copy numbers increase the risk of insertional mutagenesis we did not further raise the average copy number. GCV treatment of these high copy ES cells without pre-selection revealed strong elimination of undifferentiated cells (35.4%±2.0 ES cell survival) as shown by brightfield images ( Figure S4A) and cell counting ( Figure S4B, C). These data demonstrate that increasing the copy number significantly increases the elimination efficiency of TK/GCV system, but complete killing of undifferentiated cells was not achieved.
Importantly, when differentiating STPH-transduced ES cells with low copy number, that were pre-selected with Hygromycin for 9 days, we observed no Oct-3/4 positive cells after GCV treatment ( Figure 5F).
In conclusion, STPH-transduced ES cells of low copy number with hygromycin pre-selection for 9 days resulted in successful elimination of ES cells upon GCV treatment. Although, increasing the lentiviral copy number showed lower cell survival of undifferentiated ES cells without pre-selection the complete elimination of ES cells was not achieved and thus, pre-selection seems to be necessary.

VI: In vivo application of TK expressing ES cells
In the next step, we applied the LV mediated TK/GCV system tested so far in vitro in an in vivo mouse teratoma model to investigate tumor formation due to undifferentiated cells [42]. Therefore, we injected 1x10 6 untransduced or STPHtransduced undifferentiated ES cells after hygromycin preselection s.c. into the left and right hind limbs of 8-week old SCID® mice, respectively ( Figure 6A). Three hours after ES cell injection, we administered mice i.p. with either saline solution (0.9% (w/v) NaCl) or with 20 mg/kg GCV per day for 12 days. Teratoma will only form of undifferentiated cells that survived the selection.
Three weeks after ES cell injection mice were sacrificed for analysis of teratoma formation: No significant difference in the weight of teratoma was detected in mice treated with saline solution after injection of untransduced or STPH-transduced ES cells after pre-selection ( Figure 6B). In addition, the average size of teratomas (maximum dimension) was more than 1 cm, respectively ( Table 1). Sections of the teratomas exhibited tissues from all three germ layers as representatively shown for STPH-transduced and pre-selected ES cells without GCV treatment (Figure 6C upper). The same results for the teratoma formation were obtained after injection of untransduced cells and treatment with GCV ( Figure 6B and Table 1) indicating no influence of GCV or hygromycin alone. Most importantly, when pre-selected STPH-transduced ES cells were applied, GCV treatment of the mice did not result in tumor formation ( Figure 6B): only one out of six mice revealed the formation of a small (~0.3 cm) piece of muscle tissue that exhibited no other germ layer tissue (Table 1 and Figure 6C lower). These results clearly demonstrate that lentiviral mediated TK expression and GCV treatment can be used to eliminate even a large number of undifferentiated pluripotent cells in vivo.
In addition, we also tested, high copy number STPHtransduced ES cells without pre-selection in vivo. Compared to the control group, receiving untransduced cells, a significant reduction in teratoma weight (around 65%) was observed for high copy number STPH-transduced ES cells without preselection when treated with GCV ( Figure 6D). Furthermore, four out of five mice showed formation of only small teratomas upon GCV treatment (Table 1). In contrast, upon saline treatment all mice that were injected with STPH-transduced unselected ES cells, developed teratomas with maximum dimensions of more than 1 cm (Table 1). However, the complete absence of teratoma as observed for injection of preselected low copy number STPH-transduced ES cells was not achieved when using ES cells with high copy number LV TK (see also Figure 6B). These results confirmed our in vitro findings (see also Figure 5E and Supplementary Figure S4C). Thus, increasing the copy number of LVs does not per se result in a complete transduction of ES cell populations.
The data obtained suggested a promising future for the in vivo application of STPH-transduced ES cell lines with preselection as ES cells that have teratoma-forming potential seemed to be completely eliminated due to GCV application.   In further experiments, we started GCV treatment 7 days after ES cell injection ( Figure S5). First, we analyzed tumor growth 7 days after ES cell injection of mice that were not treated with GCV ( Figure S5A). In the absence of GCV, both untransduced as well as STPH-transduced ES cells gave rise to teratomas with no significant difference (teratomas size of approximately 4 mm). These results clearly showed the formation of teratomas within 7 days if no further treatment was applied. Hygromycin selection per se does not affect teratoma development of ES cells. Next, we compared teratoma formation of saline or GCV-treated mice that were injected with not transduced or STPH-transduced ES cells ( Figure S5B). The treatment started 7 days after cell injection, when small tumors already had developed (see also Figure S5A). No significant differences in tumor sizes between injected untransduced and STPH-transduced and pre-selected cells upon GCV treatment could be observed showing that the TK/GCV system is not able to eliminate developed tumors. Our results underline the necessity of GCV application directly after ES cell injection to prevent tumor formation caused by pluripotent cells.

Discussion
Pluripotent stem cells are characterized by their ability to self-renew and capability to differentiate into all cell types of the adult body. Due to these characteristics they are a promising cell source for regenerative medicine. However, former studies demonstrated the risk of teratoma formation in ES cell-based therapies [4] which remains a major safety obstacle that has to be solved before pluripotent cell-derived tissues can be applied in clinical trials. To avoid tumorigenesis of undifferentiated pluripotent cells different approaches have been developed. One method is the use of suicide genes such as the TK/GCV system to eliminate undifferentiated pluripotent cells [9,10,43]. Further possibilities rely on selection of residual pluripotent cells by application of specific surface markers and subsequent cell sorting [15] or on use of a specific antibody for elimination of pluripotent cells [13,14]. However, using surface markers to select residual pluripotent cells can only be applied in vitro. The TK/GCV system on the other hand, leads to elimination of undifferentiated pluripotent cells during the differentiation process in vitro as well as to removal of undifferentiated pluripotent stem cells in vivo. TK/GCV based therapies have already been used in vivo and are tested in clinical trials for e.g. the treatment of prostate cancer [44] or ovarian cancer [45].
To obtain stable transgene expression, the use of integrating LVs is a suitable tool for efficient gene transfer in a broad spectrum of stem cells including murine ES [22] and iPS cells [29] for both in vitro as well as in vivo applications [46]. In the present study, we used LVs carrying different promoters of pluripotency genes (Nanog, Oct, EOS-C3, EOS-S4) upstream of the TK cDNA in order to specifically eliminate pluripotent stem cells. In contrast, using ubiquitous promoters such as the phosphoglycerate kinase (PGK) promoter or the cytomegalovirus (CMV) promoter to drive TK expression would lead to elimination of all transduced cells and is therefore not useful [9].
Integrating viral vectors like LVs carry the risk of insertional mutagenesis [36,47]. Therefore, on the one hand low integration numbers in transduced pluripotent stem cells are desired. On the other hand, this increases the likelihood of untransduced cells that do not express TK and subsequently not all undifferentiated cells are eliminated upon GCV treatment. In the present study, we tried to find possible solutions for this dilemma by applying different approaches. In addition, motivation for this study was establishment of a protocol with ready-to-use LVs in combination with the TK/GCV system that could be applied to different pluripotent stem cell lines. This would facilitate a possible application in clinical trials.
In the first step, we applied TK expressing LVs to both ES and iPS cells to obtain infected cell populations that were achieved without further selection or clonal expansion. This resulted under undifferentiation conditions for both stem cell populations in elimination of cells after GCV treatment. We focused on ES cells that were transduced with certain amounts of LVs leading to a relatively low average copy number of approximately 1.5 per genome. Treatment of these NT-or OTtransduced mixed ES cell populations with GCV resulted in a survival of around 3 to 4% Oct-3/4-positive cells after differentiation. The presence of Oct-3/4-positive cells could be due to an insufficient TK expression level of OT-or NTtransduced ES cells or due to remaining untransduced ES cells. Analysis of different promoters to raise TK expression level revealed that EOS-S4 led to the highest TK abundance among the promoters tested. However, we could not observe higher killing efficiency of undifferentiated cells upon GCV treatment in vitro. Therefore, it is more likely that survival of undifferentiated cells is the consequence of untransduced cells that remain unaffected upon GCV treatment. To ensure that all ES cells contain the TK cDNA, we generated LV single ES cell clones by picking cell clones and further dispending the colonies. The clones were screened using Southern Blot analysis with a lentiviral specific probe and clones with a single integration were further investigated. Whereas almost no survival of undifferentiated cells was obtained for differentiated single integrant clones OT #4 and OT# 11, more than 1% of Oct-3/4 positive cells were still observed after GCV treatment of NT clones (NT #8 and NT #11). This incomplete elimination might be due to an epigenetic silencing of the transgene [29,48]. Another explanation for the presence of Oct-3/4 positive cells after GCV treatment and differentiation might be contaminations by a few untransduced ES cells (i.e. mixed cell clone). To avoid this, the single integrant clones can be further subcloned to achieve pure single integrant ES cell clones but this is more time consuming. Furthermore, in contrast to the generation of single clones, the usage of mixed ES cell populations is much simpler to handle since the additional process of clone selection and screening by Southern blot is not necessary. Irrespectively whether the TK vector is not functional (silencing) or not present (contamination with untransduced cells), a pre-selection strategy can be employed to obtain pure ES cell populations that express TK. Therefore, we incorporated hygromycin resistance gene as a pre-selection tool in the LVs (STPH vector) to select for transduced ES cells.

And indeed, treatment of STPH-transduced ES cells with hygromycin and GCV abolished undifferentiated cells in vitro.
Importantly, no teratoma formation in the in vivo studies could be observed when mice were injected with these pre-selected transduced ES cells and treated with GCV. This could be a promising application of LVs and the TK/GCV system in future clinical approaches. The GCV treatment has to be applied directly after ES cell injection to prevent tumor formation caused by pluripotent cells as it is not able to eliminate or stop the growth of established tumors. This in accordance with previously published results [49].
In addition, we tried to raise copy numbers of STPHtransduced ES cells to decline the untransduced cell population. Even with higher average copy numbers (approximately 3.8), several ES cells seemed to be not transduced because undifferentiated cells survived the GCV treatment. In conclusion, complete elimination was only achieved by hygromycin pre-selection.

Conclusions
In summary, we successfully generated pluripotent stem cells expressing the TK gene after LV transduction. Our in vitro and in vivo data clearly demonstrated an efficient killing effect of undifferentiated pluripotent stem cells by GCV treatment by using low copy number STPH-transduced ES cell pre-selected with hygromycin. This in turn eliminated the risk of teratoma formation in the hind limbs of GCV treated mice in vivo. Moreover, our in vitro data also suggests interchanging of iPS cells for ES cells averting the controversial ethical application of human ES cells. Figure S1. Immunostaining of dissociated EBs derived from mixed ES cell populations transduced with lentiviral NT or OT with low copy number. The ES cells transduced with NT or OT leading to 1.5 copy numbers per genome in average (as shown in Fig. 2B) were differentiated as EBs and treated with (+) or without (-) 20 µM GCV. After 14 days of differentiation dissociated EBs were immunostained with antibody against skeletal α-actinin (sarcomere) indicating the successful formation of cardiomyocytes (white); green fluorescence indicates eGFP expression driven by cardiac specific promoter, α-MHC. (TIF) .001 compared to without GCV treatment, respectively, Student's t-test). Data is based on images representatively shown in (A), undifferentiated cells were manually counted using three different fields of view that were counted twice. (TIF) Figure S5. Analysis of injection of untransduced or STPHtransduced ES cells with delayed GCV treatment. 1 x 10 6 ES cells (STPH-transduced with 1.5 copy numbers per genome in average and pre-selected with hygromycin for 9 days or not transduced (WT)) were injected s.c. into hind limbs of SCID/ beige mice. (A) Relative weight of teratomas was analyzed 7 days after cell injection (n=3, Mean±SEM). (B) 7 days after cell injection mice were administrated i.p. with saline solution (0.9% (w/v) NaCl; -GCV) or with GCV (20 mg/kg/day, +GCV) for consecutive 12 days and relative weight of teratomas was finally analyzed (n=6, Mean±SEM). (TIF)


s2fieldsofstudy: Biology Medicine

text:
5-Fluorouracil-Induced Changes of Intestinal Integrity Biomarkers in BALB/C Mice

Background: Intestinal mucositis is a most frequently occurring toxicity in cancer chemotherapy, and consequent malnutrition reduces tolerance to cancer therapies. Therefore it is important to lessen the severity of mucotitis and to develop complementary agents capable of reducing mucotitis-related symptoms. This study was conducted to determine 5-fluorouracil (5-FU) induced intestinal damage to understand intestinal damages due to chemotherapy and to provide information on biomarkers which can be used to screen complementary agents in future studies. Methods: BALB/c mice were divided into three experimental groups and subjected to the intraperitoneal injection of either 100 mg/kg or 200 mg/kg of 5-FU. The third group was used as PBS controls. Body weights and the consistency of the stools were recorded every day, and the animals were sacrificed on the 7th day post 5-FU administration. The expressions of intestinal tight junction proteins and mRNAs of tumor necrosis factor-alpha (TNF-α) and interleukin-1 beta (IL-1β) were determined. Results: The body weight of the animals treated with 5-FU was significantly decreased in a dose-dependent manner. However, mice given 100 mg/kg 5-FU rapidly recovered the original body weight. Symptom of diarrhea was also more severe in 200 mg/kg 5-FU treated group than that of the 100 mg/kg 5-FU treated animals. The expressions of occludin and claudin-1, not ZO-1 protein expressions in 200 mg/kg 5-FU treated animals were significantly reduced compared to those of the control group or 100 mg/kg 5-FU group. The expression of Nuclear factor-kappa B p65 (NF-κB p65) protein and TNF-α mRNA were significantly higher in 5-FU treated group compared to those of control group. No difference was observed with IL-1β expression. Conclusions: These results suggested that selected tight junction proteins and inflammatory cytokines are related to 5-FU induced mucositis, and thereby can be used as targets of developing complementary agents.

INTRODUCTION
Gastrointestinal mucositis is a common side effect of anticancer chemotherapy that involves the small intestines. Not only does mucositis decrease the quality of life in most cancer patients because of its associated intense pain, it is also a high-risk factor for hematosepsis with neutropenia and malnutrition. This association, thus, renders mucositis a clinically important disease. Moreover, this specific side effect does not allow appropriate chemotherapy and radiotherapy, and becomes an obstacle that delays the administration of effective cancer treatment. 1 5-Fluorouracil (5-FU) is an anticancer drug that is widely used in the treatment of colorectal cancer, and it possesses a chemical structure similar to that of uracil and thymine. 2 A major part of it (more than 80%) is catabolized by hepatic dihydropyrimidine dehyrogenase (DPD), whereas the rest induces cell death by inhibiting RNA and DNA syntheses through fluorodeoxyuridine monophosphate (FdUMP) and fluorouridine triphosphate (FUTP). 3,4 The major side effects of 5-FU include myelosuppression, diarrhea, cardiotoxicity, dermatitis, and mucositis. 5,6 Of these, gastrointestinal mucositis has been reported in approximately 80% of patients who have received cancer treatment by using

Sample collection and pretreatment
Once the treatment period was completed, the mice were

Diarrhea assessment
The presence or absence and degree of diarrhea were determined based on the standards established for the diarrhea assessment. 8 The diarrhea score was assessed  The extracted proteins were quantified using a Protein   quantification values and 2 -ΔΔCt were calibrated using β -actin, which is a housekeeping gene. Each experiment was performed in triplicate.

Statistical analysis
For the results of all experimental analyses, means and standard deviation in each group were calculated using a SAS program. After performing one-way ANOVA, statistical significance of the means in each group was tested using Duncan's multiple range test, at a significance level of α= 0.05.

Body weight changes and main organ weight
Body weight changes in the experimental animals are shown in Fig. 2. Compared to the control group, the extent of body weight reduction in the experimental groups manifested significant differences from Day 1 after the 5-FU administration ( Fig. 2A, P＜0.05). Significantly higher reduction in the body weight was observed in the group administered with 200 mg/kg of 5-FU compared to the 100 mg/kg group (Fig. 2A). The weights of the liver and kidney and the lengths of the small and large intestines are presented in Table 2. The weight of the spleen significantly decreased in the 5-FU treated experimental groups (P＜0.05). No significant differences in the weight of the liver and lengths of the small and large intestines were observed.

Diarrhea assessment
The  responses. NF-κB protein expression in the mucous of the small intestine was significantly higher in both the groups compared to that in the control (Fig. 4, P＜0.005). The expression level was slightly lower in the group in which high concentration of 5-FU was used compared to that in the low concentration (Fig. 4, P＜0.05).  permeability is shown to result in the inflow of endotoxins, various pathogenic microorganisms, and varied high molecular substances. 19 A previous study suggested that occludin knockout mice exhibited not only chronic inflammation and hyperplasia of the gastric epithelium but abnormal changes in testicule and salivary gland suggesting more complex role of occludin. 20 A previous study suggested that hydrogen peroxide induced TJ disruption is mediated by a tyrosine kinase dependent mechanism. 21 Claudin is a protein structurally discriminated from occuldin. Claudins are suggested as the backbone of TJ proteins, and the claudin knockout mice revealed abnormal barrier formation and paracellular permselectivity. 22 Our results suggest that 5-FU treatment significantly reduces intestinal membrane TJ protein expression possibly exhibiting stuructural and functional defects. In the case of ZO-1, which is one of the TJ proteins, no significant intergroup differences were observed. Such discrepancies in the TJ proteins may be attributable to their different interactions with other membrane compartments including actin. Although ZO-1 proteins are known to mediate early TJ protein assembly, functional role of ZO-1 has not been clearly demonstrated. Also, a recent report showed that normal TJ structures were formed with ZO-1 deficient cells suggesting further studies on the role of ZO-1 in TJ proteins. 23 Since our results indicate claudins and occuldin expressions are dramatically reduced by 5-FU treatment, it is presumable that severe functional defects are induced by 5-FU treatment and therefore appropriate interventions towards there TJ proteins may alleviate chemotherapy-induced mucositis.

TNF-α and
In conclusion, 5-FU treated animals revealed significant body weight loss in a dose-dependent manner. Clear and significant TJ protein losses are observed suggesting appropriate intervention towards these targets may alleviate chemotherapy-associated mucositis.


text:
DNA-methylation effect on cotranscriptional splicing is dependent on GC architecture of the exon–intron structure

DNA methylation is known to regulate transcription and was recently found to be involved in exon recognition via cotranscriptional splicing. We recently observed that exon–intron architectures can be grouped into two classes: one with higher GC content in exons compared to the flanking introns, and the other with similar GC content in exons and introns. The first group has higher nucleosome occupancy on exons than introns, whereas the second group exhibits weak nucleosome marking of exons, suggesting another type of epigenetic marker distinguishes exons from introns when GC content is similar. We find different and specific patterns of DNA methylation in each of the GC architectures; yet in both groups, DNA methylation clearly marks the exons. Exons of the leveled GC architecture exhibit a significantly stronger DNA methylation signal in relation to their flanking introns compared to exons of the differential GC architecture. This is accentuated by a reduction of the DNA methylation level in the intronic sequences in proximity to the splice sites and shows that different epigenetic modifications mark the location of exons already at the DNA level. Also, lower levels of methylated CpGs on alternative exons can successfully distinguish alternative exons from constitutive ones. Three positions at the splice sites show high CpG abundance and accompany elevated nucleosome occupancy in a leveled GC architecture. Overall, these results suggest that DNA methylation affects exon recognition and is influenced by the GC architecture of the exon and flanking introns.

DNA methylation is known to regulate transcription and was recently found to be involved in exon recognition via cotranscriptional splicing. We recently observed that exon-intron architectures can be grouped into two classes: one with higher GC content in exons compared to the flanking introns, and the other with similar GC content in exons and introns. The first group has higher nucleosome occupancy on exons than introns, whereas the second group exhibits weak nucleosome marking of exons, suggesting another type of epigenetic marker distinguishes exons from introns when GC content is similar. We find different and specific patterns of DNA methylation in each of the GC architectures; yet in both groups, DNA methylation clearly marks the exons. Exons of the leveled GC architecture exhibit a significantly stronger DNA methylation signal in relation to their flanking introns compared to exons of the differential GC architecture. This is accentuated by a reduction of the DNA methylation level in the intronic sequences in proximity to the splice sites and shows that different epigenetic modifications mark the location of exons already at the DNA level. Also, lower levels of methylated CpGs on alternative exons can successfully distinguish alternative exons from constitutive ones. Three positions at the splice sites show high CpG abundance and accompany elevated nucleosome occupancy in a leveled GC architecture. Overall, these results suggest that DNA methylation affects exon recognition and is influenced by the GC architecture of the exon and flanking introns.

[Supplemental material is available for this article.]
A mature mRNA is formed by the removal of introns from the mRNA precursor (pre-mRNA) and the ligation of exons through the process of splicing (Black 2003). Splicing occurs cotranscriptionally, which means that a large fraction of the introns are removed from the mRNA precursor while the transcript is still attached to the DNA by RNA polymerase II (Pol II) (Neugebauer 2002;Proudfoot et al. 2002;Luco et al. 2011). This fact sets the foundation for a cross-talk between the DNA and RNA levels, thereby providing new possibilities for splicing regulation by factors that are known to affect transcription. Pol II can affect splicing through kinetic coupling since the elongation rate of Pol II can control exon recognition by the splicing machinery (de la Mata et al. 2003;Schor et al. 2009;Ip et al. 2011).
One factor that is known to regulate transcription and could also regulate splicing is DNA methylation. DNA methylation is regarded as the ''fifth base'' and is defined as the addition of a methyl group to a cytosine base, predominantly when it is directly followed by guanine but can be found in other contexts as well (CG, CHG, CHH; where H = A, C, or T) (Bernstein et al. 2007). DNA methylation serves a role in many biological processes, including embryogenesis and development, genomic imprinting, and regulation of gene transcription (Li et al. 1992;Okano et al. 1999;Bestor 2000;Reik 2007). Due to dramatic advances in high-throughput DNA sequencing methods, several groups have recently mapped DNA methylation across the whole genome at single-base resolution (Meissner et al. 2008;Laurent et al. 2009;Lister et al. 2009;Gu et al. 2010;Li et al. 2010;Stadler et al. 2011;Smith et al. 2012;Ward et al. 2012;Xie et al. 2012). This new source of high-quality data represents a major step toward our understanding of the biological role of cytosine methylation and has already been used to establish a primary connection between DNA methylation and splicing. In their genome-wide methylation analysis, Lister et al. (2009) revealed that there is a higher level of CpG-methylation on exons than on flanking introns. However, this difference was dissolved when values were divided by cytosine composition since exons have a higher GC content, on average, than their flanking introns (Schwartz et al. 2009). A later genome-wide methylation analysis suggested that DNA methylation is a strong marker of exons and stronger yet for exon-boundaries, thus being a possible regulator of splicing (Laurent et al. 2009). In support of this suggestion, DNA methylation was found to be enriched in alternative splice sites and in splicing regulatory motifs (Anastasiadou et al. 2011), both of which are important regulatory regions for splicing. Furthermore, depletion of DNA methylation in Hox genes was found to remove Pol II stalling and facilitate transcriptional elongation and efficient splicing (Tao et al. 2010). Finally, mutually exclusive DNA methylation and transcriptional repressor CTCF protein binding were found to regulate exon inclusion by influencing Pol II elongation rate (Shukla et al. 2011).
Another epigenetic factor suspected to be involved in both transcription and splicing regulation is chromatin structure. The primary structure of chromatin involves DNA wrapped around nucleosomes. Approximately 147 base pairs of DNA are wrapped around a single nucleosome made up of an octamer of histone proteins and one linker histone (Kouzarides 2007). Positioning of nucleosomes along the genome is partially determined by the DNA sequence and GC content (Segal et al. 2006;Nikolaou et al. 2010) and is subject to modulation by chromatin remodelers (Vignali et al. 2000). DNA methylation, in itself, can act as a chromatin remodeler through two different pathways: either by inducing a rigid nucleosomal conformation that results in a more tightly wrapped nucleosome structure (Choy et al. 2010) or through the binding of methylbinding proteins (MBPs) that silence transcription and modify surrounding chromatin (Sarraf and Stancheva 2004;Klose and Bird 2006). A reverse mechanism, in which nucleosome positioning or histone modifications can affect DNA methylation levels through the recruitment of DNA methyltransferase enzymes was also suggested (Robertson 2002). Several recent publications have shown a higher level of nucleosome occupancy on exons compared to the flanking introns, suggesting a link between nucleosome positioning and splicing regulation (Andersson et al. 2009;Schwartz et al. 2009;Spies et al. 2009;Tilgner et al. 2009;Chen et al. 2010). The nucleosome is thought to act as a roadblock for Pol II elongation (Batsché et al. 2006;Hodges et al. 2009) as Pol II needs to pause and unwind the DNA double strand to release it from the nucleosome before continuing transcription elongation. This transcriptional pausing may allow cotranscriptional recognition of splicing signals in the pre-mRNA (de la Mata et al. 2003).
As nucleosome occupancy is strongly biased toward high GC content (Segal et al. 2006;), which is in high abundance in exons (Schwartz et al. 2009). DNA methylation is also positively correlated with GC content (Bernardi et al. 1985;Jabbari et al. 1997;Varriale and Bernardi 2010). That being said, the level of methylated CpGs over the available CpGs decreases with increasing GC levels ( Jabbari and Bernardi 1998;Oakes et al. 2007), thus painting a more complicated picture. If this is the case, could the roles of these epigenetic modifications in splicing be a side effect of the GC content of the exons? How could one isolate the measured levels of DNA methylation and nucleosome occupancy on exons and conclude a biological role unbiased by GC content?
To evaluate this, we use two populations of exons that are distinguished by their exon-intron GC-content pattern: differential GC, where the exon has higher GC content than do the flanking introns; and leveled GC, where there is no significant GC difference between the exon and introns. These groupings enabled us to examine DNA methylation on exon-intron architecture in an unbiased fashion and in a manner similar to our examination of nucleosome occupancy . We found DNA methylation to be a strong marker of exons regardless of GC content. Remarkably, nucleosomes are very selective in their GC-content requirements, marking exons with a differential GC content between exon and flanking introns but not exons with leveled GC content between exon and introns. In our analyses, we used mapped DNA-methylation data to calculate average DNA-methylation values: mCpG/CpG. Surprisingly, we found that when there is no GC differential between exon and introns, the mCpG/ CpG level drops significantly in the intronic regions close to the splice sites. We then extended our analysis to look at the full human exome through ''GC-content goggles'' in the form of isochore maps (Costantini et al. 2006) and found the same patterns of epigenetic preferences to be a global genomic property. Moreover, we found that methylated CpGs in a leveled GC architecture accompanied higher inclusion of exons. Overall, our findings provide further support for the role of DNA methylation as a splicing regulator and reveal dynamic patterns of DNA methylation and chromatin organization upon exon-intron structure that are dependent upon the regional GC content and the GC differential between exon and introns.

Results
The level of methylated CpGs as a marker of exons When attempting to understand the structure of DNA methylation around exons, one encounters a problem: GC content is usually higher in exons than in introns (Schwartz et al. 2009). This creates a difficulty when trying to assess relative values of DNA methylation in exons compared with flanking introns since DNA methylation requires CG dinucleotides, and those are more abundant when GC content is high. As a result, previous attempts that normalized DNA-methylation levels to cytosine abundance showed that the higher methylation marking of the exon dissolves and is similar to that in the flanking introns (Lister et al. 2009).
To assess whether higher DNA methylation in exons is due only to higher GC content in exons compared to the flanking sequences, we needed candidates that would enable us to simultaneously calculate DNA-methylation levels and control for GC content. We have recently identified two different approaches by which the splicing machinery is directed to exons and introns based on nucleotide composition of the genomic environment. In the first group, a differential GC content between exons and introns allows for better recognition of the exons by the splicing machinery. The second group is characterized by the same level of GC content in the exon and the flanking introns, and recognition by the splicing machinery is intron based ). The latter group, with a leveled GC architecture, provided us with the candidates with GC content that enabled unbiased observation of DNA-methylation patterns.
The analyses of both GC groups enabled a deeper understanding as to their splicing regulation differences. Specifically, a total of 15,874 exons with higher GC in exons compared with flanking introns (differential GC group) and 16,269 exons with leveled GC between exon and flanking introns (leveled GC group) were retrieved based on the RefSeq tracks from the UCSC Genome Browser (http://genome.ucsc.edu/) as were published previously by Amit et al. (2012). We calculated per base DNA-methylation levels for each group based on genome-wide DNA-methylation data from human H1 embryonic stem (H1 ES) cells retrieved from Lister et al. (2009).
When controlling for DNA-methylation biases, the CpG abundance also needs to be taken into account since higher methylation levels can be a direct result of a higher abundance of the CG dinucleotide. For example, the sequence ATTGGGGCAC has 60% GC content but 0% CG dinucleotides, whereas the sequence ACGCCAATCG also has a GC content of 60% GC content but has 40% CG dinucleutides. To control for this potential bias, we calculated the average DNA-methylation level only at CpG sites (mCpG/CpG). This calculation is based on the raw DNA-methylation level, which is the methylation base calls at a reference coordinate divided by total base calls at that coordinate as reported by Lister et al. (2009). Next, we calculate the mCpG/CpG value as the sum of DNA-methylation level for each position relative to the 39 or 59 splice site of each exon divided by the number of exons with a CpG occurrence in that position (see Methods). This calculation was necessary to avoid the bias of higher CpG density on the exons as we show further on. Both groups of exon-intron structures show a significantly higher level of mCpG/CpG on the exons compared with the flanking introns (t-test, P-value < 2.2 3 10 À16 ), and the overall level of methylated CpGs in the differential GC group is much higher than that of the leveled GC group ( Fig. 1, P-value < 2.2 3 10 À16 ). More specifically, in the group with differential GC content ( Fig. 1, gray line), there is an increase from ;80% methylated CpGs in the intronic sequences to ;85% in the exon. Strikingly, in the leveled GC group (black line), there is a drop in CpG methylation in the intronic regions proximal to the exon (;100 nt). If only the 100 intronic nucleotides flanking the exons are considered, there are ;55% intronic methylated CpGs and >70% exonic methylated CpGs (>30% increase). This indicates that there is significantly more DNA methylation in exons than in introns, and this difference is most apparent when exons and introns with similar GC content are compared. It was previously reported that a lesser fraction of CpGs are methylated in regions of high CG content (such as that of the leveled GC group) than in lower GC regions (differential GC group) ( Jabbari and Bernardi 1998;Oakes et al. 2007). Our results are in agreement with these findings.
In parallel with the analysis of mCpG/CpG, we also calculated the basic DNA methylation level (by base calls) and observed the expected stronger signal within exons relative to introns (Supplemental Fig. S1A). Another necessary analysis is that of the CpG abundance on the exon-intron structure since the human genome has a very low abundance of CG dinucleotides in general, but a much higher CG dinucleotides abundance in the coding sequence compared to intronic or intergenic regions (Karlin and Burge 1995;Gentles and Karlin 2001). In agreement with this, the results exhibit a strong signal of CpG abundance on the exons of both leveled and differential GC groups (Supplemental Fig. S1B). These results are a direct consequence of the CpG abundance bias in the coding sequence and prove the necessity in normalizing methylation levels to CpG abundance using the mCpG/CpG values.
We extended our analysis further to evaluate DNA-methylation levels in mouse ES cells and various human tissues based on data from Gu et al. (2010), who constructed genome-wide single-base DNA-methylation maps of the human and mouse samples using the reduced representation bisulfite sequencing (RRBS) protocol, which applies the DNA methylation-insensitive restriction enzyme MspI to map ;5% of the CpG sites in the genome. We found that DNA-methylation levels in the different human and mouse tissues exhibit the same general pattern of higher exonic values in both GC groups (Supplemental Fig. S2). The noise level is much higher in these data, however, since the methylation data set is much smaller.
In previous work , we found that when examining the differential GC group, there is a 50% increase in nucleosome occupancy on exons compared to flanking introns. However, there is a very small increase of the nucleosome occupancy signal (;10% increase) on the exons belonging to the leveled GC group. Although the two GC-content groups differ in the occupancy of nucleosomes between exons and flanking introns, both have a strong methylation signal on the exons. Moreover, the group that is weakly marked by nucleosomes is more significantly marked by methylated CpGs than the group marked more strongly by nucleosomes. This suggests that DNA methylation in exons has biological relevance. The difference in epigenetic patterns may indicate differences in how the splicing machinery recognizes these two types of exons.

DNA methylation and exon recognition
To evaluate the role of DNA methylation in exon recognition, we examined whether the level of mCpG/CpG in alternatively spliced cassette exons differs from that of constitutively spliced exons. We performed an extensive analysis to identify alternative and constitutive exons using RNA-seq data obtained from Lister et al. (2009). The mRNA reads were extracted from the same cells used to obtain DNA methylation (H1 ES cells), allowing the integration of alternative splicing effects with DNA-methylation patterns. Next, we used the SpliceTrap tool  to quantify exon inclusion ratios based on the RNA-seq data (see Methods).
We identified 37,473 exons with canonical splice sites that we could statistically ensure their inclusion levels. We then divided these exons into two groups based on GC content: (1) leveled GC, which had no significant difference in GC content between exon and flanking introns; and (2) differential GC, in which exons are significantly higher in GC than flanking introns (P-value < 0.05). There were 7413 exons in the differential GC group (5734 constitutive exons and 1679 alternative exons) and 6037 exons in the leveled GC group (4936 constitutive exons and 1101 alternative exons). GC-content maps of both groups confirmed that the expected exon-intron GC structure was present in each group ( Fig. 2A). The general pattern of mCpG/CpG signal around the exons remains the same for the differential GC constitutive exons and the leveled GC constitutive exons (multiple t-tests, P-values < 2.2 3 10 À16 ) (Fig. 2B,C). Importantly, the levels of mCpG/CpG are significantly lower in alternative exons than in constitutive ones (multiple t-tests, P-value < 2.2 3 10 À16 ) in both GC groups (Fig. 2B,C). These results further Average of methylated CpGs along exon-intron structure with a differential GC content between the intron and the exon (gray) and along exon-intron structure in which the GC content is identical between the exon and the flanking introns (black). The average value was calculated per base for exons (75 nt from each splice site) and flanking intronic regions (200 nt). A running average of 20 was applied after omitting the following positions for having no CpG occurrences: 39ss positions À4 to À1 and 59ss positions +1 and +2.
imply the conditions required for DNA methylation to play a part in the regulation of alternatively spliced exons. DNA-methylation levels by base calls, as well as CpG abundance were also calculated for these alternative and constitutive exons (Supplemental Fig. S3). Both present stronger signals on the exons than on the introns, yet the patterns are similar between alternative and constitutive exons of each group. These results emphasize, again, the need to observe the unbiased pattern of DNA methylation through the use of mCpG/ CpG. When this value is used, differential DNA methylation still clearly defines the exons compared to their flanking intronic regions and also better marks constitutive exons compared to alternative ones. Analysis of nucleosome occupancy was also done for these alternative and constitutive exons using the data used for the original GC-content groups ). We found a strong signal of nucleosome occupancy for differential GC exons and a very weak signal for leveled GC exons (Supplemental Fig. S4).

CpG abundance and chromatin organization
When examining methylation levels in regions immediately adjacent to exons, we observed extremely high methylation values near both the 39 and the 59 splice sites. This may be due to the presence of the CG dinucleotide as part of the consensus splice signal. To determine the relevance of this methylation, we first focused on deciphering the specific patterns of CpG abundance at positions À20 to +20 nt relative to each splice site. As a control, we used a data set of approximately 130,000 pseudo exons constructed by Ke et al. (2011). The GC content of these pseudo exons is not biased to any GC differential between exon and introns (Supplemental Fig. S5). In Figure 3A, the CpG abundance at the splice sites without application of a running average is shown. First, as previously observed, we found that CpG abundance is higher in the exons than in flanking introns in both GC groups but not in the group of pseudo exons. Second, we identified three positions with a significantly higher CpG percentage than the rest of the evaluated positions. All three are part of the splice site signal: position À5 of the 39 splice site, position +4 of the 59 splice site, and position À2 of the 59 splice site. When followed by G, the cytosines in all three positions are methylated in at least 70% of the cases. These same positions show very minor increase (if at all) in the group of pseudo exons. Interestingly, a CG dinucleotide is not part of the consensus splice site signal in these positions. However, the nucleotide adjacent to two of the positions (59 splice site À2 and +4) in the consensus sequence is a G, which partially accounts for the high CpG percentage in those two positions. These peaks in CpG abundance that are mostly methylated demonstrate very high methylation levels within both splice site regions and might point to a regulatory role for those positions when methylated. The fact that these peaks are diminished in the group of pseudo exons further implies on a possible role for these positions in the splicing process.
We next constructed several analyses to examine the role of these peak methylated positions as possible chromatin remodelers. It has been shown that certain DNA sequences with high affinity binding to the histone can direct nucleosome positioning (Lowary and Widom 1998), and that this in turn acts as a barrier to transcription (Bondarenko et al. 2006;Bintu et al. 2012). Furthermore, DNA methylation may have an intrinsic effect on nucleosome positioning on the DNA (Chodavarapu et al. 2010;Cedar and Bergman 2012). Thus, a change in nucleosome occupancy that accompanies DNA methylation might have an indirect influence on Pol II elongation rate and cotranscriptional splicing. We are not aware of any method that enables a quantification of the effect of a single base change in CpG-methylation upon chromatin organization. Thus, we decided to perform an analysis that was independent of the previous methylation analyses that relied on data from H1 ES cells and examined nucleosome occupancy as a result of CpG dinucleotide compositions at several positions relative to the splice sites using the data compiled by Schones et al. (2008). In H1 ES cells, a mean of 82.3% of the CpG sites are methylated (Lister et al. 2009), and 70%-88% of the CpG dinucleotides in both GC groups are methylated (see Fig. 1); therefore, we assumed that a CpG dinucleotide could act as a representative of a methylated position. Our next step was to divide the exons into subgroups based on dinucleotide composition at the three methylation peak positions: (1) CG dinucleotide at the peak position (CG composition subgroup); (2) all other dinucleotide with C and G that is not CG (GC composition subgroup); (3) all other dinucleotide compositions (i.e., AA, AT, AC, not-C composition subgroup); and (4) AG dinucleotide (AG subgroup). The GC composition subgroup was used as a control for the CG composition since the GC content is identical but CpGmethylation is not possible. For the analysis of the À2 position of the 59 splice site, we used the AG subgroup, which represents the consensus dinucleotide at position À2 of the 59ss (this subgroup is contained in the larger not-C subgroup). We constructed nucleosome occupancy data sets for the four exon subgroups based on their composition at the three peak positions for both leveled and differential GC exons.
There was a substantial increase in nucleosome occupancy when any of the three peak positions is a CG in the group of leveled GC exons (Fig. 3E-G). More specifically, we observe significantly higher levels of occupancy in the vicinity of all three positions when the composition is CG compared with any other composition examined (multiple t-tests, P-value < 1.76 3 10 À13 ). This effect is diminished in the group of differential GC exons (Fig. 3B-D), where we observe a depleted signal on the exons when position À5 of the 39 splice site is a CG (Fig. 3B, blue line) rather than GC composition (P-value < 0.006) (Fig. 3B, green line) and a small increase in signal for positions À2 and +4 of the 59 splice site (P-value < 7.8 3 10 À7 and P-value < 0.043, respectively) ( Fig. 3C-D). We extended this analysis to several other positions along the exon and flanking introns. In general, we observe a trend of increased nucleosome occupancy when a CG is present for the group of leveled GC exons, but the effect is smaller than at the peak positions (Supplemental Fig. S6). This effect was not observed in differential GC exons (Supplemental Fig. S7); in these exons, nucleosome occupancy is not significantly affected by dinucleotide composition near the splice site. We also validated that a CG at the three peak positions does not have a fundamental effect on the overall GC content of the subgroup. The GC content analysis of exons with the four dinucleotide compositions exhibit very minor changes, if at all, between the different composition subgroups (Supplemental Fig. S8) that cannot explain the strong nucleosome occupancy signal of the CG composition subgroup that is observed in Figure 3.
Two observations are apparent from these analyses: First, there is a significant change in the occupancy of nucleosomes in leveled GC exons when CG dinucleotides are present in the region of the splice sites; this difference is significant relative to any other composition at positions À5 of the 39 splice site and positions +4 and À2 of the 59 splice site and relative to the AG consensus sequence composition at position À2 of the 59ss. Since most CG dinucleotides at these positions are methylated, we assume that this may help to direct nucleosome binding to exons that can, in turn, be mediated by methyl-binding proteins. Second, these results again highlight differences between the leveled and differential GC groups with regards to the effect of DNA methylation. In the differential GC group, there is little if any effect of CG dinucleotides near the exon-intron junctions. In contrast, in the leveled GC group the effect is a strong nucleosomal signal. This difference suggests that the mechanisms of exon recognition may differ depending on the GC-content environment of a particular exon, the whole gene, or the genomic location itself, as we show further on.

Epigenetic patterns on exons based on genomic location
The major differences in epigenetic patterns between the leveled and differential GC groups led us to consider that perhaps genomic location impacted epigenetic pattern since some chromosomes are known to contain higher GC content than others (Costantini et al. 2006). Thus, if the groups belong in different chromosomes, it could explain the difference in basic GC-content architecture. It is also possible that the genomic location results in a particular epigenetic pattern and is another factor in exon recognition.
Analysis of the distribution of each GC group (leveled or differential GC) along the chromosomes shows that the groups are not equally distributed along the genome. We observed that the exons from the differential GC group (Fig. 4, light gray) are more abundant in lower GCcontent chromosomes, whereas the leveled GC exons (dark gray) are more abundant in high GC-content chromosomes. This being the case, the patterns of epigenetic modifications that are observed along the exon-intron structure (Figs. 1, 2) might be  a general property that is dependent on regional GC content. To evaluate this, we used human isochore maps that were constructed by Costantini et al. (2006). These maps divide the entire human genome into regions no smaller than 300K bases that are largely homogeneous in GC content. Costantini et al. (2006) divided those regions into two gene spaces that they called ''core'' isochores, where GC content is >46%, and ''desert'' isochores, where GC content is <46%. We divided all available internal exons into core and desert groups and obtained 53,102 core exons and 85,979 desert exons. As expected based on our previous results, DNA methylation is a strong marker of an exon when these ;140,000 human exons are considered (Supplemental Fig. S9A). The methylated CpGs were constructed for alternative and constitutive exons that are expressed sequence tag (EST)-based since the RNA-seq data provided reliable inclusion level for ;37,000 exons, which was insufficient for this analysis. We observe a general marking of the exons by a higher level of methylated CpGs (Supplemental Fig. S9B) that also exhibits a slightly larger mean increase in core exons (+24%) than in desert exons (+19%). Also, the level of methylated CpGs drops in alternative exons and their flanking introns in both desert and core exons relative to constitutively spliced exons as was observed using RNA-seq-based groupings. The differences in patterns and levels of methylated CpGs when exons are grouped using EST-based data rather than RNA-seq data are a consequence of the more heterogeneous GC differential of the isochore exons as they are grouped based on general GC content and not based on exonintron GC differential.
Analysis of nucleosome occupancy revealed that nucleosomes are more often present on exons rather than introns when regional GC content is low (Supplemental Fig. S9C, green lines). Thus, when examining sequences that reside in low GC regions, one can expect a strong occupancy of nucleosomes on the exons. However, when GC content is regionally high, nucleosomes are spread more evenly along introns and exons (Supplemental Fig. S9C, blue lines).
In these results, we observe two genome-scale exon populations that vary based on their GC content. We find them to differ in their epigenetic patterns upon exon-intron structure, and this difference could mean a different mechanism of splicing regulation that is based on genomic location.

Discussion
In this study, we analyzed DNA methylation at single-base resolution and evaluated nucleosome occupancy patterns along exonintron structures. To eliminate biases generated by differential GC content between exon and introns, we analyzed two groups of exons with and without a GC-content differential between exon and flanking introns. This enabled us to examine the pattern of DNA methylation across the exon-intron structure, taking into account, for the first time, GC-content differential. We previously showed that the GC differential between exon and intron allows better recognition of exons by the splicing machinery; yet it was unclear how exons are recognized in the leveled GC architecture, with similar GC content in exon and introns. The use of these groups for analysis of the effect of DNA methylation on splicing has two advantages: (1) it makes a control possible for a GC content bias between exon and introns (that is evident when analyzing DNA methylation); and (2) it reveals epigenetic differences  between the groups that might explain differences in exon recognition mechanisms.
Division of exons into groups characterized by GC differentials between exons and flanking introns provides a tool to analyze epigenetic modifications that are generally influenced by GC content (Bernardi et al. 1985;Jabbari et al. 1997;Segal et al. 2006;Varriale and Bernardi 2010). However, CpG density in itself could determine DNA-methylation levels (Choi et al. 2009). In order to view DNA methylation while taking into account CpG content, we compared the number of methylated CpGs to total CpGs, thus creating the mCpG/CpG ratio. We previously showed that DNA methylation levels can help predict alternative cassette exons ); a closer inspection of our current results reveals that it is not the absolute methylation value that distinguishes constitutive exons from alternative ones but the differential in the ratio of mCpG/CpG between exon and introns, which is also dependent in the general GC environment. We find a significantly higher level of methylated CpGs in exons compared with introns, regardless of GC-content differential. However, differential mCpG/CpG is much higher (;50%) in the leveled GC exon-intron architecture, where DNA methylation drops significantly in the intronic regions proximal to the exon. Remarkably, in this group there is no well-defined nucleosome marking the location of the exon relative to the flanking introns; nucleosomes mark exons significantly only when there is a substantial difference in GC content between exon and introns . Also, we previously showed that the exons exhibiting differential GC content and higher nucleosome occupancy compared to flanking introns are better recognized by the splicing machinery. Thus, CpG methylation can presumably be a part of the code that allows the splicing machinery to locate exons that have no GC differential. A possible explanation may be that the drop in DNA-methylation levels in proximal regions to the exon allows binding of certain proteins that can affect exon inclusion, such as is the case with transcriptional repressor CTCF protein (Shukla et al. 2011). On the other hand, by dividing exons based on RNA-seq data from H1 ES cells into alternative and constitutive, we identified a strong decrease in methylated CpGs in alternative exons as well as their flanking sequences compared with constitutive exons. This suggests that lower levels of DNA methylation of the whole intron-exon-intron strip are associated with suboptimal recognition of alternatively spliced exons.
The effect on splicing of a strong DNA-methylation signal in a leveled GC architecture may be indirectly through Pol II stalling. Since DNA methylation is a chromatin remodeler (Sarraf and Stancheva 2004;Klose and Bird 2006;Choy et al. 2010) and nucleosome occupancy causes Pol II stalling (Batsché et al. 2006;Hodges et al. 2009), these two mechanisms could be interlaced. Thus, it may be that peak CpG abundance at the splice sites might influence nucleosomal positioning. Through the use of the dinucleotide composition analyses, we found that in a leveled GC architecture nucleosome occupancy is strongly affected by the presence of CG dinucleotides. We consider the differences between the two GC groups as an indicator of different epigenetic interplays around exons in different GC environments.
It was previously shown that DNA sequence with high histone-binding affinity can direct nucleosome positioning and create a barrier to Pol II (Lowary and Widom 1998;Bondarenko et al. 2006). A recent paper by Bintu et al. (2012) suggested a mechanism by which there is an increase in pause density of Pol II in the central region of the nucleosomes (H3/H4 tetramer center) that is lessened in the entry and exit regions (first and last H2A/H2B dimmers, respectively). Taken together with the observed effect of CG composition on nucleosome occupancy, we propose that perhaps these findings point to a similar structure upon the exons that is governed by CpG methylation, where the peaks in CpG methylation at the 59 and 39 splice sites act as the central area binding for the nucleosomal barrier, whereas the drop in methylated CpGs in intronic flanking regions of leveled GC exons point to the entry/exit regions of the nucleosome. However, further investigation is required to better understand and tackle this hypothesis.
In conclusion, our findings provide an extensive, genomewide overview of DNA methylation and nucleosome occupancy relative to exon-intron structures. The analyses here support a role for DNA methylation in splicing and in chromatin remodeling. Our data suggest that the regional location of exons within the genome impacts their GC environment and that this, in turn, has a significant effect on how the splicesomal machinery recognizes the exons. This should be taken into account when examining the effects of other epigenetic markers, such as histone modifications, upon exon-intron structures.

Construction of EST-based exon and intron data sets
Data sets of human exons for the differential GC group and the leveled GC group were retrieved based on the RefSeq tracks from the UCSC Genome Browser (http://genome.ucsc.edu/) as described previously by Amit et al. (2012). The differential GC group contained 15,874 exons and flanking introns, and the leveled GC group contained 16,269 exons and flanking introns. GC-contentbased exonic groups for the mouse genome were also retrieved as described by Amit et al. (2012). The mouse differential GC group included 15,758 exons and flanking introns, and the mouse leveled GC group contained 16,583 exons and flanking introns. All exons examined are internal exons that are flanked on both sides by introns.
Data sets of human exons for the ''desert'' and ''core'' isochore groups were retrieved based on the RefSeq tracks from the UCSC Genome Browser (http://genome.ucsc.edu/). All internal exons of the human exome were divided into two groups based on the GC content of the isochore. Isochore maps were retrieved from Costantini et al. (2006). Exons that reside in regions of >46% GC content are considered core exons, and exons that reside in lower GC regions are considered desert exons. Overall, 85,979 exons and flanking introns construct the ''desert'' isochore group, of which 9649 are alternative exons; and 53,102 exons and flanking introns construct the ''core'' isochore group, of which 6443 are alternative exons.
Construction of RNA-seq-based leveled and differential GC exonic groups For construction of constitutive and alternative exon data sets, we used RNA-seq raw reads that were obtained from Lister et al.( 2009). The mRNA-seq reads were aligned to the human genome (hg18) using the Bowtie alignment tool (Langmead et al. 2009). Next, we applied the SpliceTrap software  to quantify exon inclusion ratios from the mapped RNA-seq reads. SpliceTrap is a statistical tool built to quantify exon inclusion ratios from RNAseq data. SpliceTrap quantifies for every exon the extent to which it is included, skipped, or subjected to size variations due to alternative 39/59 splice sites or Intron Retention. We identified 37,473 exons with canonical splice sites for which inclusion levels could be measured reliably with a minimum number of three reads per junction for each exon isoform. Next, we divided the identified exons into leveled and differential GC exons. For the differential GC group, we searched for exons with significantly higher GC content (t-test, P-value < 0.05) than their flanking introns, and for the leveled GC group we searched for exons with no significant differences in GC content between exons and flanking introns. We identified a total of 7413 exons in the differential GC group (5734 constitutive exons and 1679 alternative exons) and 6037 exons in the leveled GC group (4936 constitutive exons and 1101 alternative exons).

Construction of basic DNA methylation maps by base calls
Genome-wide single-base resolution data on DNA methylation was obtained from Lister et al. (2009), who conducted MethylCseq analyses to map genome-wide DNA methylation for two human cell lines, H1 human embryonic stem cells and IMR90 fetal lung fibroblasts. In this method, they used bisulfite-conversion to convert cytosine to uracil (which is later transformed by PCR amplification to thymidine), while methylated cytosines are not affected. Next, high-throughput sequencing was processed using the Illumina analysis pipeline. The Bowtie alignment tool (Langmead et al. 2009) was used to align the results to the human reference genome (hg18). The base calls per reference position on each strand were used to identify methylated cytosines at 1% FDR.
To map the methylated cytosines, we cast single-base resolution values upon sets of genomic intervals (in this case, exons and introns), resulting in three matrices for each set of exons: upstream intron, exon, and downstream intron. Next, we construct two databases for the two splice sites. One database we centered on the 39 slice site (intron positions are below zero and exon positions are above zero). The other database was centered on the 59 splice site. All sequences were aligned to the plus strand; thus, methylation values of an exon that is coded on the minus strand are reversed. Once all exons and introns were positioned relative to the splice site, we calculated per base statistical measures of methylation values for all intervals in a genomic set. The exon-intron strip was constructed from 75-nt exonic nucleotides and 200 intronic nucleotides. Exons shorter than 75 nt received a null value for the remaining positions, and the same is done for introns shorter than 200 nt. For exons that are shorter than 150 nt, there is a duplication of the data at 59 and 39 splice sites. For example, for an exon with the length of 100 nt, positions 1-75 will be used for display at the 39 end of the exon, and positions 25-100 will be used for display at the 59 end of the exon. We used this method in order to be faithful to the true values with regard to each splice site and not to the middle of the exon. All scripts were written using the perl script language and the R statistical computing program (Team 2006), the latter was also used for both statistical analysis and graphical display.
DNA-methylation maps for the different human tissues and mouse ES cells were constructed using data retrieved from Gu et al. (2010), who performed genome-wide single-base DNA-methylation mapping of clinical samples using the RRBS protocol. The tissues examined were primary colon tumor, normal colon, blood cells from a healthy individual, and blood cells from an individual with colon cancer.
Calculation of the mCpG/CpG score Calculation of the value of mCpG/CpG was done to take into account the original base calls of DNA methylation measured by Lister et al. (2009) and the number of CpGs in a certain position relevant to the exon-intron junction of all measured exons. This calculation is based on the base call methylation data sets of exons and introns that are positioned relative to the splice sites. For each nucleotide relative to the splice site, we extracted only the exons/ introns that have CG dinucleotide in that position. Next, we calculated the average methylation level only in those cases. For example, in a data set of 20 exons, 10 exons in the differential GC group have a CpG at position À1 of the 59 splice site. If methylation values are 0.8*4, 0.7*3, and 0.9*3, the mCpG/CpG value will be (0.8*4 + 0.7*3 + 0.9*3)/10 = 0.8. The same calculation in the leveled GC group will include more exons with methylation values (15/20) and more CpG positions (15). This method takes into account cases in which a certain exonic position is a CpG but the methylation base calls are zero (numerator can be zero). However, the method does not take into account cases in which the denominator is zero (there is no CpG at examined position of an examined exon).

CpG abundance around splice sites
The CpG abundance is the percentage of exons that have a CpG at an examined position. For example, in a data set of 100 exons, if two exons have a CpG at position +1 of the 39 splice site and three exons have a CpG at position +2 of the 39 splice site, the CpG abundances will be 0.02 and 0.03, respectively. CpG abundance was also calculated on a data set of 134,935 pseudo exons that was obtained from Ke et al. (2011). The pseudo exons were defined as intronic sequences having lengths between 50 and 250 nt and a splice site score of above 75 for 39 splice site and above 78 for 59 splice site, using the position specific scoring matrix (PSSM) method of Shapiro and Senapathy (1987). The consensus values are in the range of 0 to 100; the median value for a ''real'' exon 39 splice site is 80, and that for a ''real'' exon 59 splice site is 82. Pseudo exons that contain interspersed repeats were removed using RepeatMasker.

Construction of nucleosome occupancy maps for exons
We obtained nucleosome-occupancy levels using nucleosome score profiles generated by Schones et al. (2008). The maps were constructed using Illumina high-throughput sequencing of DNA fragments attached to nucleosomes in activated T cells, following micrococcal nuclease (MNase) digestion (Schones et al. 2008). In order to map nucleosome positioning upon the exon-intron genomic regions, we cast single-base values upon genomic intervals as was done for DNA-methylation data. Next, we applied a running average of 20 nt on the vector of average values for the exon-intron strip.
Nucleosome occupancy maps based on positional nucleotide composition were constructed using the same methods. We constructed nucleosome occupancy maps for exons based on their dinucleotide composition in 20 different positions along the exon-intron sequence. For the 39 splice site, we constructed nucleosome occupancy data sets at positions À150, À20, À10, À6, À5, +1, +2, +3, +10, and +20 relative to the exon-intron junction. For the 59 splice site, we constructed nucleosome occupancy data sets at positions À20, À10, À5, À3, À2, +4, +5, +10, +20, and +150 relative to the exon-intron junction. For each position, we constructed sub-data sets based on the following dinucleotide composition: (1) CG, all exons with a C at the position examined and G at the next position; (2) other C and G combinations: CC, GG, and GC; exons where the dinucleotide was either preceded by a CG dinucleotide (as in ''CGC'' or ''CGG,'' i.e., DGC\G) or followed by a CG dinucleotide (as in ''CCG'' or ''GCG,'' i.e., C\GCH) were excluded to prevent a case in which a methylated CpG in the vicinity to the position at hand will affect the result; (3) all other possible dinucleotide compositions at the requested and consecutive po-sitions (AA, AT, AG, AC, etc.); and (4) AG, all exons with a C at the position examined and G at the next position; this subset of group 3 was used only for evaluation of position À2 of the 59 splice site since it represents the consensus sequence at that position. We then constructed average nucleosome occupancy maps for each subgroup of exons separately. Since there is no necessity for inclusion level liability in this analysis yet there is a strong statistical need for large subgroups, we analyzed the larger exonic data sets that are based on the work by Amit et al. (2012) and not the smaller groups that originate from RNA-seq experiments.

Pictograms
Graphical representations of PSSMs were composed using a BioPerl (Stajich et al. 2002) module for generating Scalable Vector Graphics (SVG) output of Pictogram display for consensus motifs, as was described by Burge et al. (1999). The height of each letter is proportional to the frequency of the corresponding base at the given position, and bases are listed in descending order of frequency from top to bottom. Pictograms were constructed for the 59 splice site of differential and leveled GC exons and also for a subgroup of exons with a CG dinucleotide composition at position À2 of the 59 splice site.