Evaluation of DNA Methylation Profiles of LINE-1, Alu and Ribosomal DNA Repeats in Human Cell Lines Exposed to Radiofrequency Radiation

1. Introduction

In 2011, radiofrequency electromagnetic fields (RF-EMFs), ranging from 30 kHz to 300 GHz of the electromagnetic spectrum, were classified by the International Agency for Research on Cancer (IARC) as a possibly carcinogenic to humans (risk group 2B) [1]. This type of radiation is emitted by everyday devices including mobile phones, radios, televisions, medical equipment, and by many other sources in occupational and general environmental settings. There is a growing concern about possible harmful health effects, since human exposures are ubiquitous and are rapidly increasing in the world. Genotoxic events are considered to be an initial step in carcinogenesis, so during the last decade many investigations have been carried out on the possible genotoxic effects of RF-EMFs, as summarized by several reviews [2,3,4,5]. Up to date, no consistent data have emerged about the DNA damage induction, as evaluated through analysis of various endpoints such as induction of chromosomal aberrations, micronuclei, sister chromatid exchange, mutations, single and double DNA strand breaks, and DNA degradation. Interestingly, the number of published studies reporting increased genetic damage in cells exposed to RF-EMF is inversely proportional with the number of quality control measures [5], so the inconsistency of the results could be largely due to the poor quality of many studies [6]. Recent evidence indicates that carcinogenesis can also be induced by non-genotoxic agents through epigenetic alterations [7], which can affect gene function leading to cellular neoplastic transformation without changing the DNA sequence. Therefore, to evaluate the potential carcinogenicity of RF-EMF exposure, possible epigenetic effects should be investigated.

It was reported that the earth magnetic field is involved in DNA methylation (DNAm) regulation [8,9], therefore it was reasonable to suspect that the exposure to man-made electromagnetic fields, which are greater than the natural one, might interfere with DNAm causing a dysfunction. Accordingly, previous studies reported alterations of DNAm upon exposure to RF-EMF [10,11]. Mokarram et al. showed that 900 MHz RF-EMF altered DNAm of estrogen receptor in rat colon cells, while Kumar et al. reported a global decrease of DNAm levels in rat hippocampus exposed to RF-EMF. Importantly, in the last study the effect on DNAm was more pronounced with increasing frequency (900MHz-2450MHz) and exposure time (one-month to six-month exposure group) [11].

In the present study we focused on DNAm, since it is the epigenetic mechanism that has been more extensively studied in relation to environmental exposures [12]. In particular, we evaluated DNAm of repetitive DNA elements (RE-DNA), which occupy a large part of the genome and whose deregulation has been implicated in carcinogenesis [13,14].

We focused on three types of RE-DNA: (1) long interspersed nuclear elements -1(LINE-1), (2) Alu short interspersed nuclear elements (SINEs), (3) the DNA sequences encoding for RNA ribosomal (rDNA). LINE-1s are the only active retroelements, and constitute approximately 17% of the human genome [15]. A full-length LINE-1 element is about 6 kb and is composed of a 5′ untranslated region (5′UTR) containing sense and antisense promoters, two open reading frames (ORF1 and ORF2), encoding proteins involved in the retrotransposition process, and 3′ untranslated region (3′UTR) with a polyadenylation site [16]. After transcription, the LINE-1 retroelement can be inserted into another genomic site causing genetic instability [17]. The methylation of the CpG island in LINE-1 5′UTRs plays an important role in the repression of LINE-1 transcription and retrotransposition, moreover it can affect the regulation of the expression of surrounding genes [18,19,20]. Alterations of the LINE-1 5′UTR promoter methylation have been observed in human cancer cells [21,22,23] and they are considered a promising biomarker of cancer development [14]. These alterations can be a consequence of the exposure to various environmental stressors [24,25] and it has been proposed that they should be included in health risk assessment of environmental factors [26]. Alu elements are primate specific SINE, which comprise about 11% of the genome. They are 300 base pair sequences and are subdivided in three families (J, S, Y) due to different consensus sequences and different times of appearance in the human lineage. Alu elements are nonautonomous retrotransposons, however they can copy and paste themselves through the enzymatic activity of LINE-1-ORF2 [27]. Methylation represses Alu transcription and activity [28] and aberrant Alu methylation has been found in various kinds of tumours [29]. rDNAs are the most abundantly expressed genes in the eukaryotic genome. In the human genome, rDNAs (encoding 18S, 5.8S and 28S rRNA) are tandemly arrayed in a head-to-tail fashion, in nucleolar organizer regions (NORs) on the short arms of five pairs of acrocentric autosomal chromosomes (13, 14, 15, 21 and 22) [30]. Each array consists of multiple rDNA repeat units that vary in number among individuals and chromosomes, ranging from 60 to >800, with a mean copy number of 400 in a diploid human genome. Each unit contains a single transcription unit (~13 kb), encoding a 45S rRNA precursor, which is processed to form the 18S, 5.8S and 28S rRNA molecules, and a nontranscribed spacer (~30 Kb) (IGS) [31]. IGS contains promoter, multiple repetitive sequences and regions producing non-coding RNAs involved in various cellular processes, including stress reactions. A single unit contains more than 1500 CpG sites [32] and shows high levels of methylation and particular distribution of methylation on DNA strand. The methylation of promoters abolishes rDNA transcription by inhibiting the assembly of the transcription initiation complex, indeed promoters of active rRNA genes are usually hypomethylated and are associated with acetylated histones [33]. Recent evidence indicates that the methylation of gene body maintains rDNA transcription, by preventing enrichment of repressive histone modifications [34]. It was reported that the methylation profile of DNA sequences encoding ribosomal RNAs (rRNA) plays a role in aging [32,35] and in cancer onset [36,37]. Some studies have been conducted to ascertain whether the exposure to extremely low frequency magnetic fields (ELF-EMF) could affect DNAm of RE-DNA [38,39,40,41], but so far this issue has been little explored in relation to radiofrequency exposure.

The aim of the present study is to assess this endpoint in human cells exposed to 900 MHz GSM modulated RF-EMF, which is an important frequency for mobile communication electromagnetic radiation emitted from GSM mobile phones and from their base station antennas. To this purpose, DNA methylation profiles of LINE-1 5′-UTR promoter, Alu sequences, and three regions of rDNA (rDNA promoter, 18S and 28S) were investigated in three human cells lines (Hela, Be(2)C, and SHSY5Y), exposed to RF-EMF or to sham (control).

2. Results

2.1. Validation of the assay for the measurement of RE-DNAm

To measure DNAm of RE-DNA, we used a bisulfite-targeted sequencing assay as described in Materials and Methods. We checked the assay for bisulfite-PCR-bias, a phenomenon associated to preferential amplification of bisulfite-converted DNA according to its original DNAm status. To this aim we analysed samples at known DNAm percentages (0%, 25%, 50%, 75% and 100%) and then assessed the correlation between observed and expected DNAm values for all CpGs of each RE target (Supplementary Materials). A high correlation between expected and observed DNAm values was observed for all CpGs assessed in the LINE-1 target (Pearson R > 0.94, p-value<0.01) (Figure S1). Similarly, for rDNA targets (namely RiboProm1, RiboProm2, 18S1 and 28S) we observed a highly significant correlation for the vast majority of CpGs assessed (Pearson R > 0.96, p-value<0.01) (Figure S1). The correlation between expected and observed DNAm values for CpGs belonging to the Alu target tended to be lower respect to the other targets, with some of the CpG sites not reaching statistical significance (Figure S1). CpGs having a low (Pearson R < 0.7) or not-significant correlation (p-value>0.05) were considered not reliable and were excluded from the analyses.

2.2. Effects of electromagnetic radiation on RE DNAm

To study the effects of radiofrequency exposure on epigenetic regulation of repetitive elements, HeLa, BE(2)C and SH-SY5Y cell lines were exposed to 900 MHz GSM modulated RF-EMF at SAR of 1 W/kg. DNAm of LINE-1, Alu and rDNA elements was evaluated.

The analysis of LINE-1 DNAm in HeLa cell line after radiofrequency exposure showed a significant, albeit small, hypermethylation for CpGs X70 ($\Delta DNAm:0.019;$p-value:0.033), X183 ($\Delta DNAm:0.017$; p-value:0.022), X190 ($\Delta DNAm:0.022$; p-value:0.007) and X217 ($\Delta DNAm:0.015$; p-value:0.009) compared to sham-exposed samples. In SH-SY5Y and BE2C cell lines we did not find significant changes in any of the CpGs assayed in LINE-1 target, although a trend towards hypermethylation was observed in SH-SY5Y cell line (Figure 1, Table S1).

Figure 1. DNAm profiles of LINE-1 target in cell lines exposed to RF-EMF. Lineplot showing average DNAm ± standard deviation for every CpGs assayed in LINE-1-targeting amplicon for SH-SY5Y, BE(2)C and HeLa cell lines exposed to 900MHz radiation (red) or to sham (blue). CpGs with p-value<0.01 are highlighted with dark grey boxes; CpGs with p-value between 0.01 and 0.05 are highlighted with light grey boxes.

Figure 1. DNAm profiles of LINE-1 target in cell lines exposed to RF-EMF. Lineplot showing average DNAm ± standard deviation for every CpGs assayed in LINE-1-targeting amplicon for SH-SY5Y, BE(2)C and HeLa cell lines exposed to 900MHz radiation (red) or to sham (blue). CpGs with p-value<0.01 are highlighted with dark grey boxes; CpGs with p-value between 0.01 and 0.05 are highlighted with light grey boxes.

The DNAm profile of the promoter (RiboProm1 and RiboProm2) and of the body (18S1 and 28S) regions of rDNA was affected by radiofrequency exposure in different ways in the three cell lines examined (Figure 2 and Figure 4; Table S1).

In SH-SY5Y cell line, a large group of CpGs, mostly on the 5’ of RiboProm2 target, was found significantly hypermethylated following RF-EMF exposure ($\Delta DNAm$ ranging from 0.014 to 0.030 p-values ranging from 0.012 to 0.044) (Figure 2B, Table S1). A larger hypermethylation trend was observed for the amplicon targeting 18S, with CpGs X135 ($\Delta DNAm:0.079;$p-value:0.041), X179 ($\Delta DNAm:0.087;$p-value:0.042) and X190 ($\Delta DNAm:0.070;$p-value:0.042) reaching statistical significance (Figure 3A, Table S1). On the other hand, no significant DNAm changes were found for both RiboProm1 (Figure 2A, Table S1) and 28S (Figure 3B, Table S1).

In BE(2)C neuroblastoma cell line, we observed hypomethylation on both amplicons targeting the rDNA promoter following RF-EMF exposure. DNAm changes in RiboProm1 amplicon were small but significant for the majority of the target’s CpGs ($\Delta DNAmrangingfrom$ -0.026 to -0.014; p-values ranging from 0.004 to 0.045), whereas for RiboProm2, only CpG X159 ($\Delta DNAm:$-0.033; p-value:0.048) reached statistical significance (Figure 2, Table S1). The exposure of BE(2)C neuroblastoma cells to RF-EMF did not induce any significant DNAm changes for the CpGs belonging to 18S; while it promoted hypomethylation of CpGs in the 3’ end of 28S amplicon with CpGs X197 ($\Delta DNAm$:-0.018; p-value:0.009) and X205 ($\Delta DNAm$:-0.014; p-value:0.019) reaching statistical significance (Figure 3, Table S1).

In HeLa cell line, DNAm profiles did not show any significant changes for any rDNA-targeting amplicons with the exception of 28S. Indeed, radiofrequency exposure stimulated hypomethylation of CpGs on the 5’end of 28S amplicon, with CpGs x29 ($\Delta DNAm$:-0.021 p-value:0.010), x72 ($\Delta DNAm$:-0.025; p-value:0.002), x83 ($\Delta DNAm$:-0.017; p-value:0.049) and x99 ($\Delta DNAm$:-0.015; p-value:0.031) reaching statistical significance (Figure 3B, Table S1).

Figure 2. DNAm profiles of rDNA promoter targets in cell lines exposed to RF-EMF. Lineplots showing average DNAm ± standard deviation for every CpGs assayed in A) RiboProm1 and B) RiboProm2 amplicons for SH-SY5Y, BE(2)C and HeLa cell lines exposed to 900MHz radiation (red) or to sham (blue). CpGs with p-value<0.01 are highlighted with dark grey boxes; CpGs with p-value between 0.01 and 0.05 are highlighted with light grey boxes. CpGs which were discarded due to poor calibration are highlighted with red boxes.

Figure 2. DNAm profiles of rDNA promoter targets in cell lines exposed to RF-EMF. Lineplots showing average DNAm ± standard deviation for every CpGs assayed in A) RiboProm1 and B) RiboProm2 amplicons for SH-SY5Y, BE(2)C and HeLa cell lines exposed to 900MHz radiation (red) or to sham (blue). CpGs with p-value<0.01 are highlighted with dark grey boxes; CpGs with p-value between 0.01 and 0.05 are highlighted with light grey boxes. CpGs which were discarded due to poor calibration are highlighted with red boxes.

Figure 3. DNAm profiles of rDNA body in cell lines exposed to RF-EMF. Lineplots showing average DNAm ± standard deviation for every CpGs assayed in A) 18S and B) 28S amplicons for SH-SY5Y, BE(2)C and HeLa cell lines exposed to 900MHz radiation (red) or to sham (blue). CpGs with p-value<0.01 are highlighted with dark grey boxes; CpGs with p-value between 0.01 and 0.05 are highlighted with light grey boxes. CpGs which were discarded due to poor calibration are highlighted with red boxes.

Figure 3. DNAm profiles of rDNA body in cell lines exposed to RF-EMF. Lineplots showing average DNAm ± standard deviation for every CpGs assayed in A) 18S and B) 28S amplicons for SH-SY5Y, BE(2)C and HeLa cell lines exposed to 900MHz radiation (red) or to sham (blue). CpGs with p-value<0.01 are highlighted with dark grey boxes; CpGs with p-value between 0.01 and 0.05 are highlighted with light grey boxes. CpGs which were discarded due to poor calibration are highlighted with red boxes.

Finally, no DNAm changes were generally found for the CpGs of the amplicons targeting the body of Alu repeat unit in any of the cell lines assayed with the exception of a slight, but significant hypermethylation of CpG X118 ($\Delta DNAm$:-0.005; p-value:0.038) in HeLa cell lines (Figure 4, Table S1).

Figure 4. DNAm profiles of Alu repeat unit in cell lines exposed to radiofrequencies. Lineplots showing average DNAm ± standard deviation for every CpGs assayed in A) 18S and B) 28S amplicons for SH-SY5Y, BE(2)C and HeLa cell lines exposed to 900MHz radiation (red) or to sham (blue). CpGs with p-value<0.01 are highlighted with dark grey boxes; CpGs with p-value between 0.01 and 0.05 are highlighted with light grey boxes. CpGs which were discarded due to poor calibration are highlighted with red boxes.

Figure 4. DNAm profiles of Alu repeat unit in cell lines exposed to radiofrequencies. Lineplots showing average DNAm ± standard deviation for every CpGs assayed in A) 18S and B) 28S amplicons for SH-SY5Y, BE(2)C and HeLa cell lines exposed to 900MHz radiation (red) or to sham (blue). CpGs with p-value<0.01 are highlighted with dark grey boxes; CpGs with p-value between 0.01 and 0.05 are highlighted with light grey boxes. CpGs which were discarded due to poor calibration are highlighted with red boxes.

2.3. Effects of electromagnetic radiation on RE DNAm epihaplotype diversity

To further explore the effects elicited by radiofrequency radiation on the DNAm profiles of repetitive elements in the various cell lines assayed, we evaluated DNAm epihaplotype diversity.

AmplimethProfiler tool has a built-in, qiime-based analytical pipeline that allows to calculate Shannon diversity indexes via rarefaction method for each target (Materials and Methods). Filtering of samples with low coverage (Materials and Methods) resulted in the removal of one HeLa cell replicate for all the analysed regions.

HeLa cells showed a small but significant (p-value:0.01) increase in epihaplotype diversity for LINE-1-targeting amplicon upon exposure to radiofrequency radiation (Figure 5, Table S2). Similarly, a small and nearly significant (p-value:0.05) increase in epihaplotype diversity for LINE-1 was observed in radiofrequency-exposed SH-5YSY neuroblastoma cells (Figure 5, Table S2).

Figure 5. Epihaplotype diversity of DNAm profiles of LINE-1 amplicon in cell lines exposed to radiofrequencies. Boxplots reporting Shannon diversity indexes for DNAm epihaplotypes for LINE-1 amplicon for SH-SY5Y, BE(2)C and HeLa cell lines exposed to 900MHz radiation (orange) or to sham (blue). (* p-value<0.05).

Figure 5. Epihaplotype diversity of DNAm profiles of LINE-1 amplicon in cell lines exposed to radiofrequencies. Boxplots reporting Shannon diversity indexes for DNAm epihaplotypes for LINE-1 amplicon for SH-SY5Y, BE(2)C and HeLa cell lines exposed to 900MHz radiation (orange) or to sham (blue). (* p-value<0.05).

For rDNA-targeting amplicons, we found a significant decrease in epihaplotype diversity index for amplicons targeting the proximal rDNA promoter RiboProm2 (p-value:0.023) (Figure 6, Table S2) and 28S (p-value:0.025) (Figure 7, Table S2) upon radiofrequency exposure of BE(2) neuroblastoma cell line. We did not find any significant changes in epihaplotype diversity in SH-SY5Y and HeLa cell lines.

Figure 6. Epihaplotype diversity of DNAm profiles of rDNA promoter in cell lines exposed to RF-EMF. Boxplots reporting Shannon diversity indexes for DNAm epihaplotypes of RiboProm1 and RiboProm2 amplicons for SH-SY5Y, BE(2)C and HeLa cell lines exposed to 900MHz radiation (orange) or to sham (blue). (* p-value<0.05).

Figure 6. Epihaplotype diversity of DNAm profiles of rDNA promoter in cell lines exposed to RF-EMF. Boxplots reporting Shannon diversity indexes for DNAm epihaplotypes of RiboProm1 and RiboProm2 amplicons for SH-SY5Y, BE(2)C and HeLa cell lines exposed to 900MHz radiation (orange) or to sham (blue). (* p-value<0.05).

Figure 7. Epihaplotype diversity of DNAm profiles of rDNA body in cell lines exposed to RF-EMF. Boxplots reporting Shannon diversity indexes for DNAm epihaplotypes of 18S and 28S amplicons for SH-SY5Y, BE(2)C and HeLa cell lines exposed to 900MHz radiation (orange) or to sham (blue). (* p-value<0.05).

Figure 7. Epihaplotype diversity of DNAm profiles of rDNA body in cell lines exposed to RF-EMF. Boxplots reporting Shannon diversity indexes for DNAm epihaplotypes of 18S and 28S amplicons for SH-SY5Y, BE(2)C and HeLa cell lines exposed to 900MHz radiation (orange) or to sham (blue). (* p-value<0.05).

Finally, for Alu-targeting amplicons, no significant differences in epihaplotype diversity indexes were found in any of the cell lines assayed (Figure 8).

Figure 8. Epihaplotype diversity of DNAm profiles of Alu repeat unit in cell lines exposed to RF-EMF. Boxplots reporting Shannon diversity indexes for DNAm epihaplotypes of Alu amplicon for SH-SY5Y, BE(2)C and HeLa cell lines exposed to 900MHz radiation (orange) or to sham (blue).

Figure 8. Epihaplotype diversity of DNAm profiles of Alu repeat unit in cell lines exposed to RF-EMF. Boxplots reporting Shannon diversity indexes for DNAm epihaplotypes of Alu amplicon for SH-SY5Y, BE(2)C and HeLa cell lines exposed to 900MHz radiation (orange) or to sham (blue).

3. Discussion

In this work, we evaluated the impact of RF-EMF exposure on DNAm of LINE-1, Alu and rDNA repetitive elements in human neuroblastoma and cervical cancer cell lines. The effects on DNAm were evaluated in terms of both alterations in the average DNAm level and alterations of DNAm profiles (epihaplotype) richness and distribution.

Our data indicate that RF-EMF exposure can induce alterations in the average DNAm level of LINE-1 CpGs (limited to Hela cell line), and of rDNA (for all three cell lines tested). No effect was instead found for the methylation of CpGs of Alu repeat in any of the cell lines. Previously, Alu methylation was evaluated by Benassi et al. [40] in SH5H5Y cells exposed to ELF-EMF and no effect was found, similarly to our observation.

With regard to LINE-1, we observed an increase in DNAm of some CpG in Hela cells; this finding could be consistent with the observation, reported in a previous paper [42], of a significant reduction in the mRNA levels of LINE-1 in Hela cells exposed to 900 MHz. Indeed, an increased DNAm is usually associated with decreased RNA transcription.

Differently, no effect was found on LINE-1 methylation in neuroblastoma cells (SHSY5Y and BE(2)C cells). Similar investigations were previously carried out in SHSY5Y and BE(2)C cells exposed to ELF-EMF (50Hz, 1 mT) in two different studies. In agreement with our results, even though concerning a different type of electromagnetic field, no significant DNAm changes were found neither in exposed BE(2)C [39] nor in exposed SH-5YSY cells [40].

With regard to rDNA methylation, alterations were observed in all three cell lines: exposed Hela cells showed hypomethylation of some CpG in the 3’ end of 28S amplicon target, exposed SH-SY5Y cells showed hypermethylation of the 5’ of RiboProm2 and 18S targets, while exposed BE(2)C neuroblastoma cells showed hypomethylation of the 5’-end of RiboProm1, the 5’-end of RiboProm2 and the 3’-end of 28S amplicon targets. It is well known that rDNA methylation patterns contribute to the regulation of the expression of ribosomal RNA (rRNA) related to the translational rate, which depends on cell-specific homeostasis requirement [43]. Therefore, it is not surprising that the three cell lines analysed showed different changes in rDNA methylation in response to RF-EMF exposure.

To date, DNAm of rDNA has never been analysed in exposed cells to ELF-EMF or to RF-EMF but it has been investigated in both physiological and pathological conditions including studies focused on cancer [36,37], neurodegenerative diseases [44], ageing and ageing-related conditions [32,35]. For instance, changes in rDNA methylation have been described in human blood cells and hepatocytes with aging [45], in blood samples from persons with Down Syndrome [46] and in brain samples from patients with Alzheimer’s disease [44,47].

In addition to assessing the impact of radiofrequncies on the average DNAm level, we evaluated their possible effects on the organisation of methylated and unmethylted CpG sites through the analysis of epihaplotype richness and distribution.

As far as we know, the analysis of DNAm epihaplotypes, by calculating the Shannon diversity index, has never been performed in cells EMF-exposed. Our data indicate that RF-EMF exposure can induce changes in the DNAm epihaplotype diversity. In some cases, such as for LINE-1 DNAm in Hela and SH5H5Y cells, epihaplotype diversity increased in exposed samples, whereas in other cases, such as rDNA DNAm in BE(2) cells, it decreased. Changes in distribution and richness of epihaplotypes may indicate a dysfunction in the maintenance of epigenetic patterns which can contribute to epigenetic instability. Changes of Shannon diversity index have been found in aging [48,49] and in progeroid syndromes such as Down Syndrome [46].

As a whole, our data indicate that RF-EMF exposure affects DNAm of tested RE and might induce epigenetic instability. These data add further evidence to our knowledge on the effects of RF-EMF on DNAm [10,11]. A weakness of the present study is the limited number of analysed samples and that the use of cell lines may not exactly reflect what happens in vivo. Analysis on primary tissues and stem cells would provide further knowledge.

The issue of what molecular mechanisms underlie the biological effects of EMFs is still an open question [50,51,52]. One hypothesis is that EMFs induce changes in the energy levels of certain molecules affecting the concentration of free radicals, such as reactive oxygen species (ROS) [53,54]. Several studies documented the ability of RF-EMF to induce an increased production of ROS [55,56,57,58,59,60,61] and similar results have been reported for ELF-MF [62,63,64,65]. It is emerging that ROS can affect both genome-wide and site-specific DNAm patterns, causing hypomethylation or hypermethylation [66,67,68]. The action of ROS can take place through multiple mechanisms: ROS can cause the oxidation of DNA bases and enzymes, alter the regulation of DNA methyltransferase (DNMTs) and demethylase enzymes, alter the availability of metabolic cofactors [69,70]. In our research, we did not analyse ROS levels in exposed cells, however some in vitro studies showed increased ROS levels in SH5H5Y cells after RF-EMF exposure [71,72]. Therefore, the effects on DNAm reported in this study might be linked to the induction of ROS production by RF exposure. Further studies should be carried out to verify this hypothesis.

Moreover, our results indicate that the effect on DNAm of LINE-1 and rDNA strongly depends on the cellular system considered and on the genetic background; in agreement with the fact that epigenetic modifications are tissue-specific and show inter-individual variability.

Previous studies showed that the biological model is relevant in EMF exposure studies. For example, Martin et al. [73] studied the expression of three genes in three primary cultures of human keratinocytes subjected to identical exposure conditions (60.4 GHz), and found 3 different expression profiles. They speculated that this fact could be due to epigenetic modifications specific to each cellular model used. Similar results were found by Del Re et al [42] studying the expression of some repetitive elements in three different cell lines.

Collectively, these data indicate that many different models should be studied to understand the epigenetic effects of EMF exposure on the biological samples.

In conclusion, our results suggest that RF-EMF exposure can affect DNAm of LINE-1 and rDNA sequences, which are highly represented in the human genome, and that the effects depend on cell context (cell genome, cell type). Further investigations on epigenetic effects of RF-EMF would be useful both for a better understanding of RF-EMF risks and safety assessment.

4. Materials and Methods

4.1. Exposure system

The exposure system was kindly loaned us by Prof. Laura Calzà (University of Bologna). The exposure system and the dosimetry were previously described in great detail [74,75]. Briefly, the system is composed of two twin “Transverse Electro Magnetic (TEM) cells” with two independent power supply chains (two GSM signal generators and two EMF amplifiers) to obtain sham and exposed samples at the same time and it was positioned inside a Heraeus incubator (B-5060, Heraeus, Hanau, Germany) (37◦ C, 5% CO2). Each TEM cell has a squared cross-section (14 cm wide), and can contain 6 four-multiwell plates (66 mm x 66 mm external size). The RF-EMF (900 MHz; GSM basic with pulsed modulation at 217 Hz) is uniformly absorbed in the two wells of the two symmetric plates closer to the inner copper septum of the TEM cell resulting in a SAR of 1 W/Kg. Therefore, only the septum-adjacent wells were used in the experiments according to Del Vecchio et al. (2009a, 2009b). The system was controlled by LabVIEW 6.1 software (National Instruments, Austin, TX). Experiments were conducted in blind modality. In each assay, cells were plated and allowed to adhere overnight, then were exposed to 900 MHz GSM modulated RF-EMF at SAR of 1 W/kg or to sham. After 48 h exposure, DNA was extracted from each sample.

4.2. Cell culture

HeLa cells, neuroblastoma BE(2)C cells, and neuroblastoma SH-SY5Y cells were kindly provided by Prof. Perini (University of Bologna, Italy). All cell lines were grown in DMEM medium (EuroClone, Milano, Italy), containing 10% heat-inactivated fetal bovine serum (EuroClone, Milano, Italy), 100 UI/ml penicillin (Sigma, Ronkonkoma, NY, USA) and 100 μg/ml streptomycin (Sigma), under a 5% CO2 humified atmosphere at 37°C. In each assay, cells were plated in four-well plates (15 mm diameter; 1 ml/well; culture area 1.9 cm²) (Thermo Scientific Nunc, A/S Roskilde, Denmark) at a density of 80000 cells/well and allowed to adhere overnight, then the plates were randomly assigned to each Transverse Electro Magnetic (TEM) chamber and exposed to 900 MHz GSM-modulated RF-EMF at specific absorption rate (SAR) of 1 W/kg or to sham.

4.3. DNA extraction and bisulfite conversion

Genomic DNA was extracted by QIAmp DNA Mini Kit (QIAgen, Hilden, Germany) according to manufacturer’s instructions and quantified using Qubit dsDNA BR Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA). Then, 500 ng were treated with sodium bisulfite using EZ DNA Methylation Direct Kit (Zymo Research, Irvine, CA, USA) as indicated by the manufacturer. Each sample was sequenced in triplicate. DNA methylation calibration curves were prepared combining Universal Methylated and Universal Unmethylated DNA samples (Millipore, Burlington, MA, USA) to generate standards at 0, 25, 50, 75, 100% DNAm levels. Each point of the curve was sequenced in triplicate

4.4. Library Preparation and Bisulfite-targeted Sequencing

DNAm of RE-DNA was analysed using a targeted-bisulfite sequencing approach. Primers for rDNA distal (RiboProm1) and proximal (RiboProm2) promoters as well as for 5’UTR of LINE-1 were previously published [76] whereas primers targeting 18S, 28S and Alu regions were designed using MethPrimer 2.0 (Table S3). Briefly, RiboProm1 and RiboProm2 are located in the promoter region of the rDNA unit. The former is upstream of the rDNA promoter, whereas the latter encompasses the core promoter element (CP) as well as the upstream control element (UCE) of the rDNA unit. 18S and 28S targets are located at the 5’ end of their respective rDNA sequences (Figure 9A). RiboProm1, RiboProm2, 18S and 28S contain 37, 26, 27 and 30 CG sites, respectively (Figure S2). LINE-1 target is located in the 5’UTR of the repeat unit, close to the transcription start site (Figure 9B). Alu target is located in the body of the repeat unit (Figure 9C). LINE-1 and Alu repeat targets respectively contain 18 and 12 CG sites (Figure S2). Forward and reverse primers were added at each 5’ end with Nextera^TM adapter sequences, respectively TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG and GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG. In addition, a random nucleotide spacer (N) was included between Illumina adapters and primers to increase sequencing variability.

Figure 9. Structure of the rDNA, LINE-1 and Alu units and localization of target regions assessed. A) Structure and organization of the rDNA repeat unit. The position of each amplicon (colored arrows) assayed through target-bisufite sequencing is defined by its distance (in base pairs) from the transcription start site (TSS) of to the 45S pre-ribosomal RNA transcript (chr21:8202347-8222335, GRCh38/hg38). UCE = upstream control element; CP = core promoter; TSS = transcription start site; ETS = external transcribed spacer; ITS = internal transcribed spacer; 18S = 18S rRNA coding region; 5.8S = 5.8S rRNA coding region; 28S = 28S rRNA coding region. B) Structure of LINE-1 repeat unit. The amplicon (colored arrow) included in the targeted-bisulfite sequencing assay is located in the 5’UTR of the Human LINE-1 (L1.3) repetitive element (GenBank: L19088.1). UTR = untranslated region; ORF = open reading frame C) Structure structure of the Alu repeat unit. The amplicon (colored arrow) encompasses the A-rich region and covers the right arm of the repeat unit. The orientation of the arrows representing the amplicons indicates the 5’-3’ direction of the sequenced DNA molecule. bp = base pairs.

Figure 9. Structure of the rDNA, LINE-1 and Alu units and localization of target regions assessed. A) Structure and organization of the rDNA repeat unit. The position of each amplicon (colored arrows) assayed through target-bisufite sequencing is defined by its distance (in base pairs) from the transcription start site (TSS) of to the 45S pre-ribosomal RNA transcript (chr21:8202347-8222335, GRCh38/hg38). UCE = upstream control element; CP = core promoter; TSS = transcription start site; ETS = external transcribed spacer; ITS = internal transcribed spacer; 18S = 18S rRNA coding region; 5.8S = 5.8S rRNA coding region; 28S = 28S rRNA coding region. B) Structure of LINE-1 repeat unit. The amplicon (colored arrow) included in the targeted-bisulfite sequencing assay is located in the 5’UTR of the Human LINE-1 (L1.3) repetitive element (GenBank: L19088.1). UTR = untranslated region; ORF = open reading frame C) Structure structure of the Alu repeat unit. The amplicon (colored arrow) encompasses the A-rich region and covers the right arm of the repeat unit. The orientation of the arrows representing the amplicons indicates the 5’-3’ direction of the sequenced DNA molecule. bp = base pairs.

Sequencing libraries were generated using a two-step PCR approach. Briefly, 5 ng of bisulfite-converted DNA were amplified using Phusion U (ThermoFisher, Waltham, MA, USA) added with 1M Betaine (Merk, Darmstadt, Germany), 150nM forward and reverse primers, 1.75mM MgCl2 (Agena Bioscience, San Diego, CA, USA) and 200μM dNTP (ThermoFisher, Waltham, MA, USA). Thermal cycler conditions were set as follows: 1x cycle at 95°C for 1’ 40’’; 1x cycle at 98°C for 1’; 1x cycle at 58°C for 2’;1x cycle at 72°C for 1’; 36 cycles at 98°C for 10’’, 58°C for 40’’, 72°C for 20’’; 1x cycle at 72°C for 5’; hold at 4°C. Amplification products were pooled sample-wise and purified using MagSi-NGS plus beads (MagTivio BV, Nuth, The Netherlands). Afterwards, 10ul of pooled and purified samples were indexed using Illumina Nextera XT Index Set A as indicated in Nextera Library Prep Guide. Finally, indexed libraries were purified yet again and normalized to 4nM. Sequencing was performed on an Illumina Miseq Platform with Micro v2 300PE chemistry (Illumina, San Diego, CA, USA).

4.6. Sequencing Data Handling

Pair-end reads were quality-checked using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Adapter sequences were trimmed off using cutadapt [77] and finally, reads were merged together using PEAR [78] with a minimum of 20 overlapping residues and a maximum read length of 450. Assembled reads were then converted from FASTQ to FASTA using seqtk. All processing tools were compiled in anaconda2 environments.

4.7. DNA methylation data extraction

To extract CpG DNAm profiles for each target sequence, AmpliMethProfiler pipeline was applied [79]. This pipeline was developed to analyse deep bisulfite sequencing data. Briefly, for each sample, AmpliMethProfiler tool performs a quality check on sequenced reads by removing those with mean Phred < 33 (--perfQual). Subsequently, it recovers all reads with at least 80% alignment identity with bisulfite-specific primers provided (--primTresh 0.8) and retains all reads whose length differs less than 20% from that of a reference sequence provided (--threshLen 0.2). Afterwards, each read is aligned to a reference sequence provided using blastn. Alignment parameters were modified from default in order to switch off automatic low-complexity, repetitive sequences masking (--dust no). This allows to limit read loss as bisulfite-conversion and subsequent amplification reduce sequence complexity by converting unmethylated Cytosines (C) into Thymines (T). For each aligned read, AmpliMethProfiler generates a DNAm profile by assigning binary values to CpG positions that were found methylated (1) or unmethylated (0). If a CpG position contains a nucleotide other than T or C, AmpliMethProfiler assigns an unknown (2) value. Samples with a coverage lower than 400 reads were excluded from the analysis. After filtering, sequencing coverage varied between 1252-9432 (3858±2140) for 18S, 1316-11481 (5147±2462) for 28S, 809-7512 (3027±1240) for Alu, 822-7838 (3454±1477) for LINE-1, 973-7580 (3995±1539) for RiboProm1 and 864-8366 (3918±1736) for RiboProm2. Finally, average DNAm levels of each CpG units are calculated as percentage of methylated Cs on the total number of methylated and unmethylated Cs. For each cell line, average DNAm was compared between the treated and untreated group by Student t-test test using ANOVA function in R (v3.6).

4.8. Epihaplotype Analysis

AmpliMethProfiler provides information on the DNAm status of each CpG site at single DNA molecule level. The specific combination of methylated and unmethylated CpGs on a single strand of DNA defines an epihaplotype. To study epihaplotype organisation, a pipeline of analysis was developed using AmpliMethProfiler built-in functions. As the number of possible methylated and unmethylated CpG combinations largely exceeds sequencing depth, singletons, namely epihaplotypes occurring only once in all sample, were removed from epihaplotype count tables (filter_otus_from_otu_table.py –n 2). Epihaplotype composition and distribution within each sample was then evaluated using a qiime-based pipeline implemented within AmpliMethProfiler environment (qiime_analysis.py). Epihaplotype diversity was calculated through the Shannon alpha-diversity index. For each cell line, samples were grouped as either treated or control and then observations were rarefied by randomly resampling the pool of epihaplotypes at increasing depth (from 10 to median overall coverage). Shannon Indexes were calculated at a depth equal to the number of reads of the sample with the lowest coverage. However, it emerged that none of the Shannon Index rarefaction curves generated was able to reach plateau at the selected depth. Therefore, in order to increase the number of reads on which observations were rarefied, we included an additional filter and removed all samples with a maximum coverage lower than 60% of median coverage for each amplicon. Once the dataset was filtered, for every target amplicon, Shannon Indexes were compared between samples exposed to radiofrequency radiation and to sham using Student t-test.