Introduction

The brain structure aspect alters throughout life. In particular, grey and white matter volumes are known to shrink with age in diseased and normal brains (Ge et al. 2002; Fjell and Walhovd 2010; Lockhart and DeCarli 2014). Numerous cross-sectional and longitudinal studies have confirmed this trend by either studying grey matter volume changes (Lemaitre et al. 2012) or cortical sulci widening (Kochunov et al. 2005; Shen et al. 2018). Importantly, the magnitude of brain shrinkage varies across regions and individuals, and increases with age (Raz et al. 2005). Multiple factors related to the environment or to genetics likely play a role in these changes. Such genetic effect is characterized in the hippocampus atrophy with the apolipoprotein E, ε4 allele (ApoE-ε4), which is also associated with an increased risk for developing late onset Alzheimer’s disease (AD) (Moffat et al. 2000). However, the genetic underpinnings of the brain sulcal features have not been investigated, except a few studies that were interested in the heritability of sulcal depth in extended pedigrees of young adults (Le Guen et al. 2018) or non-human primates (Kochunov et al. 2010). To the best of our knowledge, no genome wide association studies (GWAS) in imaging genetics with a sample size above 10,000 subjects were conducted on the brain sulcal features. In imaging genetics, previous GWAS with such sample sizes have looked into the hippocampal and intracranial volumes (Stein et al. 2012), or the human subcortical brain structures (Hibar et al. 2015). These studies traditionally used meta-analyses which pooled subjects scanned in different centers with various scanners and from different age ranges. Such a meta-analysis initiative is best exemplified by the ENIGMA (Thompson et al. 2014) and CHARGE (Psaty et al. 2009) consortia. The UK Biobank project (Allen et al. 2012) offers a remarkable opportunity to address these issues by gathering data from a fairly homogenous population of subjects, and acquiring magnetic resonance images (MRI) on identical scanners operated at the same location. Additionally, it enables researchers to directly access a cohort with numerous participants while alleviating the uncertainty of meta-analyses.

In this paper, we consider ten prominent brain sulci that are automatically extracted and labelled using the Brainvisa cortical sulci recognition pipeline (Rivière et al. 2009; Perrot et al. 2011). These sulci are the central, the anterior and posterior cingulate, the inferior and superior temporal, the intraparietal, the subparietal, the superior and inferior frontal sulci and the Sylvian fissure (Mangin et al. 2015; Shen et al. 2018). Even though the subparietal sulcus is not as prominent as the others, it was included in our analysis because it lies in the precuneus which is a major target for atrophy in AD (Karas et al. 2007; Bailly et al. 2015). First, we replicated the known trends of cortical shrinking with age in a large sample of individuals from the UK Biobank, considering the grey matter thickness and sulcal opening. Second, we estimated the genetic influence on these features with the genome-wide complex trait analysis (GCTA) (Yang et al. 2011). In this method the genetic relationship (kinship) matrix between subjects is computed to estimate the variance of an observed phenotype, which can be explained by the single nucleotide polymorphisms (SNPs), referred to as the heritability. Finally, we performed a genome-wide association study (GWAS) of the phenotypes with the genotyped variants using PLINK (Purcell et al. 2007). A functional annotation of the phenotype-associated variants was then performed using the gene expression level published by the Genotype-Tissue Expression (GTEx) consortium (GTEx Consortium 2017), allowing the identification of expression quantitative trait loci (eQTLs).

Methods

Materials—subjects

The present analyses were conducted under UK Biobank data application number 25251. The UK Biobank is a health research resource that aims to improve the prevention, diagnosis and treatment of a wide range of illnesses. Between the years 2006 and 2010, about 500,000 people aged between 45 and 73 years old, were recruited in the general population across Great Britain. In this work, we used the data released on January 2018, consisting of 20,060 subjects with a T1-weighted MRI. We included 15,040 in our discovery sample and 5,020 in our replication sample. The subjects were separated according to their availability in NIFTI format. In January 2018, replication sample was available only in DICOM format. The imaging quality control (QC) was performed by UK Biobank following information described until and including the NIFTI conversion stage in (Alfaro-Almagro et al. 2018).

The UK Biobank genetic data underwent a stringent QC protocol, which was performed at the Wellcome Trust Centre for Human Genetics (Bycroft et al. 2017). We restrained our analysis to people identified by UK Biobank as belonging to the main white British ancestry subset (using the variable in.white.British.ancestry.subset in the file ukb_sqc_v2.txt). Additionally, we excluded from our analysis subjects with high missingness, high heterozygosity, first degree related individuals or sex mismatches. In total 12,162 subjects in the discovery cohort and 3,435 subjects in the replication cohort passed the image processing steps and the genetic criteria filtering (see study design on Figure S1). Both sets include approximately 48% of males and 52% of females.

Cortical sulci extraction

The cortical sulci were extracted from T1-weighted images via the following steps. First, the brain mask was obtained with an automated skull stripping procedure based on the SPM8 skull-cleanup tool (Ashburner 2009). Second, the images were segmented into grey matter, white matter and cerebrospinal fluid using histogram scale-space analysis and mathematical morphology (Mangin et al. 2004). Third, individual sulci were segmented and labelled using Morphologist, the sulci identification pipeline from Brainvisa [version 4.5.1, (Rivière et al. 2009)]. For segmentation, a kind of crevasse detector was used to reconstruct each fold geometry as the medial surface from the two opposing gyral banks that spanned from the most internal point of the fold to the convex hull of the cortex (Mangin et al. 2004). A Bayesian inspired pattern recognition approach relying on Statistical Probabilistic Anatomy Maps and multiscale spatial normalization was used to label the folds using a nomenclature of 125 sulci (Perrot et al. 2011; Mangin et al. 2015). For each sulcus, the average distance between both banks of the pial surface was used to quantify the sulcus width. This average distance was computed as the ratio between (1) the volume of cerebrospinal fluid (CSF) filling up the sulcus from the brain hull to the fold bottom and (2) the surface area of the sulcus estimated by half the sum of the areas of the triangles making up a mesh of the corresponding medial surface (Fig. 1a). The average thickness of the cortical mantle on both sides of the sulcus was computed using a fast marching algorithm applied to a voxel-based binary representation of the cortex grey matter (Perrot et al. 2011).

Fig. 1
figure 1

Evolution of sulcal opening and grey matter thickness with age in the main brain sulci. a Schematic definition of the opening and grey matter thickness for the Brainvisa sulci. b Age distribution in the UK Biobank sample with MR Imaging. c Mean opening (red) and mean grey matter thickness (blue) vs age. (light blue and light red represent left hemisphere values; dark blue and red represent right hemisphere values). Main sulci annotation is described on Fig. 2

Fig. 2
figure 2

Heritability estimates of the sulcal opening (a) and sulcal grey matter thickness (b). (p < 0.00125 Bonferroni corrected). Main sulci have been annotated using Brainvisa abbreviation (SC Central, FIP intraparietal, FCLp Sylvian fissure, FCMpost posterior cingulate, FCMant anterior cingulate, STs superior temporal, STiant inferior temporal, SFsup superior frontal, SFinf inferior frontal, Ssp subparietal). The sulci are displayed using the statistical probability anatomy map (SPAM) representation, which represents the average sulci shape and position on the reference base of the Brainvisa sulci extraction pipeline (Perrot et al. 2011)

Age and sulci features relationship

In our discovery sample, to quantify the influence of age, we computed the mean of each sulcal feature for each age across the subjects with the same age. It is worth noting that the age at MRI scan is provided by UK Biobank with a 1 year precision. We excluded 45 and 46 year old subjects due to their small sample sizes. We can robustly estimate the mean feature per year after 50 years, because there are more than 200 subjects per 1-year interval. Last, using linear regression, we estimated the slope of a linear model adjusted between sulcal features and age (Table S1).

Heritability estimation—genome-wide complex trait analyses (GCTA)

In our discovery sample, we used GCTA (Yang et al. 2011) that yields an estimate of the heritability (h2SNPs) in population studies with genotyped unrelated participants. We considered the genotyped SNPs variants common to the UKBiobank and UKBileve arrays (details at http://www.ukbiobank.ac.uk). To compute the kinship matrix of the population, specific SNPs were selected with PLINK v1.9 (Purcell et al. 2007) using the following thresholds: missing genotype = 5% (70.783 variants excluded), Hardy–Weinberg equilibrium test (hwe) = 10−6 (11.318), and minor allele frequency (maf) = 1.0% (102.559). We kept the SNPs in moderate linkage disequilibrium with variation inflation factor 10 within a window of 50 SNPs (92.081 variants excluded). Then, we computed the genetic relationship matrix with GCTA using the 507.515 SNPs left. The amount of phenotypic variance captured by this matrix is estimated using a mixed-effects linear model. As covariates in our genetic analyses, we systematically included the gender, the genotyping array type, the age at the MRI session and the 10 genetic principal components (PCs) provided by UK Biobank to account for population stratification (see Figure S2 for these PCs). Correction for multiple comparisons were achieved using Bonferroni correction accounting for all our phenotypes and we retained as significant p < 0.00125 = 0.05/40 (2 hemispheres × 2 cortical features × 10 sulci). Using the GCTA Power Calculator (Visscher et al. 2014), the discovery cohort sample size provides above 99% statistical power to detect heritability values above 15%, with p < 0.00215.

Genetic univariate association analyses

The genotype-phenotype association analyses were performed using PLINKv1.9 (Purcell et al. 2007), with the following thresholds: missing genotype = 10% (32.938 variants excluded), hwe = 10−6 (12.301), and maf = 1.0% (117.165), in total 621.852 variants passed the genetics QC for the discovery sample. Note that the missing genotype threshold is traditionally more stringent to estimate the kinship in GCTA analysis than to perform GWAS. As covariates in our genome-wide association analyses, we included the gender, the genotyping array type, the age at the MRI session and the 10 genetic principal components provided by UK Biobank to account for population stratification.

To estimate the multiple testing burden due to our analyses of 40 phenotypes, we used matrix spectral decomposition (matSpD, (Nyholt 2004)) which yielded 19 independent traits (see correlation matrix Figure S3). Therefore, we applied the significance threshold p < 2.6·10−9 for our discovery cohort. In addition, we also reported association that passed the suggestive genomic threshold p < 5·10−8 in our discovery sample and p < 0.05 in the replication sample for variants that passed the first threshold.

Functional annotation of loci

The functional annotations of loci were provided by FUMA (Watanabe et al. 2017) which notably annotates SNP location, expression quantitative trait (eQTL) effects [e.g., GTEx, the UK Brain Expression Consortium (http://www.braineac.org/)], and chromatin state.

Results

Age-related cortical shrinking

The UK Biobank large sample size enabled us to precisely estimate the mean of sulcal opening and grey matter (GM) thickness per age with a 1-year precision. Figure 1c underlines a strong correlation between the age and these two cortical features in the discovery sample. Between 47 and 73 years old, the sulcal opening increases on average of 0.025 mm/year, while the GM thickness decreases on average of 0.015 mm/year (Table S1). The Sylvian fissure and the subparietal sulcus show the maximum increase of opening and decrease of GM thickness, respectively, while the inferior temporal and left intraparietal have the minimum.

Heritability of sulcal opening and grey matter thickness around the sulci

In the discovery sample, we estimated the heritability (h2SNPs) of the sulcal opening and GM thickness using the GCTA method (Yang et al. 2011), i.e. the variance explained by all the SNPs (Fig. 2). Significant h2SNPs estimates for the sulcal opening range from 0.15 to 0.45, with a minimum in the left inferior temporal and a maximum in the left central sulcus (Table S2). Significant h2SNPs estimates for the GM thickness range from 0.15 to 0.37, with a minimum in the left superior frontal and a maximum in the left Sylvian fissure. Note that the sulcal opening heritability values are all higher than the ones of the GM thickness. The significant heritability estimates obtained support the fact that these two phenotypes are largely influenced by independent genetic factors and thus are good candidates for the following genome-wide association study.

Genome-wide association study (GWAS) of the cortical features

We performed a genome-wide association study on the genotyped data for the sulcal opening and GM thickness in the ten sulci. Manhattan and QQ plots are shown in the supplementary materials (Figures S4–11). Table 1 summarizes the 24 phenotype-SNP associations that were genome-wide significant in the discovery sample and nominally significant in the replication sample. Among these associations, 12 SNPs were unique at 5 different loci. The most represented locus, with 17 replicated association hits, is on chromosome 1, 27 kb before the transcription start site of the KCNK2. Within this locus, the two main associated SNPs are rs59084003 and rs864736 with 7 and 4 significant hits respectively. It should be noted that these two SNPs are not in strong linkage disequilibrium (LD) (r2 = 0.06) and other significant SNPs on this locus are either in LD with the first or the second one (Figure S12). The second most represented locus, with 4 significant associations, is on chromosome 16. This locus lies in the starting region of the non-coding RNA LOC101928708, in the vicinity of the protein coding gene C16orf95, followed by FBXO31. On this locus, the main associated SNP is rs9933149 with 3 hits. The three other loci each include a single significant replicated phenotype-SNP association. On chromosome 8, rs11774568, associated with the GM thickness in the left Sylvian fissure, is in a region with a high density of genes, between the genes DEFB136 and DEFB135. On chromosome 9, rs10980645, associated with the central sulcus opening, is an intronic variant of the LPAR1 gene. On chromosome 12, rs12146713, associated with the right STs opening, is an intronic variant of the NUAK1 gene.

Table 1 Significant genome-wide association hit SNPs (discovery p < 5·10−8 and replication p < 0.05)
Fig. 3
figure 3

GWAS hits upstream of KCNK2 regulating the sulcal opening. First and second lines correspond to rs864736 and rs59084003, respectively. Lines represents respectively: a the log10 (p value) of each SNPs mapped onto the nominally significant sulci among the ten considered; b the mean sulcal opening and standard error for each configuration of variants in the most significant sulci; c Locuszoom display (Pruim et al. 2011) of the phenotype-variants association for the region upstream of KCNK2 with the left posterior cingulate sulcus opening as a phenotype

Regional significance and direction of effect of the hit variants

DNA region upstream KCNK2 harbors significant SNP association with the sulcal opening (Fig. 3) and GM thickness (Figure S12) of the posterior cingulate, the central and intraparietal sulci. In addition, the implicated SNPs show near genome-wide significant influence on the sulcal features in the superior temporal, inferior frontal and subparietal sulci. Overall, this suggests a brain wide regulation for this genomic region. The sulcal opening is increased in carriers of the minor allele of rs864736 (maf = 46%, in our discovery sample), while it is decreased in carriers of the minor allele of rs59084003 (maf = 7%). The locus on chromosome 16 displays a more specific spatial control over the temporal and frontal lobes, with significant sulci including the inferior frontal, the superior and inferior temporal and the Sylvian fissure (Figures S13–14 a.). The sucal opening of these sulci is decreased in carriers of the minor allele rs9933149 (maf = 38%). It should be noted that the SNPs of these two loci have a significant pleiotropic influence on sulcal opening and GM thickness.

Regarding the three remaining loci, they preferentially influence one out of the two sulcal features. The intronic variant (rs12146713) of NUAK1 significantly influences the sulcal opening in the temporal region (Sylvian fissure, superior and inferior temporal sulci), as well as in the inferior frontal and intraparietal sulci (Figures S13-14 b.). The intronic variant (rs10980645) of LPAR1 significantly affects the sulcal opening of the central sulcus, and moderately affects the posterior cingulate (Figures S13-14 c.). The variant (rs11774568) chromosome 8, near DEFB136, appears to be linked with the GM thickness of the Sylvian fissure and superior temporal sulcus (Figures S7-8 d.).

Functional annotation of the loci

To obtain information on the role of the region upstream of KCNK2, we investigated the gene expression QTL using the GTEx database [GTEx Analysis Release V7 (dbGaP Accession phs000424.v7.p2)] (GTEx Consortium 2017). We found that the variants rs864736 and rs14526008, which are in LD, are significant eQTLs for the KCNK2 gene expression, with significant multi-tissue meta-analysis RE2 (random-effects model 2) (Han and Eskin 2011) p values of 9.7·10−6 and 5.7·10−9, respectively. Figure S15 presents the effect size in various tissues of rs864736 allelic configuration on KCNK2 gene expression. Even though the association barely reaches nominal significance in single brain tissue due to low sample sizes (~ 80–140 subjects), it overall emphasizes strong effects across brain tissues. The other significantly associated SNP near KCNK2 (rs59084003) was not reported as a significant eQTL in GTEx possibly for technical reasons because of its small minor allele frequency (7%). Indeed, the effect of allelic configurations cannot be well observed in the relative small GTEx sample. Table S3 summarizes the other eQTLs found in GTEx for the different loci.

In addition, rs864736 chromatin-state is annotated only for Brain Cingulate Gyrus (E069) tissue as being a recruitment site for PolyComb repressive complex (ReprPCWk, grey track, Fig. 4), which is an epigenetic system of gene silencing. This annotation in this specific tissue is particularly interesting because the strongest phenotypic associations with this variant are the ones with the posterior cingulate sulcus opening. We note also that the chromatin-state around rs504473 (in strong LD with rs864736) is annotated as enhancer (Enh, yellow track) specifically in brain tissues (Fig. 4) and thus this region seems to be crucial in the fine regulation of transcription of nearby genes, notably KCNK2, in brain tissues.

Fig. 4
figure 4

Functional annotation (eQTLs and chromatin-state in brain tissues) of the genomic region in the vicinity of variant rs864736 [with FUMA (Watanabe et al. 2017)]. Brain tissues names: E054: Ganglion Eminence derived primary cultured neurospheres, E053: Cortex derived primary cultured neurospheres, E071: Brain Hippocampus Middle, E074: Brain Substantia Nigra, E068: Brain Anterior Caudate, E069: Brain Cingulate Gyrus, E072: Brain Inferior Temporal Lobe, E067: Brain Angular Gyrus, E073: Brain Dorsolateral Prefrontal Cortex, E070: Brain Germinal Matrix, E082: Fetal Brain Female, E081: Fetal Brain Male, E125: Astrocytes Primary Cells. RoadMap Core 15 chromatin state model: TssA active transcription start site (TSS), TssFlnk flanking active TSS, TxFlnk transcription at gene 5′ and 3′, Tx strong transcription, TxWk weak transcription, EnhG genic enhancers, Enh enhancers, ZNF/Rpts ZNF genes & repeats, Het heterochromatin, TssBiv bivalent/poised TSS, BivFlnk flanking bivalent TSS/enhancer, EnhBiv bivalent enhancer, ReprPC repressed polycomb, ReprPCWk weak repressed polycomb, Quies quiescent/low

Discussion

The bulk of participants in the UK Biobank are older than 50 years old and results have to be discussed within this interval. To scrutinize certain genetic effects on brain features, it is important to understand the effect of age. Indeed, some genetic influences may only be revealed in a relatively aged cohort such as the UK Biobank. In this study, we emphasized a steady age effect on the sulcal opening and GM thickness between 45 and 75 years of age. Furthermore, we underlined that these brain features are heritable, thus their variance across individuals is significantly regulated by gene additive effects. Finally, we highlighted several causal genetic variants whose allelic configuration directly contributes to the phenotypic variance.

It is well known that with age, the GM and white matter (WM) volumes decreases (Raz et al. 2005), while the amount of cerebrospinal fluid (CSF) in the cortical folds increases (Good et al. 2001). In our analysis, we confirmed the well-established results regarding the GM thickness decrease and sulcal opening increase with age in the UK Biobank cohort. The assessment of cortical sulci opening is well described and correlates with neurocognitive decline in mild cognitive impairment and dementia disorders (Bastos Leite et al. 2004). The cortical sulci widening with age is likely related to the reduction of gyral thickness resulting in the dilatation of the sulci (Magnotta 1999; Symonds et al. 1999), but could also account for neurodegenerative processes occurring in the underlying white matter (Gunning-Dixon et al. 2009). While, the GM thickness and sulcal widening are correlated (Kochunov et al. 2008; Liu et al. 2013), the robustness of the two measures might differ. We obtained less significant variant-phenotype associations in GM thickness, because the sulcal opening might be more consistently measured across individuals than the GM thickness. This might be due to the fact that the MRI contrast between GM and CSF remains more stable across the lifespan than the GM / WM contrast (Kochunov et al. 2005). Thus, sulcal widening is commonly used by radiologists as a surrogate of cortical atrophy in clinical settings (Shen et al. 2018). It could also reflect a higher sensitivity to the primal effects of aging because of the consequences of grey and white matter atrophy. Note that the opening of larger sulci like the Sylvian valley could also be impacted by global mechanical compensation for aging processes.

The main finding of this paper is that the locus upstream of the KCNK2 transcription start site influences the sulcal opening and GM thickness. Additionally, the tissue specific gene expression (eQTL) analysis of the GTEx consortium emphasizes that overall (meta-analysis of all tissues), and particularly in the brain, this DNA region regulates the expression of KCNK2 (GTEx Consortium 2017). Thus, we can legitimately assume a link between the regulation of KCNK2 expression and the amplitude of sulcal opening. In other words, depending on the allelic configurations in the region upstream of KCNK2, an individual will have his sulci comparatively enlarged. Because sulcal widening is a marker of cortical atrophy as we pinpointed in the previous paragraph, there is a potential link between KCNK2 expression level and brain atrophy. The KCNK2 gene, also known as TREK1, is a member of the two-pore-domain potassium channel family which is expressed predominantly in the brain (Hervieu et al. 2001). Previous literature emphasized several functions for KCNK2 gene in the brain. First, the KCNK2 regulates the blood–brain barrier function and inflammation in the brain of mice (Bittner et al. 2013) and humans (Bittner et al. 2014). The inhibition or deletion of KCNK2 facilitates lymphocytes migration into the central nervous system (CNS) and promotes autoimmune CNS inflammation (Bittner et al. 2013). Second, in mice, the knockdown of KCNK2 gene impairs the neuronal migration of late-born cortical excitatory neurons, which are precursors of Layer II/III neurons (Bando et al. 2014). Third, in rat hippocampal astrocytes, the increase of KCNK2 expression mediates neuroprotection during ischemia (Banerjee et al. 2016). The mechanism might involve KCNK2 blockade, inhibiting neuronal apoptosis and protecting the brain from cerebral ischemic injury (Wang et al. 2018). Finally, KCNK2 over expression was shown to exacerbate memory impairment in middle-age mice (Cai et al. 2017). To summarize, KCNK2 controls several major cellular responses involved in memory formation and is believed to participate in neuroinflammation, cerebral ischemia and blood–brain barrier dysfunction (Bittner et al. 2014; Cai et al. 2017; Wang et al. 2018).

The first role suggests the most promising direction of future work, because previous studies have proposed that neuroinflammation is involved in cognitive decline in midlife (Marsland et al. 2015) and implicated in pathological age-related changes and AD (McGeer and McGeer 2001). Throughout life, stress, recurrent inflammation and subclinical cerebrovascular events potentially contribute to brain aging (Raz and Rodrigue 2006). The link between our findings and inflammation indicates a potential mediation role for KCNK2. Finally, it is difficult to disentangle whether or not brain inflammation has a deleterious role on cognitive functions, since there is no clear consensus. A recent study however emphasizes a slower progression of AD in patients with early neuroinflammation (Hamelin et al. 2016).

One limitation of this study is that the genetic associations were found in a British sample of elderly individuals. Thus, this association might not replicate in a sample composed of individuals with different ancestry or a younger cohort for which the brain features under scrutiny have not sufficiently evolved with age or with potentially recurrent neuroinflammatory events.

In conclusion, in a sample of 15,597 subjects representative of the general population of British ancestry, we have shown that an eQTL of KCNK2 influences sulcal widening. This appears coherent with the role of KCNK2, which affects the regulation of inflammation response in the brain.