A peripheral epigenetic signature of immune system genes is linked to neocortical thickness and memory

Increasing age is tightly linked to decreased thickness of the human neocortex. The biological mechanisms that mediate this effect are hitherto unknown. The DNA methylome, as part of the epigenome, contributes significantly to age-related phenotypic changes. Here, we identify an epigenetic signature that is associated with cortical thickness (P=3.86 × 10−8) and memory performance in 533 healthy young adults. The epigenetic effect on cortical thickness was replicated in a sample comprising 596 participants with major depressive disorder and healthy controls. The epigenetic signature mediates partially the effect of age on cortical thickness (P<0.001). A multilocus genetic score reflecting genetic variability of this signature is associated with memory performance (P=0.0003) in 3,346 young and elderly healthy adults. The genomic location of the contributing methylation sites points to the involvement of specific immune system genes. The decomposition of blood methylome-wide patterns bears considerable potential for the study of brain-related traits.

The authors then test whether the ICA associates with genetic variation using GSEA. Here two GO terms are enriched. The weighted scores from the genetic associations associates with memory performance in 5 samples. It is not clear whether this is not also a circular approach, as a number of papers describe that variance in DNA methylation is oftern carried by underlying genetic variation. Would a combined score from genetic and epigenome analysis perform even better? Are the SNPs enriched in the GSEA within mQTLs for the CPGs in the ICA2?
Overall this is an interesting paper with some robust and replicated findings. Some questions regarding the robustness of the ICA remain open and the relevance of the genetic analysis. This reviewer is not a statistics expert and would defer comments regarding ICA to reviewers more expert in the field.
Reviewer #2 (Remarks to the Author): This study tests to what degree methylomic profiles can account for age-related cortical thickness decrease. They found a significant contribution from epigenetic marks, which further point to genes involved in immune system regulation. There were also tendencies that the cortical region where thickness was related to this marker correlated with episodic memory performance. I think this is a well-conducted and interesting study, which adds to our understanding of the very pervasive phenomenon of cortical thinning with age. The study is original and will be of great interest to the research community.
I have a few questions and suggestions, especially regarding some of the analyses presented and how the sample(s) are described.
ICV and thickness are usually not highly correlated. To remove global correlations prior to EFA, maybe adjustment for mean thickness would be more appropriate?
The best strategy to handle the common influence of age, while not throwing the baby out with the water, is challenging. I think the mediation analysis used by the authors is a way of showing the mutual relationship between the variables in a clear way. In the EFA analyses, however, age was regressed out prior to the factor extraction. That may be ok, but it was not entirely clear to me how age was handled in the following analyses. I understand that ICA2 patterns were used, but were the age-adjusted thickness estimates also used, or rather the raw thickness estimates based on the ageadjusted factors? If the latter is the case, maybe this could be followed up by a mediation analysis where age is included, since it is a major source of variation that is now hidden. As I read the rational for the study, understanding age-related changes in cortical thickness is a major objective, and it would thus be informative if age was not just taken out of the equation. Similar concerns regard the memory analyses.
It is mentioned that synaptic pruning during development may account for cortical thinning. It could also be mentioned that myelination of lower cortical layers could cause the cortical mantle to appear thinner on the MR scans. This likely accounts for a substantial part of the observed cortical thinning in development.
Could the presentation of which samples that were used for the different analyses be more clearly presented, e.g. in a table, figure or flow chart? If I understand correctly, the age-thickness-methylome associations were found in a sample of you participants only? This should be more clearly stated in the ms, with age-ranges also. One easily gets the impression that this is an older-age study also, and I think this needs to be much more directly addressed, and also touched upon in the discussion, because it will have substantial impact on how the results are conceived.
Minor point: Is it possible to include some for quantitative information in the abstract also, such as sample size and effect sizes? The Fischl et al paper on whole-brain segmentation would not be the most appropriate as the first FreeSurfer-reference here, since the present paper is focused on cortical thickness and the Fischl paper on volume. The 1999 and 2000-papers from Dale and Fischl should rather be referenced first instead.
Reviewer #3 (Remarks to the Author): Freytag and colleagues analysed whole blood DNA methylation profiles generated using the 450k array in ~500 individuals and correlated these (primarily) to cortical thickness. The authors employed a wide collection of statistical techniques including state-of-the-art approaches. However, the analysis largely skips accepted straight-forward approaches for the analysis of DNA methylation data in population studies and instead heavily relies on 'black-box' types of statistical analysis, more suitable for second order analyses. It therefore is difficult to assess what data go into the analysis and even more so to interpret the outcome. It therefore remains questionable whether this study, as the authors claim, shows that the decomposition of blood methylomes bears considerable potential for the study of brain-related traits.
-The straight-forward analysis that should be presented is an epigenome-wide association analysis (EWAS) testing for association between the methylation of individual CpGs and cortical thickness. This approach allows the reader to gauge the relationship between the 'raw' data and the outcomes. -The authors instead use independent component analysis (ICA) and test the resulting components, without explaining why they jumped to such a global analysis and why specifically ICA. It seems ICA is an alternative to the more common principle component analysis (PCA) when the requirements for PCA are not met: the data are non-Gaussian and (very) noisy. The authors do not state why they used ICA but show data that indicate that the data are indeed quite noisy which raises doubts about the validity of the data sets (an may preclude an EWAS-type analysis): (1) ICA decomposition resulted in 126 independent component, 111 (88%) of which were driven by a single individual. This suggests that 111 of the 533 individuals actually represent outliers perhaps due to technical issues and probably should have been omitted from the analysis. (2) SVA identified 40 surrogate variables which is very large considering the sample size and indicates significant technical variability.
-Since the primary research question relates to the association of whole blood DNA methylation profiles with cortical thickness, SVA should be performed using cortical thickness as outcome, not age. Now, the interpretation of associations of ICA components remains unclear: is it driven by age, an age-effect on cell composition etc? -If the primary research question was to gain insight in the role of age in the decline of cortical thickness (as perhaps suggested in the introduction), a specific analysis of known age-associated CpGs (or scores based on multiple CpGs) in blood would have been in order (as mentioned in introduction).
-The 15 ICAs not driven by a single individual are not described. How much of the variability do they explain? What is their nature: with which (technical, cell count variable etc) variables do they correlate? It is elemental to know which CpGs contributed to the ICAs to contribute to our biological understanding. Are they known to be associated with age or differentially methylated between fine grained cell counts? -The main PCAs in whole blood 450k data usually capture cell counts (beyond the larger classes like the 3 measured here). It is not unlikely that the ICAs also represent cell counts. Hence, if true, the data may suggest an association between the immune system (measured through DNA methylation instead of direct cell counts) and cortical thickness. If this association survives robust control for confounding (i.e. including parameters known to be associated with cortical thickness (and immune system traits) (social economic position, medication use, etc.)). This in fact may be a very interesting observation but requires much more specific insight in the biology of the ICAs (if any).
-Did the authors validate the link between methylation and genotype data (e.g. were they from the same person)? -Using genetic variation as causal anchor is a strong approach. However, identifying robust genetic variants requires a GWAS and hence many thousands if not 10 thousands of samples (instead of 500). Were there GWASs performed for cortical thickness the authors can rely on? Gene-set enrichment approaches as applied here reduce the enthusiasm for a genetic approach and limits its interpretation.
-SWAN sometimes induces strange DNA methylation values. It will be useful to exclude such artefacts by comparing the raw DNA methylation values to those post-SWAN.
-How many methylation values were missing and imputed? What were setting of the impute package? -DNA methylation data seems to be generated on a sample taken ~1 years after the MRI. It should be discussed how this affected the study.
-Why were the measured cell counts not included in SVA? I wish to laud the authors for their additional work. Seeking and finding replication (despite the risk of falsifying their original finding) and a considerable effort to gain more insight in the impact of blood cell composition very much strengthens the manuscript. The results are intriguing and accordingly the authors present a cautious discussion on the explanation for and the implications of their findings. So, I am happy with the revisions.
Reviewer #4 (Remarks to the Author): Freytag et al. identified an epigenetic signature that was associated with cortical thickness and memory performance in 533 healthy young adults. The epigenetic signature was replicated in 596 participants with major depressive disorder and healthy controls. Further, the authors showed that this signature mediated the effect of age on cortical thickness. I consider the findings novel and the application of ICA as analytical strategy of interest to the community.
As advised by the Editor, I have specifically focused this review on the aspects of the study that related to the DNA methylation analysis. I fully share the initial concerns raised by Reviewer 3 that: (a) replication of the reported associations in an independent sample cohort is critical; and (b) the reported associations need to be independent of differential blood cell counts.
Re (a): I acknowledge the effort that was necessary to add these data and analyses to the revised manuscript. The authors replicated their primary finding, and therefore have addressed my concerns in this regard.
Re (b): The authors performed the additional analyses requested by the reviewer. The results showed that after adjustment for differential cell counts, the associations of ICA2 with both chronological age and cortical thickness remained significant. Further, the authors showed that the ICA-age correlations identified in whole-blood were also detectable in individual cell types. I note that the authors have added a paragraph to the Discussion section (p. 16), indicating the caveats of DNA methylation profiling in a heterogeneous tissue such as peripheral whole blood. Together, I think the authors have done an adequate job in replicating their findings and assessing confounding due to cellular composition.
Nonetheless, it would be useful if the authors considered the below concerns, suggestions or additions in the revised manuscript: (1) A justification as to why ICA was applied in favour of other more established analytical approaches could be added. In addition, the authors could elaborate on the limitations of ICA in the Discussion section.
(2) It would be important to add a clarification on how significance thresholds were defined, particularly for the association testing between the ICA components and traits of interest. Is there an appropriate equivalent of 'genome-wide significance'?
(3) Related to (2): Some reported correlations have highly significant p-values (p <2.2e-16), yet the corresponding correlation coefficients seem low (r= 0.32). Some associations are dismissed as speculation, despite reporting p <1e-60 as for the case of whole blood vs. monocyte ICA2 DNAm values. I am not challenging the statistics, but the authors could make a better effort in reporting the statistics more consistently and with more caution throughout the manuscript. After all, insight into molecular mechanism is not provided.
Reviewer #5 (Remarks to the Author): Increasing age is tightly linked to decreased thickness of the human neocortex. The authors identified an epigenetic signature that was associated with cortical thickness in 533 healthy young adults and replicated this finding in 596 participants with major depressive disorder and healthy controls. The authors reported that the epigenetic signature mediated the effect of age on cortical thickness and pointed to the involvement of immune system genes.
I believe the findings are robust but have several concerns about the interpretation of the results.
1. Two out of 15 ICA methylomic patterns (termed ICA1 and ICA2) were significantly correlated with age but only ICA2 predicted cortical thickness after controlling for the linear effects of age.
Supplemental table 2 suggests that ICA1 captures the linear affects of age. My question is whether ICA2 captures non-linear components of the effects of age. For this purpose, I propose an analyses where cortical thickness is first predicted from a model containing a 5th degree polynomial of age (i.e., age+age^2+…+age^5). If ICA2 is added to this model, does it significantly increase the overall R^2? I believe, this answer to this question is important because if this test is not significant, the relation between methylation and cortical thickness could be spurious and caused by age affected both cortical thickness and methylation profiles in blood.
2. The authors claim that the epigenetic signature of age on cortical thickness. (p <0.001). This seems too strong given the non-experimental nature of the data and possibility of alternative explanations. A mediator model predicts that the age-cortical thickness correlation equals the product of the ageepigenetic signature correlation times the epigenetic signature-cortical thickness correlation. This model does not seem to hold making the authors to conclude that the epigenetic signature only partly mediated the affects of age. However, I am unsure about this statement, as this is just a decomposition of a correlation and alternative explanations for such data that that do not assume mediation. (e.g., there could be 3rd variable affecting both cortical thickness and methylation). Methods related to this response: We analyzed whole-blood methylomic profiles from n=656 samples published in Hannum et al., 2013. In analogy to our methylomic dataset, multi-mapping or polymorphic probes were excluded from analysis. Raw intensities (methylated and unmethylated signals) were normalized using the lumi package (color-bias adjustment and quantile normalisation). The BMIQ algorithm was finally applied to adjust for the difference between Type I and Type II probes used in the 450K array. Given substantial non-randomness of between-plate distribution of chronological age in this sample, we suggesting that the results presented herein were not driven by global mean thickness.
In addition, F6 scores obtained from the mean thickness-adjusted EFA solution were still significantly associated with EM performance (p = 0.008) and ICA2 (p = 0.042) (see Table 2 in this response).
We added this information on page 9 of the revised manuscript. For each factor from the mean thickness-adjusted EFA (a), the factor from the ICV adjusted EFA (b) exhibiting highest loadings correlation is shown. We are particularly thankful to this reviewer for the opportunity to clarify this point. It made us realize that the manuscript's introduction failed to make completely clear to the reader that we adopted a two-stage study design with regard to the relationship between age, cortical thickness and methylation.
The first analytical step addressed the question whether global methylation patterns are possible mediators of the effect of age on cortical thinning. To this end, ICA was performed first to describe global methylomic patterns. The second stage of the analysis was initiated once such a mediator (ICA2) was identified. Importantly, the mediator exerted a significant effect on cortical thinning also after correction for age. Thus, the second stage consisted of studying the effect of the mediator (ICA2) -after correction for age-on additional relevant phenotypes, such as regional cortical thickness and related cognitive traits. Consequently, the exploratory factor analysis (EFA) was performed after systematically regressing out age effects from MRI measurements (ROIs) and methylomic patterns.
Once again, we thank the reviewer for the opportunity to clarify this and added the improved description of our approach on pages 4 and 6 of the revised manuscript.

Reviewer #2, comment 3:
I understand that ICA2 patterns were used, but were the age-adjusted thickness estimates also used, or rather the raw thickness estimates based on the age-adjusted factors? If the latter is the case, maybe this could be followed up by a mediation analysis where age is included, since it is a major source of variation that is now hidden. As I read the rational for the study, understanding age-related changes in cortical thickness is a major objective, and it would thus be informative if age was not just taken out of the equation. Similar concerns regard the memory analyses.

Authors' response:
This point is closely related to comment 2 of this reviewer (see above) and to the twostage design of our study: The mediator (ICA2) exerted a significant effect on cortical thinning also after correction for age. Thus, the second stage consisted of studying the effect of the mediator -after correction for age-on additional relevant phenotypes, such as regional cortical thickness and related cognitive traits. Importantly, the correlation between memory performance and chronological age in our healthy young population was not significant (r= -0.05, p=0.18). This precluded a mediation analysis between ICA2, memory performance and chronological age.

Reviewer #2, comment 4:
It is mentioned that synaptic pruning during development may account for cortical thinning. It could also be mentioned that myelination of lower cortical layers could cause the cortical mantle to appear thinner on the MR scans. This likely accounts for a substantial part of the observed cortical thinning in development.

Authors' response:
We followed the reviewer's suggestion and added this information on page 14 of the revised manuscript.

Reviewer #2, comment 5:
Could the presentation of which samples that were used for the different analyses be more clearly presented, e.g. in a table, figure or flow chart? If I understand correctly, the age-thickness-methylome associations were found in a sample of you participants only? This should be more clearly stated in the ms, with age-ranges also. One easily gets the impression that this is an older-age study also, and I think this needs to be much more directly addressed, and also touched upon in the discussion, because it will have substantial impact on how the results are conceived.

Authors' response:
Absolutely. As requested by the reviewer, the revised manuscript (suppl. Table 1) contains an improved description of the study populations (including age range). We also address this point on page 13 of the revised manuscript.

Reviewer #2, comment 6:
Minor point: Is it possible to include some for quantitative information in the abstract also, such as sample size and effect sizes? The Fischl et al paper on whole-brain segmentation would not be the most appropriate as the first FreeSurfer-reference here, since the present paper is focused on cortical thickness and the Fischl paper on volume. The 1999 and

2000-papers from Dale and Fischl should rather be referenced first instead.
Authors' response: The suggested papers are now mentioned on pages 9, 21 and 23 of the revised manuscript. We also put quantitative info in the abstract.

Reply to Reviewer #3
We thank this reviewer for the helpful remarks. We are pleased to inform this reviewer that his/her comments motivated us to study additional samples and to provide additional results, which corroborated our conclusions.

Reviewer #3, comment 1:
The straight-forward analysis that should be presented is an epigenome-wide association analysis (EWAS) testing for association between the methylation of individual CpGs and cortical thickness. This approach allows the reader to gauge the relationship between the 'raw' data and the outcomes.

Authors' response:
We are thankful for the opportunity to clarify this point, as it made us realize that the rationale and the two-stage approach of our study was not made entirely clear in the introduction.
Given that the methylome is a high-dimensional space, at least as much as the transcriptome is, our first goal was to search for genome-wide derived methylomic patterns, not singular data points, that might mediate the effect of age on cortical thinning. We were thus interested in gaining a systems-level, genome-wide view of the methylomic signal. Indeed, recent research has demonstrated the importance of analyzing age-associated DNA methylation patterns, as opposed to single CpG sites, when studying the impact of the methylome on physiological processes changing with age, especially when using blood as a surrogate for brain tissue (Horvath et al. Genome Biology 2012, 13:R97). Abundant evidence arising from the study of the transcriptome shows that such patterns can be successfully and robustly identified through the use of such decomposition methods as ICA, resulting in the identification of relevant biological processes (e.g. Biton et al., 2014;Rotival et al., 2011;Wexler et al., 2011; for review on ICA see Kong et al, 2008).
Here, we first addressed the question whether methylation patterns are possible mediators of the effect of age on cortical thinning. To this end, ICA was performed first to describe global methylomic patterns (The reason for choosing ICA over other decomposition methods is described in detail in our response to comment #2 of this reviewer). Using univariate analytical approaches, such as EWAS, where each individual CpG site is independently tested for association with a trait, are ill-suited for detecting methylomic patterns. Said pattern detection can be achieved by projection methods, which decompose the initial high-dimensional dataset into components representing multidimensional DNA methylation patterns amenable to association testing with agerelated traits, and amenable to mediation testing.
The second stage of the analysis was initiated once such a mediator (ICA2) was identified. Importantly, the mediator exerted a significant effect on cortical thinning also after correction for age. Thus, the second stage consisted of studying the effect of the mediator (ICA2) -after correction for age-on additional relevant phenotypes, such as regional cortical thickness and related cognitive traits. Consequently, the exploratory factor analysis (EFA) was performed after systematically regressing out age effects from MRI measurements (ROIs) and methylomic patterns.
Thus, in contrast to decomposed DNA methylation patterns, univariate approaches treating hundreds of thousands of single CpGs as solitary events are ill-suited for addressing the primary question of our study. Once again, we thank the reviewer for the opportunity to clarify this and added the improved description of our approach on pages 4 and 6 of the revised manuscript.
References related to this response:

Authors' response:
We thank the reviewer for this comment and for pointing out that the reason for choosing the specific method (i.e. ICA) was not well-described. As mentioned in the response to comment #1, we performed ICA to achieve a low-dimensional decomposition of independent age-associated DNA methylation patterns and to investigate their relationship to cortical thickness, an age-related complex trait. Both PCA and ICA represent such decomposition methods, yet they rely on different properties of the inferred components. The choice of ICA over PCA was driven by theoretical assumptions regarding the generative model of observed DNAm signals: under this model, the observed DNAm profiles are viewed as a mixture of independent biological processes. By construction, statistical independence of the inferred components implies their non-gaussianity: each component will be characterized by a restricted number of variables, i.e. CpG sites. PCA is based solely on a maximum variance criterion: the initial dataset is projected in a sub-space of successive orthogonal components, each capturing the maximum remaining data variance, irrespectively of any specific generative model underlying the investigated data. In contrast, ICA relies on statistical independency of the components: the observed signal, e.g. DNAm, is assumed to arise from a set of statistically independent sources, which can be viewed as sets of putative independent biological processes influencing the methylome. By assuming statistical independence, ICA decomposition will favor non-gaussian distribution of the inferred components, thus representing independent sources, or processes, acting on circumscribed sets of features (in our case, CpG sites). Thus, ICA decomposition represents a more realistic approach to the analysis of biological data. This approach, assuming non-gaussianity of underlying sources generating the observed signal, has been repeatedly shown to outperform variance-based decomposition approaches for the analysis of gene expression data. (Liebermeister W, 2002;Lee SI & Batzoglou S, 2003;Frigyesi et al, 2006;Teschendorff et al, 2007) and has been successfully applied to multiple DNA methylation datasets (Renard E et al, 2014 (2013)), which consists of n=656 blood DNAm profiles of participants spanning a wide age-range (19-101, mean age: 64 years).
In the first step of the analysis we identified 5 ICA patterns that were significantly associated with age (designated HCIa-e, see Table 1 in this response). In the second step we examined the overlap between ICA1 and ICA2 CpGs identified in our sample of healthy young adults and CpGs contributing to each of the Hannum age-associated IC pattern. A significant overlap with ICA1 was observed for one pattern (OR=90, p<1e-60, see Table 1 in this response). For ICA2 we observed a significant overlap with three age-associated Hannum patterns (see Table 1 in this response), with said overlap being particularly strong for Hannums' HICd pattern (OR=48.7, p < 1e-60). Importantly, the correlation of loadings between CpGs contributing to ICA2 and CpGs contributing to Hannums' HICd pattern was positive and substantial (r=0.87, p < 1e-60).
Thus, we observed highly significant between-sample overlap of ICA patterns despite the differences in age structure of these samples, supporting the robustness of the methods presented herein. We added this data on pages 7 and 8 of the revised manuscript.
The reviewer also raised an issue regarding the number of surrogate variables inferred from SVA (40 variables in our sample), considering this number to be possibly too high and indicative of possible data noisiness. We would like to stress that the number of inferred components is both a function of sample size and dimensionality of the dataset (450000 CpG sites in our case). Thus, it is expected that the number of inferred surrogate variables from the 450K array may largely exceed the number of surrogate variables determined in lower-dimensionality datasets such as gene expression data or the 27K methylation array. Importantly, we followed the same SVA approach in Hannum's dataset (450K array, 656 subjects). As expected, we determined a higher number (n=98) of surrogate variables in this sample. We added this data on page 22 of the revised manuscript.
Finally, a question was raised regarding the number of retained independent components. We applied stringent criteria and discarded from further analysis every component whose variance was explained by more than 10% by one single individual.
Such patterns may or may not represent singular modes of variation, which renders their interpretation difficult at the population level. Importantly, examination of the RLE and Multidimensional Scaling plots did not reveal outlying patterns for the 111 individuals (Figures 1 and 2 in this response) that would call for an a priori exclusion of these subjects from analysis.
In addition, we examined the association of ICA2 with age and age-adjusted cortical thickness considering the 111 individuals as a group covariate: both associations remained highly significant (age: p = 4.91e-12; age-adjusted cortical thickness: p = 4.8e-5). We added this data on page 8 of the revised manuscript. Methods related to this response: We analyzed whole-blood methylomic profiles from n=656 samples published in Hannum et al., 2013. In analogy to our methylomic dataset, multi-mapping or polymorphic probes were excluded from analysis. Raw intensities (methylated and unmethylated signals) were normalized using the lumi package (color-bias adjustment and quantile normalisation). The BMIQ algorithm was finally applied to adjust for the difference between Type I and Type II probes used in the 450K array. Given substantial non-randomness of between-plate distribution of chronological age in this sample, we

MDS5
Since the primary research question relates to the association of whole blood DNA methylation profiles with cortical thickness, SVA should be performed using cortical thickness as outcome, not age. Now, the interpretation of associations of ICA components remains unclear: is it driven by age, an age-effect on cell composition etc?

Authors' response:
This comment is tightly linked to comment #1 of this reviewer and to the fact that the two-stage approach of our study was not made entirely clear in the introduction. The objective of the study was to investigate the relationship between cortical thickness and age-associated DNAm patterns. Considering age as an outcome in the SVA allowed adjusting DNAm for unknown technical/biological confounders, while keeping the influence of age on DNAm intact.

Authors' response:
To address the reviewer's comment, we first applied Horvath's cross-tissue-and Gene-set enrichment analysis for CpGs contributing to ICA2 can be found in supplementary and 1.7% of variance of cortical thickness and spatial F6 factor score, respectively (see Table 1 in this response).
Regarding ICA1, chronological age accounted for 29.7% of variance. This methylomic pattern did not show any significant association with technical covariates nor blood cell counts (measured or estimated) (see supplementary The same analyses were conducted for the remaining 13 ICs (supplementary Tables 12 and 13 of the revised manuscript). However, without a significant association of these ICs with the phenotypes under study (i.e. age and cortical thickness), interpretation of the results is limited.  1175-1185 (2014).

Reviewer #3, comment 6:
The main PCAs in whole blood 450k data usually capture cell counts (beyond the larger classes like the 3 measured here). It is not unlikely that the ICAs also represent cell counts. Hence, if true, the data may suggest an association between the immune system (measured through DNA methylation instead of direct cell counts) and cortical thickness. If this association survives robust control for confounding (i.e. including parameters known to be associated with cortical thickness (and immune system traits) (social economic position, medication use, etc.)). This in fact may be a very interesting observation but requires much more specific insight in the biology of the ICAs (if any).

Authors' response:
We thank the reviewer for this important comment. Examination of the relationship of ICA2 with known technical confounders (e.g. plate, sentrix position), or gross cell count measures did not show any significant association (p >0.01, Supplementary Importantly, after adjusting ICA2 for these blood sub-cell types, the associations with both chronological age and cortical thickness remained highly significant (age: r=0.29, p=2e-11; cortical thickness: r=-0.22, p =8.3e-7).
Gene-set enrichment analysis for CpGs contributing to ICA2 can be found in supplementary These results suggest that association between ICA2 and cortical thickness is not driven by major or subtle changes in blood composition. Nevertheless, ICA2 harbors primarily CpGs known to influence immune function. We comment on this important feature of ICA2 on pages 7, 15 and 16 of the revised manuscript.

Reviewer #3, comment 7:
Did the authors validate the link between methylation and genotype data (e.g. were they from the same person)?

Authors' response:
Absolutely. A per-subject crosscheck between phenotypic data, methylation data and genetic data was performed using the reported sex and sex-predictions based on the array data, as well as matching of all SNPs represented on the Illumina 450K array to the corresponding Affymetrix SNP 6.0 genotype calls. This crosscheck allowed an unambiguous assignment of each methylation dataset to the corresponding genetic and phenotypic dataset. We added this information on page 25 of the revised manuscript.

Reviewer #3, comment 8:
Using genetic variation as causal anchor is a strong approach. However, identifying robust genetic variants requires a GWAS and hence many thousands if not 10 thousands of samples (instead of 500). Were there GWASs performed for cortical thickness the authors can rely on? Gene-set enrichment approaches as applied here reduce the enthusiasm for a genetic approach and limits its interpretation.

Authors' response:
We would like to stress that we did not perform a GWAS on cortical thickness, nor was this the aim of the study. After having identified a methylomic pattern (i.e. ICA2) mediating the effect of age on cortical thinning, we searched for possible genetic patterns that might explain in part this epigenetic phenotype's variability. To this end we capitalized on the power of gene-set enrichment analysis (GSEA) to detect genetic pathways that are associated with the trait of interest. As shown repeatedly by us (e.g.

Reviewer #3, comment 10:
How many methylation values were missing and imputed? What were setting of the impute package?

Authors' response:
The average per CpG and per sample missing rates were < 1e-04; CpGs with missing rate > 0.05 were excluded from analysis. The impute package was used by taking default settings (k=10 nearest neighbors).

Reviewer #3, comment 11:
DNA methylation data seems to be generated on a sample taken ~1 years after the MRI. It should be discussed how this affected the study.

Authors' response:
DNA methylation profiles were obtained on average 1 year after imaging acquisition. We performed a sensitivity analysis examining the association between ICA2 and cortical thickness after regressing out the difference (delta age) between age at blood sampling and age at MRI assessment from the methylomic pattern. The association remained significant (p = 6.4e-5, r = -0.18 after adjustment for age) indicating that delta age did not affect the results of the study. We added this information on page 19 of the revised manuscript.

Reviewer #3, comment 12:
Why were the measured cell counts not included in SVA?

Authors' response:
SVA is recommended as a robust method for cell-type mixture adjustment (McGregoret al, 2016), allowing correction also beyond gross cell composition measures. Importantly, we did not observe significant association of measured blood cell counts with age nor cortical thickness (suppl.

Reviewer #3 (Remarks to the Author):
I genuinely appreciate the considerable additional work the authors put into the manuscript. I believe this has strengthened the manuscript and improved the clarity of their reasoning. I feel, however, that when adopting non-standard analysis approaches (ICA instead of EWAS) particularly in a novel field like epigenetic epidemiology, the burden of proof that the outcomes are robust is on the authors. I have three major concerns, of which certainly the latter two can be addressed easily.
1. The finding that ICA2 is associated with cortical thickness is based on a single study of 553 individuals. Replication is a key aspect of any genomic association study.

Response to Reviewer # 3-1:
We are particularly thankful to the reviewer for this comment. Motivated by his/her suggestion to address the issue of replication, we reached out to our colleagues at the Max Planck Institute for Psychiatry in Munich. Analyses of the relationship between DNAm and cortical thickness in this sample (termed herein "Munich sample") were done entirely by the Munich team, i.e. we were not involved in order to ensure full independence of the analytical steps.
The Munich sample comprised N=627 participants who underwent MRI assessment and whole-blood methylomic profiling (423 MDD patients, mean age 47.9 ± 13.8; 204 controls, mean age 49.5 ± 13.3). After pre-processing of methylomic data and MRI-QC based exclusion, N=596 subjects remained for association testing of cortical thickness and ICA2 (estimated from whole-blood samples).
The ICA2 pattern was estimated as the linear combination between ICA2 loadings (as inferred from the Swiss DNAm sample) and individual DNAm values of the Munich sample. In this independent sample, we observed a significant positive correlation between ICA2 and chronological age (r=0.48, p < 1e-10; as a reminder: the corresponding value in the Swiss sample was r=0.29) and a negative correlation with global cortical thickness (r= -0.31, p < 1e-10; as a reminder: the corresponding value in the Swiss sample was r= -0.24). After adjustment for chronological age and controlling for potential confounders (diagnosis, sex, intracranial volume, MR-batch effects, time difference between MRI examination and blood drawing), the association between ICA2 and cortical thickness remained significant (r= -0.094, p = 0.011; as a reminder: the corresponding value in the Swiss sample was r= -0.18). Importantly, the same analysis in a sub-sample of N=163 participants younger than 40 years (thus, within an age range similar to that of the Swiss participants) revealed an almost identical effect size (r= -0.19, p = 0.009) compared to that observed in the Swiss sample. Thus, we are pleased to report the independent replication of our primary finding and are therefore particularly thankful to this reviewer for motivating us to put this additional effort into our study.
We added these data on page 8 of the revised manuscript.

Methods related to this response:
Munich sample:  (Stein et al., 2012;Schmaal et al., 2016). Other than in the flagship study (Stein et al., 2012), no bipolar patients were included for reasons of clinical homogeneity (Schmaal et al., 2016).
MDD diagnoses were based on clinical consensus in addition to M-CIDI or SCAN interviews, depending on the original study protocols. After preprocessing of methylomic data, and MRI-QC based exclusions, combined data of N=596 subjects was available for statistical analysis.  (Aryee et al., 2014). Beta values for the pre-selected 397,947 autosomal probes from the Swiss sample were calculated from SWAN normalized intensities. After pre-processing of methylomic data, and MRI-QC based exclusions, combined data of N=596 subjects was available for statistical analysis.
Association with ICA2 pattern: ICA2 patterns were calculated separately for the whole Munich sample and a subsample of <40-year-old subjects (N=163).
DNA methylation values were first adjusted for sex using linear regression.
ICA2 patterns were then calculated as linear combination between the scaled residuals and the inverse ICA2 loadings inferred from the Swiss sample.
Separate Pearson's correlation analyses were performed between ICA2scores and biographical age, and ICA2-scores and CT scores. In addition,

The authors claim that the association of ICA2 with cortical thickness and
other phenotypes is not driven by blood cell counts. In general, cell counts are well-known to be main drivers of major components explaining methylomic variance. Accounting for main cell types is not sufficient to exclude this possibility. The fact that ICA2-CpGs are near genes that are primarily involved in processes related to inflammation and leukocyte differentiation is compatible with a cell count effect. As mentioned in my first review, an immune component related to cortical thickness is of interest but we should know its nature. The authors may use the approach adopted in a recent paper on age-related changes in DNA methylation (Slieker et al. Genome Biol 2016;17:191), where the authors took public 450k data from many blood cell subtypes and then assess whether identified CpGs are differentially methylated between cell subtypes. The authors can do this for their 970 ICA2 CpGs.
3. The authors do use a public data set (Hannum et al) to confirm the presence of ICA1 and 2 and their correlation with age. This data-set was, however, on whole blood. They should extend this analysis to public data on purified blood cells like the one published in Nat Commun covering monocytes and T cells (Reynolds et al. Nat Commun. 2014;5:1-8). This would further substantiate the independence of blood cell counts.

Response to Reviewer # 3-2 & 3-3:
We thank the reviewer for suggesting to provide further evidence for the independence of the results on blood cell counts and for proposing to look deeper into the putative nature of the relation between immunity and cortical thickness by using the data reported in Slieker et al. Here we followed these suggestions and respond to both comments in common because they are closely related.
In addition to the originally reported (Supplementary were observed for regulator and memory CD4+ T-cells ( Figure 3). Generally, the correlation coefficients might suggest high concordance of the cortical thickness-related blood DNAm patterns with DNAm of B lymphocytes and of the common myeloid progenitor lineage, and relatively less concordance with DNAm of natural killer cells and T lymphocytes. Given the correlative nature of all these analyses, we suggest being humble in our interpretation and report this analysis in the main body of the manuscript by highlighting the cautiousness, which should be applied to such conclusions and stating that no further inference can be drawn towards the contribution of a specific immune cell type to the reported associations.
Altogether, we are truly thankful to this reviewer for his/her input, which led to significant improvement of the manuscript. We now provide independent replication of the primary finding and show that the results are independent of blood cell composition-related bias. The new results can be found on pages 7, 11 and 12 of the revised manuscript.

References:
Absher   In each dataset, ICA1 and ICA2 patterns were estimated as the linear combination between the inverse of genome-wide ICA1 and ICA2 loadings (inferred from the Swiss sample) and scaled SV-adjusted DNAm values. This score was subsequently tested for association with chronological age. variance criterion: the initial dataset is projected in a sub-space of successive orthogonal components, each capturing the maximum remaining data variance, irrespectively of any specific generative model underlying the investigated data. In contrast, ICA relies on statistical independency of the components: the observed signal, e.g. DNAm, is assumed to arise from a set of statistically independent sources, which can be viewed as sets of putative independent biological processes influencing the methylome. By assuming statistical independence, ICA decomposition will favor non-gaussian distribution of the inferred components, thus representing independent sources, or processes, acting on circumscribed sets of features (in our case, CpG sites). Thus, ICA decomposition represents a more realistic approach to the analysis of biological data.

Comparison of whole-blood and cell-specific
This approach, assuming non-gaussianity of underlying sources generating the observed signal, has been repeatedly shown to outperform variance-based decomposition approaches for the analysis of gene expression data. (Liebermeister W, 2002;Lee SI & Batzoglou S, 2003;Frigyesi et al, 2006;Teschendorff et al, 2007) and has been successfully applied to multiple DNA methylation datasets (Renard E et al, 2014).
On the limitations side, decomposition of genome-wide methylomic profiles comes at the cost of specificity of the inferred solution towards the genomic localization of CpG markers. The detection of CpGs contributing to the methylomic signature relies on a fixed threshold on the distribution of the components' loadings. In our case, this approach allowed relating ICA2 broadly to genes involved in immune system function.
However, the specific relationships between the identified marker sets and the phenotypes of interest can be studied only in downstream experiments focusing on single CpG sites.
It is also important to stress that the ICA model relies on the assumption that methylomic signals arise from a fixed set of independent sources. In absence of a priori knowledge about the source signal, the number of inferred components must be determined empirically, which might impact negatively on generalizability. Indeed, generalizability is a general challenge of systems-level genome-wide approaches (Ritchie et al., 2015).
Integration of multiple-layers of molecular traits, such as genotypic data used in this study, is therefore important to address whether the identified patterns represent relevant features of the dataset.
subjects, thus resulting in lower p-values, as compared to statistics reported for our N=533 methylomic sample. Further, the p-value < 1e-60 relates to correlation tests between average CpG-centered DNAm values across tissues.
To improve clarity, we now report sample sizes together with statistical values for any test involving different sample sizes than our methylomic sample.

Reviewer #4, comment 4:
Supplementary Figure 1: I assume the y-axis indicates proportion and not percentage?

Reply to Reviewer #5
We highly appreciate this reviewer's constructive remarks and important suggestions, which helped corroborating the conclusions of our statistical analyses.

Authors' response:
We are particularly thankful for this important remark and for the reviewer's proposal to investigate the impact of non-linear effects of age on the reported associations.
As suggested by the reviewer, we performed an F-test analysis to compare the fit of a model predicting cortical thickness from a 5 th degree polynomial of age (age+age^2+…+age^5) to the fit of the same model augmented by ICA2. We observed a highly significant increase in adjusted R 2 with the addition of ICA2 to the model (F(1,507)= 15.6, p =8.8e-05).
Thus, these results indicate that the association between ICA2 and cortical thickness is not driven by non-linear effects of age. These results have been added on page 7 of the revised manuscript.
Methods related to this response: An F-test was performed to compare a full model between cortical thickness as dependent variable and a 5th degree polynomial for age and ICA2 as independent variables to a restricted model without ICA2. reversible upon cessation. Epigenetics. 2014 Oct; 9(10):1382-96.

Reviewer #5, comment 4:
The authors state that DNA methylation age values were significantly correlated with actual participants' age. As DNA methylation age is essentially the deviation from methylation predicted age and chronological age, this should not be the case. It makes me wonder whether this is caused by possible non-linear effects of age remaining in the DNA methylation age values.

Authors' response:
We are thankful for this comment because we realized that the meaning of the term 'DNA methylation age' was not made completely clear to the reader.
By using the term 'DNA methylation age', we refer to the direct outcome of publicly available predictors for chronological age based on methylomic markers. These predictors have been shown to yield estimates for age that are highly correlated with actual participants' age (Hannum et al., 2013;Marioni et al., 2015), as also observed in our methylomic sample. We clarify this now in the manuscript (p. 8).