A Study of Alterations in DNA Epigenetic Modifications (5mC and 5hmC) and Gene Expression Influenced by Simulated Microgravity in Human Lymphoblastoid Cells

Cells alter their gene expression in response to exposure to various environmental changes. Epigenetic mechanisms such as DNA methylation are believed to regulate the alterations in gene expression patterns. In vitro and in vivo studies have documented changes in cellular proliferation, cytoskeletal remodeling, signal transduction, bone mineralization and immune deficiency under the influence of microgravity conditions experienced in space. However microgravity induced changes in the epigenome have not been well characterized. In this study we have used Next-generation Sequencing (NGS) to profile ground-based “simulated” microgravity induced changes on DNA methylation (5-methylcytosine or 5mC), hydroxymethylation (5-hydroxymethylcytosine or 5hmC), and simultaneous gene expression in cultured human lymphoblastoid cells. Our results indicate that simulated microgravity induced alterations in the methylome (~60% of the differentially methylated regions or DMRs are hypomethylated and ~92% of the differentially hydroxymethylated regions or DHMRs are hyperhydroxymethylated). Simulated microgravity also induced differential expression in 370 transcripts that were associated with crucial biological processes such as oxidative stress response, carbohydrate metabolism and regulation of transcription. While we were not able to obtain any global trend correlating the changes of methylation/ hydroxylation with gene expression, we have been able to profile the simulated microgravity induced changes of 5mC over some of the differentially expressed genes that includes five genes undergoing differential methylation over their promoters and twenty five genes undergoing differential methylation over their gene-bodies. To the best of our knowledge, this is the first NGS-based study to profile epigenomic patterns induced by short time exposure of simulated microgravity and we believe that our findings can be a valuable resource for future explorations.


Introduction
During space flight, astronauts are exposed to powerful environmental assaults such as microgravity, cosmic radiation and magnetic fields that have the potential to impinge upon cellular ontogeny through epigenetic modifications [1]. Throughout the evolutionary history, gravity has been a constant factor in defining the architecture and morphology of living beings [2]. Hence a broader understanding of gravity's influence on biological functions is important for an accurate evaluation of risks associated with the health of astronauts in spaceflights and should be of enormous interest to the scientific community. The effects of microgravity in altering gene expression have been documented in mammalian cells [3,4] and other model organisms, such as yeast and bacteria [5][6][7]. Microgravity associated pathological alterations include reduction in bone mass and calcium concentrations [8], alterations in hormonal levels [9], impairment of immunocompetence [10] and apoptosis signaling [11]. Studies of human lymphoblast and lymphoblast cell cultures following periods of simulated microgravity have demonstrated alterations in metabolic processes and DNA repair pathways which could in turn signify an increased susceptibility to malignancy [12,13]. Collectively, these studies indicate exposure to microgravity during space flight alters gene expression patterns and subsequently cellular physiology.
DNA methylation is regarded as a major epigenetic mechanism and play key roles in regulating cellular processes in living organisms [14,15]. Biochemically, DNA methylation refers to the addition of a methyl group (CH 3 ) to the 5' carbon on the pyrimidine ring of cytosine nucleotides (commonly abbreviated as 5mC). Aberrations in genome-wide 5mC patterns are widely prevalent in cancer and other diseases [14,[16][17][18]. Traditionally DNA methylation marks have been associated with "transcriptionally silent" genes, however the revelations of global methylation studies facilitated by recent advances in Next Generation Sequencing (NGS) tools have established that the role of 5mC in regulating gene expression is complex, varies according to the genomic context and warrants extensive explorations [19][20][21][22][23][24][25]. Discovered in 2009, DNA hydroxymethylation (5hmC) is a relatively new epigenetic modification occurring on Cytosine [26,27] generated by Ten-Eleven Translocation (TET) protein-mediated oxidative catalysis of 5mC [26]. Though, potential roles of 5hmC at promoter and gene bodies are not clearly understood, it is shown to play some role in maintaining and/or promoting gene expression [14,[16][17][18]28]. Microgravity induced alteration in DNA methylation patterns have been reported previously [29][30][31] but the effect of microgravity on 5hmC is virtually unknown. During the time period this study was being conducted there were no reports of a NGS based study documenting the effects of microgravity on the epigenomic landscape.
The goal of our study was to profile genome-wide effects of "simulated" microgravity on 5mC, 5hmC and gene expression patterns employing Next Generation Sequencing (Fig 1). The TK6 lymphoblastoid cell line, derived from T cell blast crisis of a patient with chronic myelogeneous leukemia [32], served as our model organism. TK6 cells are well characterized and have been extensively used as a substitute for peripheral blood lymphocytes for immunological and epidemiological studies [13]. The limited availability of biological specimen subjected to conditions of microgravity in spaceflights makes ground based "simulated" microgravity studies critical in determining thresholds and thorough testing of the model organism before conducting the experiments during space missions [33]. A High Aspect Ratio Vessel (HARV) based rotary cell culture system (initially developed by NASA) was used in our study to "simulate" microgravity in the TK6 cells as has been described previously [34] and compared to a control static cell-culture system under the influence of earth's gravity. While assessing the merits of ground based "simulation" studies, it has to be appreciated that the effects of gravity cannot be completely negated but reduced to near zero to achieve a state of "functional near weightlessness" [33].

Materials and Methods
Cell culture TK6 human lymphoblastoid cells (ATCC, Manassas, VA) were maintained in the log phase of cell growth by culturing in RPMI-1640 (Life Technologies, Grand Island, NY) medium supplemented with 10% Fetal Bovine Serum (Atlanta Biologicals, Flowery Branch, GA) and 1% Penicillin/Streptomycin (Life Technologies, Grand Island, NY) at 37°C in 5% CO 2 and 95% air. For ground-based simulation of microgravity, HARV Rotary Cell Culture System (Synthecon, Houston, TX) was used. Actively growing TK6 cells were seeded in the bioreactor at 2 X 10 5 cells/ml and rotated at 12 rpm/min. In parallel, cells (at the same cellular density i.e. 2 X 10 5 cells/ml) were maintained in bioreactors in normal gravity (static) condition as controls. The bioreactors were maintained in an incubator at 37°C, with 5% CO 2 and 95% air for 48 hours.

DNA isolation, sonication and adapter ligation
Genomic DNA was isolated from the TK6 cells cultured under microgravity and control static conditions using the DNeasy Blood &Tissue kit (Qiagen Inc., Valencia, CA) following manufacturer's instructions. 2.5 μg of genomic DNA from each sample was sheared using Covaris S2 Device (Covaris Inc., Woburn, MA). Sheared DNA was purified by binding to AmPure beads (Beckman Coulter Inc.) and End-repair performed by incubating sonicated DNA and End repair solution (New England Biolabs Inc., Ipswich, MA) as per manufacturer's specifications. A-tailing was obtained by incubating the end repaired DNA with dA-tailing mix (New England Biolabs Inc., Ipswich, MA) at 37°C for 30 minutes. At this stage, to facilitate multiplexing each sample was equally divided in two parts (one half for MeDIP and the other half for hMeDIP respectively). Blunt end ligation was performed by incubating the A-tailed DNA samples (1 μg) with unmethylated versions of adapters (IDT Inc., Coralville, IA) based on sequences of the methylated Truseq adapters (Illumina Inc., San Diego, CA) for multiplexing. Thus 8 libraries were prepared as described above. Samples were assayed by qPCR in duplicate and standard curve constructed to establish the molarities of the libraries.

MeDIP-seq /hMeDIP-seq
MeDIP and hMeDIP were performed using the methylated/hydroxymethylated DNA enrichment kits (Diagenode Inc., Denville, NJ) following the manufacturer's protocol. Briefly, to 1.2 μg of adapter ligated sonicated genomic DNA, three DNA controls (known sequences bearing unmethylated, methylated or hydroxymethylated Cytosines respectively to assess the efficiency of immunoprecipitation reactions) were spiked-in. The concentration of genomic DNA was adjusted to incorporate the addition of the adapter sequences, preserving the appropriate molar ratio between the genomic DNA and anti-5mC/anti-5hmC antibody during MeDIP/ hMeDIP as described by Butcher  conc.gDNA ➔ adjusted genomic DNA concentration in the adapter ligated libraries, conc.adapter ligated gDNA ➔ concentration of the adapter ligated gDNA libraries, bpsonicated gDNA ➔ average size of the pre-ligation sonicated gDNA and bpadapter DNA ➔ average size of the adapters After incubation at 95°C to denature the double stranded DNA, immunoprecipitation was performed by incubation with monoclonal antibody directed against 5mC/5hmC (Diagenode Inc., Denville, NJ) and secondary antibody with magnetic bead conjugates (Diagenode Inc., Denville, NJ) overnight at 4°C while being spun continuously at 40 rpm. The captured 5mC/ 5hmC bearing DNA fragments were separated from the others by magnetic pulldown. After repeated cleanups, the captured DNA was isolated from the magnetic beads bearing antibody using the IPure kit (Diagenode Inc., Denville, NJ). The enrichment of 5mC/5hmC bearing DNA was assessed by performing qPCR on the pre and post immunoprecipitated samples. As a control, an identical immunoprecipitation reaction with mouse IgG instead of monoclonal 5mC/5hmC antibody was performed. The methylated/hydroxymethylated DNA immunoprecipitated libraries were amplified by PCR and submitted to the Purdue Genomics Core Facility for high-throughput sequencing by Hi Seq 2000 (Illumina Inc., San Diego, CA).

MeDIP-seq and hMeDIP-seq data processing
FastQC v 0.10.1 [36] was used to assess the quality of the reads and to generate graphical representations of numerous quality metrics (per base sequence quality, GC content and sequence duplication/size distribution levels). The reads were aligned to human reference genome hg19 using BWA v 0.6.2 [37], with default parameters and a maximum insert size of 400 bp. The resulting SAM files were converted to BAM format and sorted using Samtools v0.1.18 [38] as illustrated in (Fig 1). PERL script from the MeDUSA package [39] was used to convert the BAM files to BED format. Since the MEDIPS v1.0 [40] package requires only selective fields as input, the BED format was then reduced to four fields using the UNIX cut option. The MeDUSA pipeline utilizes the Bioconductor package MEDIPS v1.0 and custom R scripts to calculate quality metrics for the MeDIP-seq data were designed. The data was normalized for the size of the sequence libraries by calculating reads per million (RPM) in tiled windows across the genome. Wig files obtained for the normalized read depth following alignment and filtering were presented as RPM. Quality check on the MeDIP-seq data was also performed by calculating CpG enrichment values, saturation plots and coverage plots. Genome-wide correlations between the replicates were performed as a quality check for consistency among the replicates using QCSeqs from the Useq package (v8.40) [41] using a window size of 500 bp, increasing in 250 bp increments and a minimum number of 5 reads in a window.

Identification of DMRs and DHMRs
Peak calling software SPP v1.10 [42], was used to call peaks and rank them based on significance of enrichment (p-values and false discovery rates). IDR (Irreproducible Discovery Rate) framework was used to measure experiment quality in terms of reproducibility [43] and to select the reproducible, consistent peaks (overlapped significant peaks from both replicates) determined based on IDR values. The threshold of 0.05 IDR was used for truncating the peak list as suggested by the developers. The differentially methylated/hydroxymethylated regions identified by IDR analyses were then annotated with their chromosomal locations and feature types for further biological interpretation using custom Perl scripts of MeDUSA package, BED-Tools [44] and feature annotation files (GFF files from UCSC) as illustrated in Fig 1. Further annotation (plots for enrichment of 5mC/5hmC) was done using CEAS (v1.02) [45] and proximity of the peaks to the TSS was determined using PeakAnalyzer [46]. The FDR value of 0.05 was used as cut-off for all peak association studies. The complete MeDIP-seq and hmeDIP-seq data was submitted to NCBI GEO (GSE65944) and available in the database http://www.ncbi. nlm.nih.gov/geo/query/acc.cgi?acc=GSE65944).

RNA-seq
Total RNA was extracted from cells subjected to simulated microgravity and static control using RNA-STAT-60 (Tel-Test,Inc., Friendswood, TX) using the manufacturer's instructions. Briefly, 1ml of RNA-STAT solution was added per 10 6 cells and homogenized for 5 minutes over ice. 1ml of chloroform was added, contents shaken vigorously and centrifuged at 12,000g at 4°C for 15 minutes. The aqueous solution was transferred to corex tube (Corning Inc., Lowell, MA) and 0.8 ml isopropanol added. After incubation of 10 minutes, the contents were centrifuged at 12,000g for another 10 minutes to precipitate the RNA. The RNA pellet was washed with 75% ethanol and centrifuged at 7,500g for 5 minutes at 4°C. The ethanol was aspirated and the RNA pellet dried. The RNA pellet was finally resuspended in DEPC water and submitted to the Purdue Genomic Center for conversion into cDNA, sonication, adapter ligation and sequencing as described previously. The reads (fastq files) were aligned to human reference genome hg19 using Tophat v2.1.0 [47], with default parameters and known transcriptome as illustrated in Fig 1. Alignment results were filtered by Bamutils v0.5.0 [48] to remove reads with multiple mappings. Statistics data of the resulting alignment files were created using Samtools v0.1.18 [49] and Bamutils v0.5.0. The counts of aligned reads mapping to known genes were calculated using bamutils v0.5.0. EdgeR v2.11 [50] was used to compute the differentially expressed genes. Pathway analysis on the set of differentially expressed genes was done using the GeneCodis3 software designed at the Complutense University of Madrid [51].The complete RNA-seq data was submitted to NCBI GEO (GSE65944) and available in the database (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE65944).

Effect of simulated microgravity on cell growth and viability
The effect of simulated microgravity on cell growth and viability 48 hours after the cells were seeded in bioreactors in either the rotating or static condition was determined using Trypan Blue staining method by Automated Cell Counter (Nexcelom Bioscience LLC., Lawrence, MA) in S1 Fig. No significant differences in the percentage of viable cells between the two cell culture conditions after 48 hours was observed. Specifically, 95.1 ± 2.12% of the cells subjected to simulated microgravity were viable, while 91.1 ± 3.54% of the static control cells were determined to be viable. The average cellular diameter (μm) was determined to be 12.6 ± 0.42 and 12.8 ± 0.28 in TK6 cells subjected to microgravity and control respectively. Similar cellular growth rates between the rotating and static culture conditions facilitated ruling out the possibility of cell growth being the major contributor to the changes in the methylation and gene expression patterns.

Changes in 5mC profile following simulated microgravity
We applied MeDIP coupled with high-throughput sequencing to identify the differences in the genome-wide patterns of 5mC upon simulated microgravity on TK6 cells. 2.8x10 8 and 1.8x10 8 reads were obtained during MeDIP-seq from TK6 cells subjected to static and simulated microgravity respectively and more than 90% of these reads aligned to the human genome GRCh37/ hg19, 2009 Assembly (S1 Table and S2 Fig). Quality assessment generated by FastQC [36] showed satisfactory sequence quality for all measures except for GC content. As GC rich regions of the genome are enriched in MeDIP-seq datasets, this result was not unexpected. The depth of sequencing for MeDIP-seq samples ranged from 2.8X to 6.1X (S1 Table). Cross-correlation analysis was performed as per the ENCODE consortium guidelines [42,52,53] and all the samples displayed Normalized Strand Correlation (NSC) and Relative Strand Correlation (RSC) values (S1 Table) characteristic of "high-quality data sets". The similarity between the replicates was evident as hierarchical dendrogram displayed distinct clustering of biological replicates in two groups (S3 Fig) and sequence coverage analyses displayed that MeDIP-seq reads generated from the samples covered similar number of bases of the reference genome (S4 Fig).
Differentially methylated region (DMRs) were defined as genomic regions in TK6 cells under simulated microgravity that showed alteration in methylation (either increase or decrease) compared to TK6 cells under static conditions. 3204 DMRs (S2 Table & S3 Table) were detected using the IDR pipeline having an IDR cutoff value of 0.05 or less. Of the total DMRs, 1286 (40.14%) were associated with hypermethylation (gain-of-5mC) (S2 Table) and 1918 (59.86%) with hypomethylation (loss-of-5mC) (S3 Table) upon simulated microgravity respectively. The DMRs were further analyzed to determine the overlap of DMR regions with different genomic features by the methylome analysis pipeline described in details by Wilson et al. [39]. Functional genomic distribution analyses indicated that 969 and 1381 genes associated with DMRs have undergone gain-of-5mC and loss-of-5mC respectively ( Table 1). Also, 105 hypermethylated and 193 hypomethylated DMRs were observed around -1500 to 1500 bps of Transcription Start Sites (TSS) as demonstrated in Tables 2 & 3. The distribution of the genomic repeat sequences (LINE, SINE and LTR) located within the DMRs has been represented in Table 1. Metadata describing features such as genes, transcripts, Pseudogene, non-coding RNA and other regulatory features present on each DMR has been included in S2 Table & S3 Table. Investigation of annotations from 20 different ontologies from genomic coordinates of DMRs was generated by utilizing Stanford University's Genomic Regions Enrichment of Annotations Tool (GREAT) version 3.0.0 [54] and included in S4 Table & S5  Table. Gain-of-5mC DMRs induced by simulated microgravity were found to enrich GO Biological Processes like regulation of metabolic process (GO: 0019222), primary metabolic process (GO: 0044238) and cellular metabolic process (GO: 0044237) (Fig 2A). PANTHER Pathway Analysis implicated genes involved in p53 pathway (P00059), PI3 kinase pathway (P00048), T cell activation (P00053) and B cell activation (P00010) to be associated with hypermethylated DMRs (Fig 2B). On the other hand, loss-of-5mC DMRs were observed to enrich GO Biological Processes like cellular metabolic process (GO: 0044237) and primary metabolic process (GO: 0044238) ( Fig 2C). PANTHER Pathway Analysis, revealed that these hypomethylated DMRs were associated with genes involved in EGF receptor signaling (P00018), Apoptosis signaling (P00006) and FGF signaling (P00021) pathways among others ( Fig 2D).

Changes in 5hmC profile upon simulated-microgravity
We applied hMeDIP analyses coupled with high-throughput sequencing to identify the differences in the genome-wide patterns of 5hmC upon simulated microgravity on TK6 cells. 2.7x10 8 and 1.4x10 8 reads were obtained during hMeDIP-seq from TK6 cells under static and simulated microgravity respectively and more than 90% of these read uniquely aligned to the human genome GRCh37/hg19, 2009 Assembly (S1 Table and S2 Fig). The depth of sequencing for the hmeDIP-seq samples ranged from 1.8X to 4.6X depending on the sample (S1 Table). Cross-correlation analysis was performed as per the ENCODE consortium guidelines [42,52,53] and all the samples displayed Normalized Strand Correlation (NSC) and Relative Strand Correlation (RSC) values greater than the minimum threshold (S1 Table). The consistency of reads in the biological replicates were observed through the cluster analysis (S3  Table & S7 Table) generated at IDR < 0.05, 154 (92.2%) were associated with Table 1. Genome annotation Summary. The number of genomic features such as CpG islands, CpG shores, ENSEMBL Genes and DNA Repeats (LINE, SINE and LTR) associated with regions undergoing gain-of-5mC/5hmC and loss-of-5mC/5hmC DMRs or DHMRs in TK6 cells cultured under simulated microgravity compared to static condition.

Features DMR DHMR
Gain-of-5mC Loss-of-5mC Gain-of-5hmC Loss-of-5hmC  hyper-hydroxymethylation (gain-of-5hmC) (S6 Table) and 13 (7.8%) with hypo-hydroxymethylation (loss-of-5hmC) (S7 Table) upon simulated microgravity respectively. The overlap of DHMRs with different genomic features indicated that 86 and 7 genes were associated with gain-of-5hmC and loss-of-5hmC DHMRs respectively (Table 1). Also, 5 gain-of-5hmC ( Table 4) and2 loss-of-5hmC ( Table 5) DHMRs were observed around -1500 to 1500 bps of Transcription Start Sites (TSS). The distribution of DNA repeat regions present within the DHMRs was represented in Table 1. Metadata describing each DHMR was included in S6 Table & S7 Table. Investigation of GREAT version 3.0.0 ontology annotation [54] was included in S8 Table. Gain-of-5hmC DHMRs induced by simulated microgravity were found to be associated with genes that enriched in GO Biological Processes like positive regulation of B cell activation (GO: 0050871), positive regulation of B cell proliferation (GO: 0030890) and positive regulation of cell-cell adhesion (GO: 0034116) among others (Fig 3A). Panther Pathway Analysis of these gan-of-5hmC DHMRs implicated the muscarinic acetylcholine receptor signaling (P00042), insulin/IGF pathway-protein kinase B signaling cascade (P00033) and Fas signaling (P00020) among others (Fig 3B), Due to an extremely small gene set associated with loss-of-5hmC DHMRs, significant p-values of pathway associations were not obtained and have not been reported here.

Changes in the transcriptome upon simulated-microgravity
In TK6 cells, simulated microgravity induced differential expression of 370 transcripts out of 22,376 transcripts analyzed (FDR<0.1) compared to static control (S9 Table). 271 (73.24%) differentially expressed transcripts were associated with a decrease in expression, while 99 (26.76%) differentially expressed transcripts were associated with an increase in gene expression. 17 (4.59%) genes were associated with a drastic change of differentially expression (greater than 2 fold increase or decrease), while the vast majority were associated with intermediate (0-2 fold) change in differential expression. Furthermore, the pathway analysis (S10 Table) of transcriptionally upregulated genes showed enrichment of GO Biological Processes such as response to oxidative stress (GO:0006979) and ion transport (GO:0006811) (Fig 3C), while the downregulated genes could be linked to regulation of DNA-dependent transcription (GO:0006355) and carbohydrate metabolic processes (GO:0005975) (Fig 3D). Some of the top upregulated genes include CHAC1, TRPA1, ATAD3C, INHBE, CTH, HMOX, HBD, SPG20,       Table).

Correlation between simulated microgravity induced DMRS/DHMRs and gene expression
A comparison of the simulated microgravity induced differentially expressed genes (S9 Table) with DMRs located at gene promoters (Tables 2 & 3) revealed that two transcriptionally upregulated genes (TSPAN5 and SPG20) were associated with loss-of-5mC at their promoter and three transcriptionally downregulated genes (PLIN2, MAP3K13 and FBXO1) were associated with loss-of-5mC at their promoter. None of the gene promoters linked to DHMRs (Tables 4  & 5) were found to be differentially expressed (S9 Table). Similarly, the comparison of simulated microgravity induced differentially expressed genes with DMRs/DHMRs located at gene bodies revealed that 25 differentially expressed associated with DMRs at their gene bodies and none of the differentially expressed genes associated with DHMRs at their gene bodies. The relationship between methylation status at gene bodies and their respective transcriptional activity of these 25 differentially expressed genes did not show any significant correlation by Fisher's Exact Test (Fig 4A) and could be divided into five distinct groups, (i) five transcriptionally upregulated genes with loss-of-5mC DMRs at their gene bodies (CTH, CACNA1D, SPG20, PLS1 and SLC39A14), (ii) eleven transcriptionally downregulated genes with loss-of-5mC DMRs at their gene bodies (FBXO17, AKAP6, RIT1, GTF2IRD2P1, MSTO1, PMS2CL, MAP3K13, ST3GAL1, NCKIPSD, MAST1 and MSTO2P), (iii) three transcriptionally Effect of Simulated Microgravity on 5mC and 5hmC Patterns of TK6 Cells downregulated genes associated with gain-of-5mC DMRs at their gene bodies (CACNB2, WDR45B and CABLES1), (iv) three transcriptionally upregulated genes with gain-of-5mC DMRs at their gene bodies (CASZ1, VCL and ATF3) and (v) two transcriptionally upregulated genes with gain-of-5mC as well as loss-of-5mC DMRs at their gene bodies (ARID5B and TSPAN5) (Fig 4B). The comparison of DMRs and DHMRs located at gene bodies (S6 Fig) yielded six overlapping groups namely (i) 140 genes were associated with gain-of-5mC and loss-of-5mC DMRs at their gene bodies, (ii) eight were genes were associated with loss-of-5mC DMRs and gain-of-5hmC DHMRs, (iii) five gene were associated with gain-of-5mC and lossof-5mC DMRs as well as gain-of-5hmC DHMRs at their gene bodies, (iv) seven genes were associated with gain-of-5mC DMRs and gain-of-5hmC DHMRs at their gene bodies, (v) one gene was associated with gain-of-5mC DMR and loss-of-5hmC DHMR at its gene body and  (vi) two genes were associated with gain-of-5hmC DHMRs and loss-of-5hmC DHMRs at their gene bodies.

Discussion
The objective of this ground-based study was to map the genome-wide effects of simulated microgravity on DNA methylation, hydroxymethylation; and gene expression patterns in TK6 lymphoblastoid cells by a powerful Next Generation Sequencing pipeline. Although on the basis of numerous studies reporting microgravity-induced physiological changes in living organisms ranging from prokaryotes to humans, it has been speculated that microgravityinduced changes may occur in the methylome, very little is known about the effects of microgravity on DNA methylation. In 2009, Ou et al reported hypermethylation of a set of genes and transposable elements in rice (Oryza sativa L.) plants germinating from space-flown seeds [29].
Ou et al also reported that the spaceflight-induced hypermethylated genes did not generally correlate with alterations in their gene expression status [29]. In 2010, Singh et al reported that   Effect of Simulated Microgravity on 5mC and 5hmC Patterns of TK6 Cells consider: (i) differences in mechanisms establishing and maintaining DNA methylation patterns in plants and animals (for a comprehensive review refer to [55]), (ii) while Ou et al's investigation was based on spaceflight-induced "epigenetic memory" being transmitted from the seeds to the sapling, Singh et al had investigated the simulated microgravity-influenced changes in DNA methylation in immortalized T-lymphocyte cell cultures that might not be inheritable and (iii) while Singh et al had investigated the effects of only simulated microgravity on DNA methylation, Ou et al was investigating the effects of numerous factors like cosmic radiation, microgravity and space magnetic fields encountered during spaceflight. These reports therefore provided a strong basis for us to perform this study with advanced methods such as MeDIP-seq, hMeDIP-seq and RNA-seq to explore the relationship between the methylome and the transcriptome in microgravity exposed cells. To the best of our knowledge, this is the first report profiling the effects of simulated microgravity on the epigenomic landscape of human cells. 3204 DMRs and 2116 DHMRs distributed throughout the genome were identified in TK6 cells subjected to simulated microgravity. The majority of the DMRs (59.86%) were identified to undergo hypomethylation, which was consistent with the findings of Singh et al [30]. On the other hand the majority of DHMRs (92.2%) were associated with hyper hydroxymethylation. Additionally, we have been able to perform ontology based annotations to obtain information about the biological processes that might be affected by genes associated with simulated microgravity induced changes occurring in the methylome. In particular, genes involved in primary metabolic processes, immune functions and the p53 pathway seems to be undergoing changes in their methylation/hydroxymethylation status under the influence of simulated microgravity. An early study on lymphoblastoid cells subjected to 48 hours of simulated microgravity by Degan et al. reported decrease of cellular ATP content, suggesting a simulated microgravity induced alteration in cellular metabolism [56]. It remains to be seen how simulated microgravity induced changes over the methylation levels of p53 effector genes play in TK6 which expresses the wild-type p53 [57], a tumor suppressor functioning extensively in the DNA repair pathway. Reduction of global methylation has been proposed to be a hallmark of genomic instability [14,58] and it remains to be seen if the extensive loss-of-5mC induced by simulated microgravity reported in this study has any functional implications. Another finding in TK6 cell line, which was originally derived from a patient with T-blast crisis [32,59,60] and could potentially harbor progenitor forms of lymphocytes, pertains to changes in methylation/ hydroxymethylation patterns over genes involved in lymphocyte development and activation cultured under simulated microgravity conditions. Interestingly whole-exome sequencing has revealed similarities in the genomic content of lymphocytes and lymphoblastoid cells [61], and thus in light of our findings TK6 lymphoblastoid cells may emerge as a good model to study B and T-lymphocyte development and activation in in vitro genomic studies.
Our study also revealed that simulated microgravity could alter the expression of 370 transcripts, however only 17 of these underwent greater than 2-fold change of up/downregulation. The transcriptionally upregulated genes showed enrichment of pathways involving response to oxidative stress and negative regulation of gene expression, while the downregulated genes could be linked to pathways responsible for glucose metabolism and transcription regulation. While our study illustrated that there was no direct relationship between differentially expressed genes and changes in 5mC/5hmC over its promoters/gene bodies, we have been able to determine the methylation status of individual genes implicated in earlier studies to be affected in transcriptional or translational activities on exposure to simulated microgravity in ground based studies or in spaceflights. For instance, the voltage-dependent calcium channel L-type, alpha 1D (CACNA1D) gene transcript was observed to be differentially expressed in human T-lymphocytes subjected to microgravity conditions during spaceflight compared to ground static controls [49]. While, we observed a nominal increase at the transcript level for CACNA1D, we observed a decrease in 5mC levels over its gene body under simulated microgravity.
In another study, the Activating Transcription Factor-3 (ATF3) has been implicated to be differentially expressed upon being subjected to microgravity during spaceflight in cultured HUVEC cells [62]. Our RNA-seq data illustrated an increase of 1.3 fold in the transcript level of ATF3 and a decrease in the 5mC levels over its gene body under the influence of simulated microgravity. Interestingly, ATF3, a member of the ATF/CREB family of transcription factors, has been observed to be upregulated when cells are exposed to stress conditions [50]. On the other hand, Integrin alpha-6 (ITGA6) which is an integral cell surface protein has been observed to be down-regulated at the transcriptional scale during short-term weightlessness produced by parabolic maneuvers in human cells [51]. While RNA-seq revealed a decrease of ITGA6 transcript by more than 2-fold, we were not able to observe changes in the 5mC and 5hmC profile over its gene body or promoter, implying that possibly mechanisms other than DNA methylation might be involved in its regulation.
Some of the novel gene functions that we have linked with DNA methylation status include the F-Box Protein 17 (FBXO17), which constitutes one of the four subunits of the ubiquitinprotein-ligase complex called SKP1-cullin-F-box (SCFs) and mediates substrate specificity [63,64]. While the transcript level of FBXO17 was observed to be downregulated by 2.47 fold, the 5mC levels over the gene body of FBXO17 (chr19:39437782-39438198) decreased in TK6 cells subjected to simulated microgravity. Recently it has been demonstrated that the recruitment of F-box motif bearing homologous protein in yeast Met30 is regulated by a complex mechanism and has been implicated in stress response [65,66]. In sync with these observations, reduction of 5mC levels over gene bodies of other F-Box motif containing proteins such as FBXO31 (chr16:87421262-87421678) and FBXO42 (chr1:16674945-16675361) and promoter of FBXO5 (2024 bps upstream of TSS; chr6:153306530-153306946), was also observed though their transcripts were not differentially expressed. Interestingly, genes which function as molecular mechano-sensors like Vinculin (VCL) or mediate stress-signal transduction events like Tetraspan-5 (TSPAN5) were also seen to undergo changes in its gene methylation levels and expression. Similarly, other genes involved in the Metabolic process (GO:0008152) like Cystathionine gamma-lyase (CTH), Phospholipid scramblase-1 (PLS1) Microtubule-associated serine/threonine-protein kinase-1 (MAST1), Zinc finger protein castor homolog-1 (CASZ1), CMP-N-acetylneuraminate-beta-galactosamide-alpha-2,3-sialyltransferase-1 (ST3GAL1), AT-rich interactive domain-containing protein-5B (ARID5B), Mitogen-activated protein kinase-13 (MAP3K13) and Perilipin-2 (PLIN2) were implicated in this study contributing to mechano-stress response. Though our study does not show a global correlation between methylation status and transcriptional activity, the simulated microgravity induced changes over SPG20 (a gene implicated in endosomal trafficking and mitochondrial functions) recapitulates the conventional theory of decrease in promoter methylation corresponding to elevated gene activity. This novel finding suggests that methylation-dependent transcriptional activity is not a genome-wide phenomenon, instead it may be applicable for specific genes.
Thus, in conclusion we believe that 48 hours of treatment with simulated microgravity triggered changes in the transcriptome particularly involving biological processes such as negative regulation of transcription, response to stress and reduction in carbohydrate metabolic processes. This study revealed that simulated microgravity influenced alteration of genome-wide 5mC and 5hmC patterns, however no correlation was found between DMRs/DHMRs situated at gene bodies and promoters and their transcriptional status. While it has been long held that genes with methylated promoters are transcriptionally silent, recent studies have uncovered the association of methylated gene promoters with both transcriptionally active and inactive genes [20,21,[67][68][69][70]. On the other hand, gene body methylation has been observed to be positively correlated with gene expression in some studies [71,72] and no such correlation has been found in others [22,[73][74][75]. Recent deep-sequencing based explorations have challenged the traditional paradigm and illustrated complexities of the nature of relationship between DNA methylation and gene expression [19][20][21][22][23][24][25]. It is also conceivable that pronounced alterations in epigenetic patterns may take place in cells subjected to prolonged microgravity environments. The ground-based microgravity simulators like the one used in our study have undoubtedly enhanced our understanding of microgravity but it has to be pointed out that the principle of "simulating" microgravity involves changing the direction of Earth's gravity subjected to the samples through continuous rotation and represent "functional near weightlessness". While this is the first study to profile the simulated microgravity induced changes in 5mC/5hmC patterns and gene expression simultaneously providing a perspective of epigenetic alterations we could expect during short-term exposures, our understanding is far from complete. We believe that genes involved in altered biological processes identified in this study will be of considerable interest and provide a valuable resource for future investigations. Finally, in the interest of astronauts who are exposed to microgravity for prolonged periods of time, future studies should focus on performing time course experiments monitoring the influence of "real" and "simulated" microgravity exposure on a variety of models to determine the precise effects of microgravity on the epigenome Supporting Information of genes whose gene body was found to be associated with DMRs and DHMRs (gain-of-5mC/hmC and loss-of-5mC/hmC) (TIF) S1 Table. Sequencing summary quality statistics. (XLSX) S2 Table. List of simulated microgravity-induced DMRs undergoing hypermethylation. For every DMR identified, a description of the genomic features found in this region has been provided. The columns represent the following information: (A) Genomic coordinates of the region defined as a DMR and (B) ENCODE IDs of features (such as gene, transcript, pseudogene, non-coding RNA or other regulatory feature) present in the region. (XLSX) S3 Table. List of simulated microgravity-induced DMRs undergoing hypomethylation. For every DMR identified, a description of the genomic features found in this region has been provided. The columns represent the following information: (A) Genomic coordinates of the region defined as a DMR and (B) ENCODE IDs of features (such as gene, transcript, pseudogene, non-coding RNA or other regulatory feature) present in the region. (XLSX) S4 Table. GREAT Ontology Summary Statistics for hypermethylated DMRs. The columns represents the respective ontology name, term name / identifier, term description, binomial rank, binomial p-value (uncorrected), binomial Bonferroni corrected p-value, binomial FDR q-value and names of gene hits generated by GREAT version 3.0.0; Species assembly: hg19 and association rule: Basal+extension: 5000 bp upstream, 1000 bp downstream, 1000000 bp max extension, curated regulatory domains included. (XLSX) S5 Table. GREAT Ontology Summary Statistics for hypomethylated DMRs. The columns represents the respective ontology name, term name / identifier, term description, binomial rank, binomial p-value (uncorrected), binomial Bonferroni corrected p-value, binomial FDR q-value and names of gene hits generated by GREAT version 3.0.0; Species assembly: hg19 and association rule: Basal+extension: 5000 bp upstream, 1000 bp downstream, 1000000 bp max extension, curated regulatory domains included. (XLSX) S6 Table. List of simulated microgravity-induced DHMRs undergoing hyper-hydroxymethylation. For every DHMR identified, a description of the genomic features found in this region has been provided. The columns represent the following information: (A) Genomic coordinates of the region defined as a DHMR and (B) ENCODE IDs of features (such as gene, transcript, pseudogene, non-coding RNA or other regulatory feature) present in the region. (XLSX) S7 Table. List of simulated microgravity-induced DHMRs undergoing hypo-hydroxymethylation. The columns represent the following information for each identified DHMR: (A) Genomic coordinates of the region defined as a DHMR and (B) ENCODE IDs of features (such as gene, transcript, pseudogene, non-coding RNA or other regulatory feature) present in the region. (XLSX) S8 Table. GREAT Ontology Summary Statistics for hyperhydroxymethylated DHMRs. The columns represents the respective ontology name, term name / identifier, term description, binomial rank, binomial p-value (uncorrected), binomial Bonferroni corrected p-value, binomial FDR q-value and names of gene hits generated by GREAT version 3.0.0; Species assembly: hg19 and association rule: Basal+extension: 5000 bp upstream, 1000 bp downstream, 1000000 bp max extension, curated regulatory domains included. (XLSX) S9 Table. List of Differentially Expressed Genes induced by simulated microgravity. The columns represent geneID, name of gene from UCSC Genome Browser (duplicates exist because multiple geneID can map to same gene), chromosome location, strand: + or-, transcription start position, transcription end position, the log2-fold-change of gene expression, the average log2-counts-per-million of comparison, p-value of comparison, false discovery rate (corrected p-value) of comparison. (XLSX) S10 Table. Pathway Analysis of simulated microgravity induced differentially up and downregulated genes. (XLSX)