Increased DNA methylation variability in type 1 diabetes across three immune effector cell types

Abstract

The incidence of type 1 diabetes (T1D) has substantially increased over the past decade, suggesting a role for non-genetic factors such as epigenetic mechanisms in disease development. Here we present an epigenome-wide association study across 406,365 CpGs in 52 monozygotic twin pairs discordant for T1D in three immune effector cell types. We observe a substantial enrichment of differentially variable CpG positions (DVPs) in T1D twins when compared with their healthy co-twins and when compared with healthy, unrelated individuals. These T1D-associated DVPs are found to be temporally stable and enriched at gene regulatory elements. Integration with cell type-specific gene regulatory circuits highlight pathways involved in immune cell metabolism and the cell cycle, including mTOR signalling. Evidence from cord blood of newborns who progress to overt T1D suggests that the DVPs likely emerge after birth. Our findings, based on 772 methylomes, implicate epigenetic changes that could contribute to disease pathogenesis in T1D.

ARTICLE

Received 17 Apr 2016 | Accepted 4 Oct 2016 | Published 29 Nov 2016

Dirk S. Paul1,2,*, Andrew E. Teschendorff3,4,*, Mary A.N. Dang5,*, Robert Lowe5,*, Mohammed I. Hawa5, Simone Ecker1, Huriya Beyan5, Stephanie Cunningham5, Alexandra R. Fouts6, Anita Ramelius7, Frances Burden8,9, Samantha Farrow8,9, Sophia Rowlston8,9, Karola Rehnstrom8,9, Mattia Frontini8,9,10, Kate Downes8,9,Stephan Busche11,12, Warren A. Cheung11,12, Bing Ge11,12, Marie-Michelle Simon11,12, David Bujold11,12,Tony Kwan11,12, Guillaume Bourque11,12, Avik Datta13, Ernesto Lowy13, Laura Clarke13, Paul Flicek13,Emanuele Libertini1, Simon Heath14,15, Marta Gut14,15, Ivo G. Gut14,15, Willem H. Ouwehand8,9,10,16,Tomi Pastinen11,12, Nicole Soranzo8,16, Sabine E. Hofer17, Beate Karges18,19, Thomas Meissner19,20,Bernhard O. Boehm21,22,23, Corrado Cilio7, Helena Elding Larsson7, ke Lernmark7, Andrea K. Steck6,

Vardhman K. Rakyan5,*, Stephan Beck1,* & R. David Leslie5,*

The incidence of type 1 diabetes (T1D) has substantially increased over the past decade, suggesting a role for non-genetic factors such as epigenetic mechanisms in disease development. Here we present an epigenome-wide association study across 406,365 CpGs in 52 monozygotic twin pairs discordant for T1D in three immune effector cell types. We observe a substantial enrichment of differentially variable CpG positions (DVPs) in T1D twins when compared with their healthy co-twins and when compared with healthy, unrelated individuals. These T1D-associated DVPs are found to be temporally stable and enriched at gene regulatory elements. Integration with cell type-specic gene regulatory circuits highlight pathways involved in immune cell metabolism and the cell cycle, including mTOR signalling. Evidence from cord blood of newborns who progress to overt T1D suggests that the DVPs likely emerge after birth. Our ndings, based on 772 methylomes, implicate epigenetic changes that could contribute to disease pathogenesis in T1D.

1 Medical Genomics, UCL Cancer Institute, University College London, London WC1E 6BT, UK. 2 Cardiovascular Epidemiology Unit, Department of Public Health and Primary Care, University of Cambridge, Strangeways Research Laboratory, Cambridge CB1 8RN, UK. 3 CAS Key Lab of Computational Biology, CAS-MPG Partner Institute for Computational Biology, Shanghai Institute for Biological Sciences, Chinese Academy of Sciences, Shanghai 200031, China. 4 Statistical Cancer Genomics, UCL Cancer Institute, University College London, London WC1E 6BT, UK. 5 The Blizard Institute, Barts and The London School of Medicine and Dentistry, Queen Mary University of London, London E1 2AT, UK. 6 Barbara Davis Center for Childhood Diabetes, University of Colorado School of Medicine, Aurora, Colorado 80045, USA.

7 Department of Clinical Sciences, Lund University, Skne University Hospital, SE-20502 Malm, Sweden. 8 Department of Haematology, University of Cambridge, Cambridge Biomedical Campus, Cambridge CB2 0PT, UK. 9 National Health Service Blood and Transplant, Cambridge Biomedical Campus, Cambridge CB2 0PT, UK.

10 British Heart Foundation Centre of Excellence, Cambridge Biomedical Campus, Cambridge CB2 0QQ, UK. 11 Department of Human Genetics, McGill University, Montreal, Qubec, Canada H3A 0G1. 12 McGill University and Genome Quebec Innovation Centre, Montreal, Qubec, Canada H3A 0G1. 13 European Molecular Biology Laboratory, European Bioinformatics Institute, Wellcome Genome Campus, Hinxton, Cambridge CB10 1SD, UK. 14 CNAG-CRG, Centre for Genomic Regulation, Barcelona Institute of Science and Technology (BIST), Baldiri Reixac 4, 08028 Barcelona, Spain. 15 Universitat Pompeu Fabra, Plaa de la Merc 10, 08002 Barcelona, Spain. 16 Human Genetics, Wellcome Trust Sanger Institute, Wellcome Genome Campus, Hinxton, Cambridge CB10 1SA, UK. 17 Department of Pediatrics, Medical University of Innsbruck, 6020 Innsbruck, Austria. 18 Division of Endocrinology and Diabetes, RWTH Aachen University, 52074 Aachen, Germany. 19 German Center for Diabetes Research (DZD), 85764 Neuherberg, Germany. 20 Department of General Pediatrics, Neonatology and Pediatric Cardiology, University Childrens Hospital, Heinrich Heine University of Dsseldorf, 40225 Dsseldorf, Germany. 21 Division of Endocrinology, Department of Internal Medicine I, Ulm University Medical Centre, 89081 Ulm, Germany. 22 Lee Kong Chian School of Medicine, Nanyang Technological University, Singapore 636921, Singapore. 23 Imperial College London,London SW7 2AZ, UK. * These authors contributed equally to this work. Correspondence and requests for materials should be addressed to D.S.P.(email: mailto:[email protected]

Web End [email protected] ) or to R.D.L. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 7:13555 | DOI: 10.1038/ncomms13555 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

DOI: 10.1038/ncomms13555 OPEN

Increased DNA methylation variability in type 1 diabetes across three immune effector cell types

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13555

Type 1 diabetes (T1D) is a common, organ-specic autoimmune disease that results from the progressive loss of insulin-producing b-cells in the pancreas. Genetic

predisposition and environmental factors contribute to the disease onset1. The incidence of T1D has dramatically increased in recent years (34% per annum), with the most rapid upsurge seen in children younger than ve years of age2. The increasing rate of T1D, along with disease discordance in monozygotic (MZ) twins, suggest that non-genetic factors play a major role3,4. Such factors, including viral and bacterial infections, diet, and potentially epigenetic and stochastic events, may affect disease predisposition either in utero or in early childhood when predictive autoantibodies emerge3. However, conclusive evidence about causal environmental factors in T1D pathogenesis has not been obtained to date.

Epigenetic modications, including DNA methylation, are cell type-specic and induce stable changes in gene expression that are heritable during cell division. DNA methylation occurs at cytosine residues mainly in the context of CpG dinucleotides, and is generally associated with transcriptional silencing5. It can contribute to disease development and progression through its inuence on gene expression, and function as mediator in response to environmental stimuli6. In systematic epigenome-wide association studies (EWASs), DNA methylation levels are typically measured at hundreds of thousands of CpG sites across individuals in a case-control cross-sectional cohort. CpG sites are then associated with disease status, and differences in DNA methylation levels between cases and controls are recorded7,8. However, it has to be noted that the meaningful interpretation of EWAS ndings is impeded by several confounding factors, in particular cellular heterogeneity in accessible sample material (for example, peripheral blood) and genetic heterogeneity between individuals7,9.

Disease-associated CpG sites can be identied using different analytical approaches (Fig. 1). Most EWASs have sought to quantify differences in mean DNA methylation at CpG sites

between cases and controls, that is, differentially methylated CpG positions (DMPs). In recent years, DMPs have been shown to associate with a multitude of complex traits and diseases, including blood pressure10, triglyceride levels11, pain sensitivity12, schizophrenia13, rheumatoid arthritis14 and T1D (refs 15,16). However, the difference in mean DNA methylation at these CpGs is often small (o5%), raising challenges to their biological interpretation.

In parallel, the potential importance of increased DNA methylation variability has been noted in cancerous tissue1722. Differentially variable CpG positions (DVPs) are heterogeneous outlier events that occur mainly, if not exclusively, in disease cases (Fig. 1). DVPs usually involve larger shifts in DNA methylation (410%), albeit in a smaller number of cases. For example, DVPs have recently been identied in precursor cervical cancer lesions that are predictive of progression to neoplasia when compared with matched control tissue21. The contribution of such DNA methylation outliers in non-cancerous tissues has not yet been evaluated. Further, the distinct functional characteristics of DVPs compared with those of DMPs have not been fully appreciated.

In this study, we determine differential DNA methylation in 52 MZ twin pairs discordant for T1D. In these twin pairs, we perform an EWAS in immune cells known to act as key drivers in the disease process, namely CD4 T cells, CD19 B cells and

CD14 CD16 monocytes, using Illumina Innium Human-Methylation450 BeadChips (450K arrays). Importantly, our experimental design reduces the impact of all major confounding factors in EWASs, due to the proling of puried, primary cells derived from MZ twins, who share virtually all somatic variation and early-life environmental exposure23. With the exception of one DMP in T cells, we do not identify signicant T1D-associated DMPs in any of the investigated immune cell types. However, we nd a strong enrichment of DVPs in T1D twins relative to their healthy co-twins. We also observe a cell type-specic enrichment when compared with healthy, unrelated individuals. These

Discovery set Analytical strategy

Validation sets

CD14

n

T1D-discordant twins n = 52 pairs

CD19 n

Annotation

Regulatory circuits

CD4 n

Cord blood

n = 139 n = 201 Unrelated controls

n = 12 n = 12

Temporal stability Disease onset

n = 50 n = 48

Figure 1 | Overview of the study design and analytical approach. We performed an EWAS in 52 MZ twin pairs discordant for T1D in three immune effector cell types: CD4 T cells, CD19 B cells and CD14 CD16 monocytes. We used two different approaches to determine differential DNA methylation associated with T1D status in disease-discordant twin pairs. First, we identied DMPs between T1D and healthy co-twins, which correspond to differences in mean DNA methylation levels. Second, we determined DVPs, which reect heterogeneous epigenetic outliers in T1D twins compared with their healthy co-twins. To assess the biological signicance of our ndings, we analysed three additional, genome-wide DNA methylation data sets in CD14 monocytes and CD4 T cells from 12 T1D-discordant MZ twin pairs; CD14 and CD4 cells from 201 and 139 unrelated, healthy individuals;

and cord blood from 98 newborns of whom 50 had progressed to overt T1D during childhood. Finally, we characterized T1D-associated DVPs using cell type-specic gene regulatory circuits. Credits: The immune response, Big Picture (http://bigpictureeducation.com/

Web End =http://bigpictureeducation.com/).

2 NATURE COMMUNICATIONS | 7:13555 | DOI: 10.1038/ncomms13555 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13555 ARTICLE

T1D-associated DVPs are temporally stable; not under genetic control; enriched at gene regulatory elements; and located at genes involved in immune cell metabolism and the cell cycle.

ResultsDNA methylation proles of immune effector cell types. In 52 T1D-discordant MZ twin pairs, we isolated three immune effector cell types that play a pivotal role in T1D pathobiology: CD4 T cells, CD19 B cells and CD14 CD16 monocytes1.

Cells were isolated and puried from collected peripheral blood mononuclear cells using magnetic activated cell sorting (MACS). Cell purity of each preparation was evaluated using uorescence-activated cell sorting (FACS) analysis (Supplementary Fig. 1).

For the discovery stage, we generated a total of 302 genome-wide DNA methylation proles on the 450K array platform (Fig. 1). The array platform allows the assessment of DNA methylation status at 4485,000 CpG sites at single-nucleotide resolution, and covers 99% of RefSeq genes with an average of 17 CpG sites per gene region and 96% of CpG islands24. Array data preprocessing and quality control were performed using established analytical tools (see the Methods section), leaving 406,365 CpG sites for subsequent statistical analysis (Supplementary Fig. 2a).

Multidimensional scaling and hierarchical clustering revealed that most of the variation in the data was captured by variation between twin pairs (for example, genetic effects) and cell types (Supplementary Figs 2b and 3). In addition, we performed singular value decomposition to determine principal components of variation in DNA methylation proles. In our analysis, no principal component was found to correlate with T1D status (Supplementary Fig. 2c).

Identication of T1D-associated DMPs. We rst measured differences in mean DNA methylation levels between T1D twins and their healthy co-twins in each cell type using a pair-wise analysis. We identied a single DMP at genome-wide signicance, cg01674036 in T cells (P 2.2 10 9, false

discovery rate (FDR)-corrected P 9.1 10 4, paired t test;

Fig. 2a). This DMP demonstrated a mean DNA methylation difference of 2.3% between T1D and healthy co-twins, and mapped to an intergenic region 24.3 kb downstream of the DDIT4 gene (also known as REDD1) encoding DNA-damage-inducible transcript 4 (Fig. 2b). Notably, we did not detect any additional DMPs at an FDR of o0.05 in any of the three cell types.

The DMP cg01674036 co-located with an active gene regulatory region in T cells (Fig. 2c). Chromatin interaction data obtained from a lymphoblastoid cell line provided experimental evidence that this region binds to the promoter region of DDIT4. DDIT4 functions as an inhibitor of the mammalian target of rapamycin (mTOR) complex 1; activation of mTOR complex 1 is controlled by anabolic hormones including insulin25.

The 450K array platform has a xed set of CpG sites, covering o2% of all annotated CpGs. While this platform is scalable to large sample sizes, the complementary application of sequencing-based approaches is required to comprehensively capture disease-associated DNA methylation loci on a genome-wide level11,12. To this end, we further measured DNA methylation levels in CD4

T cells using whole-genome bisulte sequencing (WGBS-seq) in four MZ twin pairs, who were originally proled on the 450K array. In total, we obtained 4500 million reads per sample resulting in a mean coverage of between 12.6 and 15.1 reads per CpG site. This allowed us to investigate over 8.7 million CpGs with a minimum coverage of 10 reads across all eight samples (Supplementary Table 1). This analysis was sufciently powered to detect differentially methylated regions (DMRs) that consist of

at least ve CpGs and exhibit a mean DNA methylation difference of 430% at an FDR of o0.05. We did not identify such DMRs to be associated with T1D, irrespective of FDR values.

In conclusion, with the exception of the DMP cg01674036, we did not identify mean DNA methylation differences between T1D twins and their healthy co-twins in any of the three immune cell types using the 450K array platform (Fig. 3a). At genomic loci not covered by the array, results based on WGBS-seq data indicate that mean DNA methylation differences of large effect size are unlikely to exist.

Identication of T1D-associated DVPs. Next, we explored whether DNA methylation variability between T1D-discordant MZ twins can shed light on the phenotypic discordance. A recent comparative study26 demonstrated that current algorithms for DVP detection can substantially differ in terms of their sensitivity and type-1 error rate (see the Methods section). Established algorithms typically assume frequent alterations in the disease phenotype, and thus lack the sensitivity to detect outlier events26. Instead, the novel algorithm iEVORA22, which is based on a regularized version of Bartletts test, improves the sensitivity to detect DVPs (see the Methods section).

Using iEVORA, we identied 10,548 DVPs in B cells, 4,314 in T cells and 6,508 in monocytes at a stringent FDR of o0.001 (Fig. 3b). Strikingly, in each cell type we found strong enrichment of DVPs that are hypervariable in T1D twins compared with their healthy co-twins (Po1 10 100, binomial test; Fig. 3c). These

T1D-associated DVPs represent epigenetic outliers that often occur in individual twin pairs and cell types (Fig. 4a and Supplementary Fig. 4). At DVPs, the DNA methylation differences between the T1D twin and its healthy co-twin were found to be comparatively large in many cases (Supplementary Fig. 4).

We next assessed a range of potential confounding factors that could lead to increased variability in DNA methylation levels, including cellular heterogeneity and differences in cell purication efciency (as quantied by FACS), age of twins at both disease diagnosis and sample collection, medication use (statins and thyroxine), as well as presence of other autoimmune diseases (thyroiditis, as characterized by thyroid peroxidase autoantibodies). We calculated the fraction of DVPs in T1D twins exhibiting a signicant deviation from the healthy co-twins, and then correlated this fraction with different potential confounding variables. For all tested variables, these correlations were not statistically signicant (P40.05;

Supplementary Fig. 5).

DNA methylation levels can associate with genetic variants in cis, leading to an increase in interindividual DNA methylation variability. Although post-zygotic somatic mutations may occur and give rise to mosaicism in identical twins (with a controlled genetic background)23, due to the heterogeneous nature of DVPs, we anticipated a modest (if any) genotypic effect on DNA methylation levels at DVPs. To nd a denitive answer, we genotyped all 52 twin pairs on Illumina HumanOmni2.58 BeadChips, and mapped methylation quantitative trait loci (meQTLs) using a linear-additive modelling approach (see the Methods section). For this analysis, only single-nucleotide polymorphism (SNPs) with a minor allele frequency of 45%

and located o50 kb up- and downstream of each CpG site were considered. We found that T1D-associated DVPs (FDRo0.001) were depleted at meQTLs compared with random sets of CpG sites in all three cell types (permutation Po1 10 4).

In summary, in all three immune cell types we discovered statistically signicant DVPs that correlate with T1D status. We provided evidence that these CpG sites are unlikely to be a

NATURE COMMUNICATIONS | 7:13555 | DOI: 10.1038/ncomms13555 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13555

a

b

10

8

1.0

cg01674036

0.8

0.6

Log 10(observed P-value)

[afii9826]-value

6

4

2

0.4

0.2

0.0

B cellsT cells Monocytes

IGR

CGI

2 4 6 8 10

Log10 (expected P-value)

cg20080702

cg26292232

cg01354384

cg15993841

cg18899730

cg02190572

cg04973995

cg01674036

cg22998476

cg14003655

cg05347948

cg07828833

c

74,034,000

74,039,000

74,044,000

74,049,000

74,054,000

74,059,000

74,064,000

74,069,000

74,074,000

74,079,000

74,084,000

74,089,000

74,094,000

74,099,00

E037 Primary T helper memory cells from peripheral blood 2 PrimaryHMM E037 Primary T helper memory cells from peripheral blood 2 ImputedHMM

E039 Primary T helper naive cells from peripheral blood PrimaryHMM E039 Primary T helper naive cells from peripheral blood ImputedHMM E040 Primary T helper memory cells from peripheral blood 1 PrimaryHMM E040 Primary T helper memory cells from peripheral blood 1 ImputedHMM

E041 Primary T helper cells PMA-I stimulated PrimaryHMM E041 Primary T helper cells PMA-I stimulated ImputedHMM E042 Primary T helper 17 cells PMA-I stimulated PrimaryHMM E042 Primary T helper 17 cells PMA-I stimulated ImputedHMM E043 Primary T helper cells from peripheral blood PrimaryHMM E043 Primary T helper cells from peripheral blood ImputedHMM E044 Primary T regulatory cells from peripheral blood PrimaryHMM E044 Primary T regulatory cells from peripheral blood ImputedHMM E045 Primary T cells effector/memory enriched from peripheral blood PrimaryHMM E045 Primary T cells effector/memory enriched from peripheral blood ImputedHMM

CpG island RefSeq genes

ChIA-PET H3K4me3 in GM12878 lymphoblastoid cells

DDIT4

Figure 2 | Assessment of the functional signicance of the T1D-associated DMP cg01674036. (a) QQ plot for the identication ofdifferentially methylated CpG positions (DMPs) between T1D-discordant MZ twin pairs in different immune effector cell types. Only the DMP cg01674036 reached genome-wide signicance in T cells, with P 2.2 10 9 (FDR-corrected P 9.1 10 4) and a mean DNA methylation difference of 2.3%. (b)

Regional plot of the locus harbouring the T-cell-specic DMP cg01674036. The statistically signicant DMP is indicated with a black arrow. Data points represent the DNA methylation b-values (y axis) at the indicated CpGs (x axis) in one individual. For each CpG site, we calculated the mean DNA methylation value (indicated with a larger data point). Every CpG site is annotated with regards to epigenomic feature and gene element using the 450K array annotation manifest. (c) Annotation of the genomic locus using epigenomic reference data sets. The genomic locus on chromosome 10q22.1 (position 74,028,00074,100,000; genome build hg19) harbouring the DMP cg01674036 (chr10:74,058,002) is shown using the WashU Epigenome

Browser v40.0.0 (http://epigenomegateway.wustl.edu/browser/

Web End =http://epigenomegateway.wustl.edu/browser/). The T1D-associated DMP is located at a CpG island (indicated with a red arrow).

A total of 16 epigenomic reference tracks provided by the Roadmap Epigenomics project are displayed. Specically, we show both the primary and imputed chromatin state maps in eight distinct primary T cell populations. The highlighted CpG island overlaps with an active transcription start site (red)or enhancer (orange/yellow) in all available T cell populations. In addition, H3K4me3 ChIA-PET data in the lymphoblastoid cell line GM12878 revealed a long-range chromatin interaction between the active regulatory element and the gene promoter region of DDIT4. CGI, CpG island; ChIA-PET, chromatin interaction analysis by paired-end tag sequencing; IGR, intergenic region.

consequence of confounding factors, including cellular heterogeneity, and that they act independently of genetic variation.

Temporal stability of T1D-associated DVPs. Following the discovery of DVPs on 450K arrays, we reassessed the

T1D-associated DNA methylation hypervariability phenotype in twins after ve years, using a second assay platform. We retrieved genome-wide DNA methylation proles of CD14 monocytes and CD4 T cells from 12 T1D-discordant MZ twin pairs generated on Illumina Innium HumanMethylation27 BeadChips

4 NATURE COMMUNICATIONS | 7:13555 | DOI: 10.1038/ncomms13555 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13555 ARTICLE

a b c

2.5

2.5

10,000

B cellsT cells Monocytes

B cellsT cells Monocytes

Healthy T1D

Hypervariable CpG sites

2.0

2.0

8,000

1.5

1.5

6,000

Density

Density

1.0

1.0

4,000

0.5

0.5

2,000

0.0

0.0

0.05 0.25 0.45 0.65 0.85

0.05 0.25 0.45 0.65 0.85

0 B cells T cells Monocytes

Pvalue Pvalue

d e

5

0.30

4

0.25

0.20

3

Odds ratio

0.15

2

0.10

1

0.05

0

0.00

MM TT TM MT TM MT

450K27K 450K450K

450K27K 450K450K

Positive predictive value

Figure 3 | DNA methylation variation in identical twin pairs discordant for T1D. (a) Histogram of P-values for the identication of differentially methylated CpG positions (DMPs) between T1D-discordant MZ twin pairs in different immune effector cell types. DMPs were determined using a paired t test. (b) Histogram of P-values for the identication of T1D-associated differentially variable CpG positions (DVPs). DVPs were determined at an FDR of o0.001 using the algorithm iEVORA. (c) Bar plots showing the enrichment of DVPs in T1D twins compared with their healthy co-twins. While this hypervariability phenotype was found in all cell types (Po1 10 100, binomial test), it was particularly pronounced in B cells. (d) Bar plots showing

the odds ratios of the assessment of temporal stability of T1D-associated DVPs in an external data set of CD14 and CD4 cells derived from12 disease-discordant MZ twin pairs generated on 27K arrays. Importantly, the identied DVPs in CD14 and CD4 cells replicated in a cell type-specic context. Stars denote statistical signicance assessed using a one-tailed Fishers exact test: *Po1 10 2 and **Po1 10 4. (e) Positive predictive values

for the analyses shown in d. B, CD19 B cells; M, CD14 CD16 monocytes; T, CD4 T cells.

(27K arrays)15. These twins belonged to the same twin registry used for the discovery cohort of this study, but provided new DNA samples for reassessment after ve years.

We conrmed directionality of the T1D-associated DVPs (FDRo0.001), indicating robust technical detection across assay platforms in both CD14 cells (P 7.7 10 5, one-tailed

Fishers exact test; Fig. 3d) and CD4 cells (P 7.8 10 3;

Fig. 3d). Consistently, the attained positive predictive values were higher in the direct cell type comparison (Fig. 3e).

Taken together, we showed that DNA methylation levels at T1D-associated DVPs are temporally stable over at least ve years in patients with established diabetes, and can be observed across two assay platforms.

Evaluation of T1D-associated DVPs in unrelated individuals. We further assessed T1D-associated DVPs using independent genome-wide DNA methylation proles retrieved from the BLUEPRINT Consortium. Specically, we reasoned that DVPs hypervariable in T1D ought to be hypervariable when compared with an external set of healthy controls. We obtained 450K array data sets of CD14 and CD4 cells derived from 201 and 139 unrelated, healthy individuals, respectively. These individuals were drawn from a population of blood donors, and thus are unlikely to have strong genetic susceptibility to T1D.

DVPs that were found to be hypervariable in T1D twins compared with their healthy co-twins, were also hypervariable

when compared with unrelated individuals with limited genetic susceptibility markers (Supplementary Fig. 6). As demonstrated before, DVPs showed cell type specicity (P 1.3 10 60 and

P 4.5 10 107, for monocytes and T cells, respectively).

In conclusion, our analysis provided further evidence that the identied DVPs represent relevant, cell type-specic markers for T1D.

Assessment of T1D-associated DVPs in cord blood. To explore whether the identied DVPs emerged before the onset of T1D, we generated genome-wide DNA methylation proles of umbilical cord blood obtained from newborns. These newborns were part of the DiPiS cohort, a population-based prospective study of T1D in children27. We selected samples from 98 newborns of whom 50 had progressed to overt T1D during childhood, while 48 did not. We hypothesized that if the T1D-associated DVPs (that are independent of genetic risk factors) were already observed in cord blood before disease onset, they could potentially contribute to T1D pathogenesis or be an early indicator of disease.

We correlated DNA methylation levels at T1D-associated DVPs identied in puried immune cell types with those in cord blood tissue. This assessment did not reach statistical signicance (P40.05, one-tailed Fishers exact test).

Based on these ndings, we conclude that the discovered DVPs occur post-birth and are likely associated with the pathogenesis of T1D either before or after the clinical diagnosis. Nonetheless,

NATURE COMMUNICATIONS | 7:13555 | DOI: 10.1038/ncomms13555 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13555

b c d e

[afii9826] > 0 in T1D [afii9826] > 0 in T1D [afii9826] < 0 in T1D [afii9826] < 0 in T1D

30

30

30

30

20

20

20

20

Enrichment (%)

Enrichment (%)

10

10

10

10

0

0

Enrichment (%)

0

Enrichment (%)

0

10

10

10

10

20

20

20

20

30

30

30

30

Island (n=4,468)

SShore (n=980)

NShore (n=1,166)

SShelf (n=83)

NShelf (n=87)

Other (n=920)

TSS200 (n=1,776)

TSS1500 (n=1,686)

5UTR (n=1,021)

1st Exon (n=684)

3UTR (n=92)

IGR (n=957) Body (n=1,488)

Other (n=5,196)

NShelf (n=781)

SShelf (n=727)

SShore (n=701)

NShore (n=952)

Island (n=908)

Body (n=4,293)

IGR (n=2,531)

3UTR (n=554)

5UTR (n=586)

1st Exon (n=137)

TSS1500 (n=894)

TSS200 (n=270)

a

Figure 4 | Functional annotation of T1D-associated DVPs. (a) Venn diagram showing the overlap of T1D-associated DVPs (FDRo0.001) across cell types. Although many of the identied DVPs were found to be cell type-specic, B cells and monocytes showed a substantial proportion of overlap.

(b,c) Enrichment of T1D-associated DVPs at different epigenomic features and gene elements. Here, only DVPs at which the DNA methylation level was increased (hypermethylated; Db40) in T1D twins compared with their healthy co-twins are shown. The enrichment is shown in relation to all 450K array probes that passed quality control. (d,e) The same analyses as shown in b and c, but for DVPs at which the DNA methylation level was reduced (hypomethylated; Dbo0) in T1D twins. (f) Integration of T1D-associated DVPs with gene regulatory circuits in CD19 B cells. The network was constructed using the corresponding genes of all T1D-associated hypomethylated DVPs that map to gene promoters and hypermethylated DVPs at gene bodies identied in B cells. The resulting network consisted of 297 genes connected via 906 regulatory edges. Three network modules were identied and are highlighted in different colours: Module 1 (n 61 genes) is shown in purple, module 2 (n 69) in green and module 3 (n 167) in orange. These

modules were further characterized using functional enrichment analysis (Supplementary Table 2 and Supplementary Table 3). IGR, intergenic region; N, north, that is, upstream; S, south, that is, downstream; TSS200/1500, 200/1500 bp upstream of a transcription start site; UTR, untranslated region.

6 NATURE COMMUNICATIONS | 7:13555 | DOI: 10.1038/ncomms13555 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13555 ARTICLE

it should be noted that cord blood, similar to peripheral whole blood, is a substantially heterogeneous tissue that hampers the precise measurement of DNA methylation levels. Our study is by far the most powered in this area, but future studies with increased statistical power will be needed to corroborate our conclusions.

Functional signicance of T1D-associated DVPs. Next, we investigated whether T1D-associated DVPs exhibit a certain level of functional organization and whether the nearby genes cluster in biological pathways. First, we performed enrichment analyses with regards to gene elements and epigenomic features as dened in the 450K array annotation manifest. We distinguished between T1D-associated DVPs (FDR o0.001) at which the DNA methylation level is either increased (hypermethylated) or decreased (hypomethylated) in T1D twins compared with their healthy co-twins. This discrimination is important, because the effect of DNA methylation on the regulation of gene expression is distinct at different gene elements5,28. Across the three immune cell types, we found an enrichment of T1D-associated hypermethylated DVPs at CpG islands (P 1.5 10 265,

hypergeometric test; Fig. 4b) and proximal gene promoters (P 2.3 10 131; Fig. 4c). In contrast, hypomethylated DVPs

were depleted at CpG islands (Po2.2 10 308; Fig. 4d) and

enriched at gene bodies (P 1.0 10 47; Fig. 4e).

We corroborated these enrichment patterns using cell type-specic chromatin state maps. We retrieved reference chromatin state data based on ve chromatin marks in primary cells from peripheral blood, provided by the Roadmap Epigenomics project29 (see the Methods section). For all three immune cell types, we found cell type-specic enrichment of T1D-associated hypermethylated DVPs at chromatin states marking active transcription start sites proximal to gene promoters (Po2.2 10 308 in B cells, P 4.2 10 38 in T cells,

and P 5.2 10 262 in monocytes, hypergeometric tests;

Supplementary Fig. 7a). We also observed depletion at states marking enhancers (P 6.2 10 269 in B cells, P 3.4 10 22

in T cells, and P 1.3 10 258 in monocytes; Supplementary

Fig. 7a). Hypomethylated DVPs showed inverted enrichment patterns (Supplementary Fig. 7b).

We then tested for enrichment of biological process ontology terms attributed to genes in proximity to T1D-associated DVPs. We adjusted for the differing number of CpGs per gene present on the 450K array to reduce bias in the gene set analysis30. Across all cell types, the T1D-associated DVPs cumulatively clustered at genes involved in molecular metabolic processes and the cell cycle (Supplementary Fig. 8). The enrichment was particularly pronounced in B cells (Supplementary Fig. 8).

In summary, these ndings showed that T1D-associated DVPs localize at specic gene regions and active chromatin states implicated in the regulation of gene expression, and highlighted gene pathways related to cell metabolism and the cell cycle.

Integration of T1D-associated DVPs with regulatory circuits. Finally, using a gene regulatory network approach we further explored the T1D-associated DVPs in B cells that may lead to transcriptional regulation of relevant pathways. We obtained CD19 B-cell-specic regulatory circuits31 that consist of interactions between transcription factors and genes derived from genome-wide promoter and enhancer activity maps presented by the FANTOM5 project32,33. We dened two sets of DVPs that may lead either to transcriptional repression or activation in B cells, and assigned these to their corresponding genes (see the Methods section). Then, we intersected the resulting gene lists with the regulatory circuits.

The regulatory network created using gene-repressing DVPs consisted of 1,465 genes and 16,712 regulatory edges. The corresponding network constructed using gene-activating DVPs consisted of 297 genes connected via 906 edges. While we were unable to resolve the structure of the gene-repressing network and could not identify signicant modules within this network, the gene-activating network showed three network modules (Fig. 4f). We further characterized these modules using gene enrichment analyses (see the Methods section), and the results of all gene set analyses are shown in detail in Supplementary Table 2.

Module 1 contained 61 genes, including NRF1 encoding nuclear respiratory factor 1 (Fig. 4f; shown in purple). NRF1 is a transcription factor that regulates the expression of genes encoding key enzymes in oxidative metabolism and mitochondrial function34. The module showed enrichment in ontology terms related to glucose-6-phosphate transmembrane transporter activity. Further analysis revealed overrepresentation of genes involved in mTOR signalling, a central pathway in the regulation of cell metabolism, growth and proliferation (Supplementary Table 3)35. Module 2 contained 69 genes (Fig. 4f; shown in green), and was enriched for genes connected to interleukin-1 receptor binding and receptor antagonist activity. This pathway is implicated in T1D-associated altered innate immunity36. Module 3 contained 167 genes, including the FOXP1 gene hub (Fig. 4f; shown in orange). FOXP1 encodes forkhead box P1, an important transcriptional regulator of B cell, T cell and monocyte differentiation. Recent studies in mice also demonstrated that Foxp1 is essential for islet a-cell proliferation and function37, and plays a key role in the regulation of systemic glucose homeostasis38.

The integration of T1D-associated DVPs with gene regulatory circuits in CD19 B cells conrmed our initial ndings (Supplementary Fig. 8), and further implicated signalling pathways related to immune cell metabolism. While alteration in these pathways could be secondary to the systemic metabolic abnormalities associated with diabetes, we note that these pathways could also predispose to autoimmune diseases including T1D.

DiscussionIn this study, we investigated whether differential epigenetic variation can explain discordance of T1D in identical twins. We measured genome-wide DNA methylation levels in 52 twin pairs across disease-relevant immune effector cell types. Our unique study design allowed us to reduce confounding factors that have impeded many previous EWASs, namely cellular heterogeneity (by using multiple, sorted, primary cell types) and genetic heterogeneity, age and early-life environmental effects (by using disease-discordant MZ twins).

Notably, with the exception of a single T-cell-specic DMP, we did not detect convincing differences in mean DNA methylation associated with T1D in our MZ twin cohort using the 450K array platform (FDRo0.05; Fig. 2a). The DMP that did reach statistical signicance, cg01674036 in CD4 T cells, is not contained on the 27K array platform and therefore could not be technically replicated in our data set. Annotation using epigenomic reference data sets revealed that the DMP maps to an active gene regulatory region in T cells and interacts with the gene promoter of DDIT4 (Fig. 2c). The corresponding protein is involved in the mTOR signalling pathway, which has been implicated in the gene network analysis of T1D-associated DVPs (Fig. 4f and Supplementary Table 3).

In addition, we did not nd DMRs of large effect size (430%; Z5 CpGs) in four disease-discordant MZ twin pairs

NATURE COMMUNICATIONS | 7:13555 | DOI: 10.1038/ncomms13555 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13555

using WGBS-seq. It is possible that T1D-associated DMPs and DMRs could be discovered in much larger cohorts or more highly selected cell populations using either Innium arrays or bisulte sequencing39. In particular, the recent availability of the Illumina Innium MethylationEPIC BeadChip covering over 850,000 CpG sites40, of which many are located at enhancer regions identied by the ENCODE and FANTOM5 projects, may allow for discovery of additional T1D-associated loci. However, if they were present, such loci would most likely be of small effect size. Indeed, this notion would be consistent with ndings from genome-wide association studies (GWASs) of T1D and other complex traits and diseases. In GWASs, many hundreds of trait-associated genetic variants have been identied, the vast majority of which possess small effect sizes41.

In contrast, we identied a substantial number of CpGs that are hypervariable in T1D twins compared with their healthy co-twins. The DNA methylation differences at DVPs were found to be comparatively large in many cases (Supplementary Fig. 4). The sensitive yet robust identication of DVPs is challenging, and is characterized by a high type-1 error rate22. Thus, replication of our ndings in independent sample cohorts in future studies is paramount.

Here, for the rst time, we detected and functionally annotated DVPs in a common disease phenotype other than cancer. DVPs have been shown to correlate with the early stages of carcinogenesis. Consistently, our data suggest that T1D-associated DVPs are associated with T1D after clinical diagnosis. In this regard, causal inference analysis may be applied to further characterize and quantify the extent of the relationship between genetic variants, epigenetic variants and phenotypic discordance14,42. Future longitudinal studies of pre-diabetic individuals will establish whether the epigenetic changes antedate the clinical diagnosis.

Our ndings have important implications for the future application of the EWAS approach to elucidate human disease mechanisms. First, the use of puried, primary cell populations likely reduced the overall number of association signals typically detected in case-control EWASs conducted in peripheral blood, but with the critical advantage of yielding genuine disease-relevant signals, if present. Second, we propose the complementary assessment of DNA methylation variability in parallel to mean DNA methylation for any future EWAS. Analytical tools are now readily available to identify DVPs for other complex traits and diseases19,22.

Our results showed that DVPs differ considerably between the proled cell types (Fig. 4a and Supplementary Fig. 4), suggesting that the response of each cell type is specic. The importance of the three immune effector cell types used in this study in the development of T1D has been recognized through experimental evidence1,3. Furthermore, genetic variants associated with T1D are enriched at enhancer sequences active in T and B cells, as well as CD34 stem cells and thymus tissue43. However, it is plausible that other cell types not assayed here may contain relevant DNA methylation differences. Alternatively, it may also be possible that rare sub-populations of the three immune cell types, such as regulatory T cells (CD25 FOXP3 cells) or T helper 17 cells (TH17 cells), harbour epigenetic signals that remain undetected in the broad population of CD4 T cells1,44,45. Future EWASs in subsets of T cells may be conducted to address this possibility.

Of note, a recent report found an increase in DNA hydroxymethylation levels at gene promoters in CD4 T cells in patients with systemic lupus erythematosus, an autoimmune disease, compared with healthy controls46. DNA hydroxymethylation remodelling has also been observed in CD4

T cell differentiation47. As our experimental approach did not allow the discrimination between methylated and

hydroxymethylated cytosine bases, DNA hydroxymethylation could thus contribute to the observed differential variability, potentially providing a general mechanism underlying the pathogenesis of autoimmune diseases.

Studies have reported the co-localization of meQTL at genetic risk loci of complex traits and common diseases identied through GWASs, including schizophrenia13, blood pressure10, and several cancer types48. Consequently, we overlapped our T1D-associated DVPs with 59 T1D genetic susceptibility loci retrieved from T1DBase, a curated web resource (http://www.t1dbase.org

Web End =http://www. http://www.t1dbase.org

Web End =t1dbase.org ; v4.19). We did not nd a statistically signicant enrichment of DVPs at these loci (P40.05, hypergeometric test).

A specic enrichment test of T1D-associated DVPs mapping to the major histocompatibility complex (MHC) locus also did not achieve statistical signicance compared with all assessed CpG sites (Supplementary Fig. 9). The MHC locus is key in conferring genetic risk of T1D and other autoimmune diseases, as it harbours many genes encoding cell surface molecules that orchestrate components of the immune system. This analysis provided further evidence that T1D-associated genetic and epigenetic variants appear to act independently.

We have identied T1D-related DVPs in immune effector cells that associate with genes involved in cell metabolism and the cell cycle (Supplementary Fig. 8 and Supplementary Table 2). Specically, by integrating T1D-associated DVPs with gene regulatory circuits in CD19 B cells, we pinpointed key transcriptional regulators such as NRF1 and FOXP1 (Fig. 4f), and pathways such as mTOR signalling (Supplementary Table 3). Indeed, the same signalling pathways have been implicated in differentiation, proliferation and metabolism of both T cells and monocytes4952. For example, deletion of Foxp1 in nave CD8

T cells leads to activation of the mTOR signalling cascade53, indicating a relationship between gene modules 1 and 3 of the regulatory network we identied here (Fig. 4f). Therefore, DVPs could modulate disease activity through the regulation of immune effector cell gene expression either before or after the induction of the disease process.

However, it remains possible that other T1D-associated DVPs result from disease-associated metabolic disturbances. Previous studies reported DMPs (including cg19693031) at the TXNIP gene to be inversely correlated with both type 2 diabetes and sustained hyperglycaemia (for example, haemoglobin A1c levels)54,55. In our data set of T1D patients, we also found cg19693031 to be a DVP in monocytes (P 9.1 10 4); this

observation suggests that a proportion of DVPs result from the diabetes-associated metabolic effect. In either case, it is likely that the impact of epigenetic changes on T1D-associated immune effector cells would adversely affect the natural history of the disease3.

The exact mechanism by which epigenetic instability in T1D is manifested, its timing in relation to induction of islet auto-immunity, as well as its impact on disease progression, remains to be explored. However, these questions can now be addressed through the study of individuals at high T1D-risk and those with variable disease severity. In this way, we can achieve our ultimate aim of identifying diagnostic and prognostic epigenetic biomarkers that can improve the management of T1D.

Methods

Ethics statement. This study was approved by the Northern and Yorkshire Research Ethics Committee (REC reference number: 06/MRE03/22) and the NRES Committee East of England-Hertfordshire (12/EE/0040). All participants gave informed consent either personally or by parental consent, as appropriate.

Study samples. MZ twin pairs were ascertained by referral through their physicians to the British Diabetic Twin Study, the Barbara Davis Center for

8 NATURE COMMUNICATIONS | 7:13555 | DOI: 10.1038/ncomms13555 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13555 ARTICLE

Childhood Diabetes and Diabetes Prevention TrialNet (USA) and the BMBF Pediatric Diabetes Biobank (Germany). T1D-associated autoantibodies were analysed by radioimmunoassay56,57. We established monozygosity by means of DNA ngerprinting using an AmpFLSTR Identiler PCR Amplication Kit (Life Technologies) and consultation of clinical data. T1D status was established by standard criteria58. T1D patients have been treated from diagnosis with insulin and take highly puried human insulin at least twice daily. We excluded twins who were pregnant and twins with signicant co-morbidities including severe macrovascular and microvascular complications of diabetes. Umbilical cord blood was obtained from 35,000 newborns enroled in the DiPiS Study (Sweden)27 between the years 2000 and 2004. From this cohort, we selected 98 neonates of whom 50 progressed to T1D and 48 did not. Children were followed for 15 years to monitor if they develop markers of islet autoimmunity and T1D. The samples consisted of dried cord blood spots dotted onto cards.

Cell sorting and purity analysis. Peripheral blood mononuclear cells were prepared from 50 ml of heparinized blood using Percoll density gradient separation. CD4 T cells, CD19 B cells and CD14 CD16 monocytes were isolated using MACS according to the manufacturers instruction. First, CD19

B cells were separated with CD19 MicroBeads (130-050-301, Miltenyi Biotec). The negative fraction was then washed and incubated with CD16 MicroBeads (130-045-701, Miltenyi Biotech). The fraction depleted of CD16 cells was selected for CD14 monocytes using CD14 MicroBeads (130-050-201,

Miltenyi Biotech). Finally, the resulting negative fraction was further incubated with CD4 MicroBeads (130-045-101, Miltenyi Biotech) to obtain CD4 T cells.

Based on the number of isolated peripheral blood mononuclear cells, we used 50 ml of CD19 MicroBeads, 20 ml of CD14 MicroBeads, and 20 ml of CD4 MicroBeads per 10 million total cells. We assessed the puried cell populations with FACS. The following antibodies were used at a dilution of 1:11 per 10 million total cells for each cell type: 20 ml of FITC-conjugated mouse anti-human CD14 clone MjP9 (345784, BD Biosciences) and 10 ml of CD4 clone M-T466 (130-080-501, Miltenyi Biotech); 10 ml of phycoerythrin (PE)-conjugated mouse anti-human CD19 clone LT19 (130-091-247, Miltenyi Biotech) and 20 ml of CD16 clone B73.1/leu11c (332779, BD Biosciences); 5 ml of PerCP-Cy5.5-conjugated mouse anti-human CD64 clone 10.1 (561194, BD Biosciences); and 5 ml of PE-Cy7-conjugated mouse anti-human CD45 clone HI30 (MHCD4512, Invitrogen). Cells were incubated with antibodies at 4 C for 15 min, washed with 2 ml of phosphate-buffered saline (PBS) and ethylenediaminetetraacetic acid (EDTA), and re-suspended in a volume of 500 ml for FACS analysis. Across all cell types, the mean cell purity was 90%.

DNA extraction. We extracted genomic DNA from MACS-enriched cell populations using a QIAamp DNA Blood Mini Kit (QIAGEN) according to manufacturers instructions. DNA was extracted from cord blood using a GenSolve DNA Recovery Kit (Labtech) according to the manufacturers instructions. DNA concentration was determined using a Qubit dsDNA HS Assay Kit (Invitrogen) and DNA integrity visually inspected on a 2% agarose gel.

Illumina Innium HumanMethylation450 assay. Genomic DNA was bisulte-converted using an EZ-96 DNA Methylation MagPrep Kit (Zymo Research) according to the manufacturers instructions. We applied 500 ng of genomic DNA to bisulte treatment, and eluted puried, bisulte-converted DNA in 20 ml of

M-Elution Buffer (Zymo Research). DNA methylation levels were measured on Innium HumanMethylation450 BeadChips (Illumina) following the manufacturers protocol. In brief, 4 ml of bisulte-converted DNA was isothermally amplied, enzymatically fragmented and precipitated. Next, precipitated DNA was resuspended in hybridization buffer and dispensed onto the BeadChips. To limit batch effects, samples were randomly distributed across slides and arrays. The hybridization was performed at 48 C for 20 h using a Hybridization Oven (Illumina). After hybridization, BeadChips were washed and processed through a single-nucleotide extension followed by immunohistochemistry staining using a Freedom EVO robot (Tecan). Finally, the BeadChips were imaged using an iScan Microarray Scanner (Illumina).

Illumina Innium HumanMethylation450 data preprocessing. The DNA methylation fraction at a specic CpG site was calculated as b M (M U 100) 1,

for which M and U denote methylated and unmethylated uorescent signal intensities, respectively. The b-value statistic ranges from absent (b 0) to

complete DNA methylation (b 1) at a specic CpG. We normalized the

450K array data using BMIQ (Beta MIxture Quantile dilation), an intra-array normalization method that adjusts the b-values of type-2 design probes into a statistical distribution characteristic of type-1 probes59. Next, we ltered (1) probes with median detection P-value Z0.01 in one or more samples; (2) probes with bead count of o3 in at least 5% of samples; (3) probes mapping to sex chromosomes; (4) non-CG probes; (5) probes mapping to ambiguous genomic locations60; and (6) probes harbouring annotated SNPs within 2 bp of the probed CG irrespective of allele frequency in the European populations, as reported by dbSNP v135 (ref. 60). Finally, we adjusted for known batch effects using an empirical Bayesian framework61, as implemented in the ComBat function of the R package SVA62. The nal data matrix consisted of b-values across 406,365 CpG

sites 302 samples, that is, 49, 50 and 52 MZ twin pairs in T cells, B cells and

monocytes, respectively.

Identication of DMPs and DVPs. To identify DMPs, we applied a pairedt test and estimated the FDR using the R package q-value63. DVPs were identied using iEVORA22, an algorithm based on a regularized version of Bartlettstest. The algorithm is freely available as an executable R script from the Supplementary Information of the publication at http://www.nature.com/ncomms/

Web End =http://www.nature.com/ncomms/ . A disadvantage of Bartletts test is that single outliers can drive the DVP ranking. Therefore, iEVORA uses a novel procedure to regularize Bartletts test, by selecting CpGs based on signicant Bartletts test P-values, but ranking these selected features according to t test P-values22. This heuristic method guarantees (1) that selected CpGs are signicant DVPs; and (2) that the ranking favours DVPs that are either DMPs at genome-wide signicance or as close to being DMPs as possible. This regularization step favours DVPs that are driven by more frequent outliers compared with DVPs driven by single outliers. Bartletts test P-values from iEVORA are corrected for multiple testing using the FDR method implemented in the R package q-value. To keep the number of false positives as small as possible, avoiding any impact on the top-ranked features, a stringent FDR of o0.001 was used. Of note, application of an alternative approach, DiffVar (ref. 64), did not reveal signicant DVPs at an FDR of o0.05. DiffVar compares the absolute deviations from the respective group means using a (moderated) t test, as the method assumes that the differential variability is driven by numerous outliers within a disease phenotype. This algorithm offers improved control of the type-1 error rate at the expense of reduced power26. Thus, iEVORA can be seen as a compromise between DiffVar (which ignores differential variability driven by few outliers resulting in a low type-1 error rate and low sensitivity) and EVORA19 (which favours differential variability driven by single outliers resulting in much greater sensitivity albeit at the expense of a higher type-1 error rate).

WGBS-seq data preprocessing and DMR calling. Sample preparation and preprocessing of WGBS-seq data were conducted using previously established protocols and pipelines65. Sequencing statistics are provided in Supplementary Table 1. Counts of unmethylated and methylated cytosine in the context of CpG sites were extracted from the mapped BAM les using a publicly available algorithm (https://bitbucket.org/lowelabqmul/bs-seq-dmr-caller

Web End =https://bitbucket.org/lowelabqmul/bs-seq-dmr-caller). In brief, the algorithm uses a windowless approach that progresses along the genome and determines groups of CpG sites that have the same directional difference between cases and controls. The method requires each CpG to be located within 1,000 bp of its neighbouring CpG. To determine the signicance of the DMR, the w2 statistic is calculated for the pooled counts across each of the CpGs at the locus and across all the samples. The sample identities are then permutated and a new statistic is calculated; this is repeated 1,000 times, and the original statistic is compared with the permutated statistics to produce a P-value for each CpG. Then, the P-value for each CpG is combined into a single P-value using Fishers method. Finally, the FDR is estimated for each of the DMRs using the R package q-value63.

Whole-genome genotyping and meQTL mapping. The quantity and integrity of DNA samples were assessed using a NanoDrop spectrophotometer (Thermo Scientic). Samples were normalized to a concentration of 50 ng ml 1 before amplication. Then, DNA was hybridized to Innium HumanOmni2.58 v1.2

BeadChips (Illumina), according to the manufacturers instructions. Following genotyping, raw data were imported into GenomeStudio (Illumina), and genotypes called using the standard cluster le provided by the arrays. Quality checks, including comparisons with called versus reported sex and genotype consistency between twins, were performed using GenomeStudio. We excluded all SNPs with a minor allele frequency of o5% and HardyWeinberg equilibrium o1 10 6,

leaving 609,587 SNPs for subsequent meQTL analysis. Further, we conrmed matching DNA methylation and genotype data sets by comparison of genotype calls across the Innium platforms. To investigate whether DNA methylation levels at DVPs are correlated with genotypes, we mapped meQTLs genome-wideusing the software Matrix eQTL66. We applied standard parameters except the P-value output threshold was set to 1 10 8 and the maximum distance between

interactions of CpGs and SNPs was set to 100,000 bp. We included the following covariates in the linear-additive model: age, sex, batch and T1D status. The analysis identied 13,579 CpG sites for T cells, 11,790 for B cells and 15,531 for monocytes that correlated with at least one SNP. Then, we determined whether T1D-associated DVPs are enriched at meQTLs compared with random sets of CpGs (n 10,000).

Assessment of DVPs in additional data sets. We retrieved DNA methylation proles of CD14 monocytes and CD4 T cells from 12 T1D-discordant MZ twin pairs generated using 27K arrays15. In addition, we used 450K array DNA methylation proles of CD14 and CD4 cells from 201 and 139 unrelated, healthy individuals, respectively, obtained from the BLUEPRINT Consortium.

From the DVPs identied using the 450K array in the discovery stage (FDRo0.001), we selected all probes that were also present in the external data set.

First, we computed the log-ratio of the variances in T1D twins versus healthy co-twins. To assess congruence between the discovery and validation sets, we then

NATURE COMMUNICATIONS | 7:13555 | DOI: 10.1038/ncomms13555 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 9

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13555

calculated the log-ratio of the variances in each set against each other. Finally, we counted the number of selected DVPs with signicant P-values in the external set and the subset of those that were hypervariable and hypovariable in T1D cases. This resulted in a 2 2 table, with a subsequent Fishers exact test allowing us to

statistically assess whether the selected DVPs validate in the external set.

Functional annotation of T1D-associated DVPs. For the enrichment analyses with regards to gene elements and epigenomic features, we used the annotation provided by the 450K array annotation manifest. For the analyses with regards to chromatin states, we retrieved data generated using the core 15-state ChromHMM model based on ve chromatin marks (H3K4me3, H3K4me1, H3K36me3, H3K27me3 and H3K9me3) from http://egg2.wustl.edu/roadmap/web_portal/

Web End =http://egg2.wustl.edu/roadmap/web_portal/ . We selected chromatin states maps of primary B cells, T cells and monocytes from peripheral blood, corresponding to the reference epigenome identiers E032, E034 and E029, respectively. Enrichment was assessed by repeated random sampling (n 1,000) using all probes that passed quality control. T1D-associated DVPs were

linked to genes using the 450K array annotation manifest24. Then, by applying the function gometh implemented in the R package missMethyl30, genes were associated with ontology terms and enrichment of these terms was calculated in relation to all CpG sites on the 450K array platform that passed quality control. This method takes account of the differing number of probes per gene present on the 450K array.

Analysis of CD19 B-cell-specic regulatory circuits. We retrieved the CD19 B cell regulatory network from http://regulatorycircuits.org

Web End =http://regulatorycircuits.org . The network consists of 11,997 nodes (genes) and 1,148,319 edges (interactions between transcription factors and regulatory elements of target genes). We selected all T1D-associated DVPs that were genome-wide signicant (P 1.2 10 7).

Gene-activating DVPs were dened as CpG sites that were either hypomethylated in T1D twins compared with their healthy co-twins and annotated as TSS1500, TSS200, 50-UTR or 1stExon on the 450K array annotation manifest; or hyper-methylated and annotated as Body or 30-UTR. Accordingly, gene-repressing DVPs were dened as CpGs that showed hypermethylation at gene promoters or hypomethylation at gene bodies. We only considered genes that directly interacted with other genes of the dened gene set. Network modules were identied using Gephi (http://gephi.org

Web End =http://gephi.org) and the Lovain method67. Then, we performed functional enrichment analyses of these modules using the R packages GOstats68 and ReactomePA69. We tested for overrepresentation of gene ontology (GO) molecular function terms using the following parameters: conditional TRUE and

FDRo0.25 (Benjamini and Hochberg method70). We performed further functional enrichment tests of the network modules at an FDR of o0.01 using Cytoscape71 and ClueGO72. For these tests, we specied the following ontologies: GO Biological

Process, GO Immune System Process, GO Molecular Function, KEGG, REACTOME and WikiPathways. We applied GO Term Fusion and a minimum number of three genes or 4% of all genes for the corresponding GO category or pathway. The kappa score was set to 0.4. All enrichment analyses of network modules were contrasted to all genes in the whole regulatory network that were also associated with CpG sites passing quality control on the 450K array platform (n 10,660).

Data availability. All 450K array and WGBS-seq data sets that support the ndings of this study have been deposited in the European Genome-phenome Archive (EGA) with the accession code EGAS00001001598 (https://www.ebi.ac.uk/ega/studies/EGAS00001001598

Web End =https://www.ebi.ac.uk/ https://www.ebi.ac.uk/ega/studies/EGAS00001001598

Web End =ega/studies/EGAS00001001598 ). We retrieved 450K array data sets of CD14 and CD4 cells from EGA with the accession code EGAS00001001456 (https://www.ebi.ac.uk/ega/studies/EGAS00001001456

Web End =https://

https://www.ebi.ac.uk/ega/studies/EGAS00001001456

Web End =www.ebi.ac.uk/ega/studies/EGAS00001001456 ), and 27K array data sets of CD14 and CD4 cells from Gene Expression Omnibus (GEO) with the accession code

GSE56606 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE56606

Web End =http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE56606).

References

1. Bluestone, J. A., Herold, K. & Eisenbarth, G. Genetics, pathogenesis and clinical interventions in type 1 diabetes. Nature 464, 12931300 (2010).

2. Patterson, C. C., Dahlquist, G. G., Gyrs, E., Green, A. & Soltsz, G. Incidence trends for childhood type 1 diabetes in Europe during 1989-2003 and predicted new cases 2005-20: a multicentre prospective registration study. Lancet 373, 20272033 (2009).

3. Atkinson, M. A., Eisenbarth, G. S. & Michels, A. W. Type 1 diabetes. Lancet 383, 6982 (2014).

4. Redondo, M. J. et al. Heterogeneity of type I diabetes: analysis of monozygotic twins in Great Britain and the United States. Diabetologia 44, 354362 (2001).

5. Schbeler, D. Function and information content of DNA methylation. Nature 517, 321326 (2015).

6. Bird, A. DNA methylation patterns and epigenetic memory. Genes Dev. 16, 621 (2002).

7. Michels, K. B. et al. Recommendations for the design and analysis of epigenome-wide association studies. Nat. Methods 10, 949955 (2013).

8. Rakyan, V. K., Down, T. A., Balding, D. J. & Beck, S. Epigenome-wide association studies for common human diseases. Nat. Rev. Genet. 12, 529541 (2011).

9. Paul, D. S. & Beck, S. Advances in epigenome-wide association studies for common diseases. Trends Mol. Med. 20, 541543 (2014).

10. Kato, N. et al. Trans-ancestry genome-wide association study identies 12 genetic loci inuencing blood pressure and implicates a role for DNA methylation. Nat. Genet. 47, 12821293 (2015).

11. Allum, F. et al. Characterization of functional methylomes by next-generation capture sequencing identies novel disease-associated variants. Nat. Commun. 6, 7211 (2015).

12. Bell, J. T. et al. Differential methylation of the TRPA1 promoter in pain sensitivity. Nat. Commun. 5, 2978 (2014).

13. Hannon, E. et al. Methylation QTLs in the developing brain and their enrichment in schizophrenia risk loci. Nat. Neurosci. 19, 4854 (2016).

14. Liu, Y. et al. Epigenome-wide association data implicate DNA methylation as an intermediary of genetic risk in rheumatoid arthritis. Nat. Biotechnol. 31, 142147 (2013).

15. Rakyan, V. K. et al. Identication of type 1 diabetes-associated DNA methylation variable positions that precede disease diagnosis. PLoS Genet. 7, e1002300 (2011).

16. Stefan, M., Zhang, W., Concepcion, E., Yi, Z. & Tomer, Y. DNA methylation proles in type 1 diabetes twins point to strong epigenetic effects on etiology.J. Autoimmun. 50, 3337 (2014).17. Issa, J.-P. Epigenetic variation and cellular Darwinism. Nat. Genet. 43, 724726 (2011).

18. Hansen, K. D. et al. Increased methylation variation in epigenetic domains across cancer types. Nat. Genet. 43, 768775 (2011).

19. Teschendorff, A. E. & Widschwendter, M. Differential variability improves the identication of cancer risk markers in DNA methylation studies proling precursor cancer lesions. Bioinformatics 28, 14871494 (2012).

20. Teschendorff, A. E. et al. Epigenetic variability in cells of normal cytology is associated with the risk of future morphological transformation. Genome Med. 4, 24 (2012).

21. Teschendorff, A. E. et al. The dynamics of DNA methylation covariation patterns in carcinogenesis. PLoS Comput. Biol. 10, e1003709 (2014).

22. Teschendorff, A. E. et al. DNA methylation outliers in normal breast tissue identify eld defects that are enriched in cancer. Nat. Commun. 7, 10478 (2016).

23. Castillo-Fernandez, J. E., Spector, T. D. & Bell, J. T. Epigenetics of discordant monozygotic twins: implications for disease. Genome Med. 6, 60 (2014).

24. Bibikova, M. et al. High density DNA methylation array with single CpG site resolution. Genomics 98, 288295 (2011).

25. Dennis, M. D., McGhee, N. K., Jefferson, L. S. & Kimball, S. R. Regulated in DNA damage and development 1 (REDD1) promotes cell survival during serum deprivation by sustaining repression of signaling through the mechanistic target of rapamycin in complex 1 (mTORC1). Cell Signal. 25, 27092716 (2013).

26. Teschendorff, A. E., Jones, A. & Widschwendter, M. Stochastic epigenetic outliers can dene eld defects in cancer. BMC Bioinformatics 17, 178 (2016).

27. Larsson, K. et al. Genetic and perinatal factors as risk for childhood type 1 diabetes. Diabetes Metab. Res. Rev. 20, 429437 (2004).

28. Jones, P. A. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat. Rev. Genet. 13, 484492 (2012).

29. Roadmap Epigenomics Consortium. Integrative analysis of 111 reference human epigenomes. Nature 518, 317330 (2015).

30. Phipson, B., Maksimovic, J. & Oshlack, A. missMethyl: an R package for analysing methylation data from Illuminas HumanMethylation450 platform. Bioinformatics 32, 286288 (2016).

31. Marbach, D. et al. Tissue-specic regulatory circuits reveal variable modular perturbations across complex diseases. Nat. Methods 13, 366370 (2016).

32. The FANTOM Consortium and the RIKEN PMI and CLST (DGT).

A promoter-level mammalian expression atlas. Nature 507, 462470 (2014).

33. Andersson, R. et al. An atlas of active enhancers across human cell types and tissues. Nature 507, 455461 (2014).

34. Patti, M. E. et al. Coordinated reduction of genes of oxidative metabolism in humans with insulin resistance and diabetes: potential role of PGC1 and NRF1. Proc. Natl Acad. Sci. USA 100, 84668471 (2003).

35. Zoncu, R., Efeyan, A. & Sabatini, D. M. mTOR: from growth signal integration to cancer, diabetes and ageing. Nat. Rev. Mol. Cell Biol. 12, 2135 (2011).36. Moran, A. et al. Interleukin-1 antagonism in type 1 diabetes of recent onset: two multicentre, randomised, double-blind, placebo-controlled trials. Lancet 381, 19051915 (2013).

37. Spaeth, J. M. et al. The FOXP1, FOXP2 and FOXP4 transcription factors are required for islet alpha cell proliferation and function in mice. Diabetologia 58, 18361844 (2015).

38. Zou, Y. et al. Forkhead Box P1 (FOXP1) transcription factor regulates hepatic glucose homeostasis. J. Biol. Chem. 290, 3060730615 (2015).

10 NATURE COMMUNICATIONS | 7:13555 | DOI: 10.1038/ncomms13555 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13555 ARTICLE

39. Tsai, P.-C. & Bell, J. T. Power and sample size estimation for epigenome-wide association scans to detect differential DNA methylation. Int. J. Epidemiol. 44, 14291441 (2015).

40. Moran, S., Arribas, C. & Esteller, M. Validation of a DNA methylation microarray for 850,000 CpG sites of the human genome enriched in enhancer sequences. Epigenomics 8, 389399 (2016).

41. Visscher, P. M., Brown, M. A., McCarthy, M. I. & Yang, J. Five years of GWAS discovery. Am. J. Hum. Genet. 90, 724 (2012).

42. Olsson, A. H. et al. Genome-wide associations between genetic and epigenetic variation inuence mRNA expression and insulin secretion in human pancreatic islets. PLoS Genet. 10, e1004735 (2014).

43. Onengut-Gumuscu, S. et al. Fine mapping of type 1 diabetes susceptibility loci and evidence for colocalization of causal variants with lymphoid gene enhancers. Nat. Genet. 47, 381386 (2015).

44. Herold, K. C., Vignali, D. A. A., Cooke, A. & Bluestone, J. A. Type 1 diabetes: translating mechanistic observations into effective clinical outcomes. Nat. Rev. Immunol. 13, 243256 (2013).

45. Bluestone, J. A. et al. Type 1 diabetes immunotherapy using polyclonal regulatory T cells. Sci. Transl. Med. 7, 315ra189 (2015).

46. Zhao, M. et al. Increased 5-hydroxymethylcytosine in CD4 T cells in

systemic lupus erythematosus. J. Autoimmun. 69, 6473 (2016).47. Nestor, C. E. et al. 5-Hydroxymethylcytosine remodeling precedes lineage specication during differentiation of human CD4 T cells. Cell Rep. 16,

559570 (2016).48. Heyn, H. et al. Linkage of DNA methylation quantitative trait loci to human cancer risk. Cell Rep. 7, 331338 (2014).

49. Feng, X. et al. Foxp1 is an essential transcriptional regulator for the generation of quiescent naive T cells during thymocyte development. Blood 115, 510518 (2010).

50. Araki, K. et al. mTOR regulates memory CD8 T-cell differentiation. Nature 460, 108112 (2009).

51. Cheng, S. C. et al. mTOR- and HIF-1alpha-mediated aerobic glycolysis as metabolic basis for trained immunity. Science 345, 1250684 (2014).

52. Ray, J. P. et al. The interleukin-2-mTORc1 kinase axis denes the signaling, differentiation, and metabolism of T helper 1 and follicular B helper T cells. Immunity 43, 690702 (2015).

53. Wei, H. et al. Foxp1 controls naive CD8 T cell quiescence by simultaneously

repressing key pathways in cellular metabolism and cell cycle progression.J. Immunol. 196, 35373541 (2016).54. Chambers, J. C. et al. Epigenome-wide association of DNA methylation markers in peripheral blood from Indian Asians and Europeans with incident type 2 diabetes: a nested case-control study. Lancet Diabetes Endocrinol. 3, 526534 (2015).

55. Soriano-Trraga, C. et al. Epigenome-wide association study identies TXNIP gene associated with type 2 diabetes mellitus and sustained hyperglycemia. Hum. Mol. Genet. 25, 609619 (2016).

56. Hawa, M. I. et al. Antibodies to IA-2 and GAD65 in type 1 and type 2 diabetes: isotype restriction and polyclonality. Diabetes Care 23, 228233 (2000).

57. Wenzlau, J. M. et al. The cation efux transporter ZnT8 (Slc30A8) is a major autoantigen in human type 1 diabetes. Proc. Natl Acad. Sci. USA 104, 1704017045 (2007).

58. Expert Committee on the Diagnosis and Classication of Diabetes Mellitus. Report of the expert committee on the diagnosis and classication of diabetes mellitus. Diabetes Care 26, S5S20 (2003).

59. Teschendorff, A. E. et al. A beta-mixture quantile normalization method for correcting probe design bias in Illumina Innium 450k DNA methylation data. Bioinformatics 29, 189196 (2013).

60. Nordlund, J. et al. Genome-wide signatures of differential DNA methylation in pediatric acute lymphoblastic leukemia. Genome Biol. 14, R105 (2013).

61. Johnson, W. E. & Li, C. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 8, 118127 (2007).

62. Leek, J. T., Johnson, W. E., Parker, H. S., Jaffe, A. E. & Storey, J. D. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics 28, 882883 (2012).

63. Storey, J. D. & Tibshirani, R. Statistical signicance for genomewide studies. Proc. Natl Acad. Sci. USA 100, 94409445 (2003).

64. Phipson, B. & Oshlack, A. DiffVar: a new method for detecting differential variability with application to methylation in cancer and aging. Genome Biol. 15, 465 (2014).

65. Kulis, M. et al. Whole-genome ngerprint of the DNA methylome during human B cell differentiation. Nat. Genet. 47, 746756 (2015).

66. Shabalin, A. A. Matrix eQTL: ultra fast eQTL analysis via large matrix operations. Bioinformatics 28, 13531358 (2012).

67. Blondel, V. D., Guillaume, J.-L., Lambiotte, R. & Lefebvre, E. Fast unfolding of communities in large networks. J. Stat. Mech. 10, P10008 (2008).

68. Falcon, S. & Gentleman, R. Using GOstats to test gene lists for GO term association. Bioinformatics 23, 257258 (2007).

69. Yu, G. & He, Q.-Y. ReactomePA: an R/Bioconductor package for reactome pathway analysis and visualization. Mol. BioSyst. 12, 477479 (2016).

70. Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B 57, 289300 (1995).

71. Shannon, P. et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 24982504 (2003).

72. Bindea, G. et al. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics 25, 10911093 (2009).

Acknowledgements

This work was funded by the EU-FP7 project BLUEPRINT (282510) and the Wellcome Trust (99148). We thank all twins for taking part in this study; Kerra Pearce and Mark Kristiansen (UCL Genomics) for processing the Illumina Innium HumanMethylation450 BeadChips; Rasmus Bennet for technical assistance; and Laura Phipps for proofreading the manuscript. The BMBF Pediatric Diabetes Biobank recruits patients from the National Diabetes Patient Documentation System (DPV), and is nanced by the German Ministry of Education and Research within the German Competence Net Diabetes Mellitus (01GI1106 and 01GI1109B). It was integrated into the German Center for Diabetes Research in January 2015. We thank the Swedish Research Council and SUS Funds for support. We gratefully acknowledge the participation of all NIHR Cambridge BioResource volunteers, and thank the Cambridge BioResource staff for their help with volunteer recruitment. We thank members of the Cambridge BioResource SAB and Management Committee for their support of our study and the NIHR Cambridge Biomedical Research Centre for funding. The Cardiovascular Epidemiology Unit is supported by the UK Medical Research Council (G0800270), BHF (SP/09/002),and NIHR Cambridge Biomedical Research Centre. Research in the Ouwehand laboratory is supported by the NIHR, BHF (PG-0310-1002 and RG/09/12/28096) and NHS Blood and Transplant. K.D. is funded as a HSST trainee by NHS Health Education England. M.F. is supported by the BHF Cambridge Centre of Excellence (RE/13/6/ 30180). A.D., E.L., L.C. and P.F. receive additional support from the European Molecular Biology Laboratory. A.K.S. is supported by an ADA Career Development Award (1-14-CD-17). B.O.B. and R.D.L. acknowledge support from the Deutsche Forschungsgemeinschaft (DFG) and European Federation for the Study of Diabetes, respectively.

Author contributions

R.D.L., S.B., V.K.R. and D.S.P. designed and supervised the study. M.A.N.D.and D.S.P. performed the experiments. D.S.P, A.E.T., R.L. and S.E. analyseddata and performed statistical analyses. All other authors provided samples or analysis tools. D.S.P. wrote the manuscript. All authors read and approved the nal manuscript.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: Paul Flicek is a member of the Scientic Advisory Board for Omicia, Inc. All other authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Paul, D. S. et al. Increased DNA methylation variability in type 1 diabetes across three immune effector cell types. Nat. Commun. 7, 13555 doi: 10.1038/ncomms13555 (2016).

Publishers note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional afliations.

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the articles Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

Web End =http://creativecommons.org/licenses/by/4.0/

r The Author(s) 2016

NATURE COMMUNICATIONS | 7:13555 | DOI: 10.1038/ncomms13555 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 11