Genetic contributions to brain serotonin transporter levels in healthy adults

The serotonin transporter (5-HTT) critically shapes serotonin neurotransmission by regulating extracellular brain serotonin levels; it remains unclear to what extent 5-HTT levels in the human brain are genetically determined. Here we applied [11C]DASB positron emission tomography to image brain 5-HTT levels and evaluated associations with five common serotonin-related genetic variants that might indirectly regulate 5-HTT levels (BDNF rs6265, SLC6A4 5-HTTLPR, HTR1A rs6295, HTR2A rs7333412, and MAOA rs1137070) in 140 healthy volunteers. In addition, we explored whether these variants could predict in vivo 5-HTT levels using a five-fold cross-validation random forest framework. MAOA rs1137070 T-carriers showed significantly higher brain 5-HTT levels compared to C-homozygotes (2–11% across caudate, putamen, midbrain, thalamus, hippocampus, amygdala and neocortex). We did not observe significant associations for the HTR1A rs6295 and HTR2A rs7333412 genotypes. Our previously observed lower subcortical 5-HTT availability for rs6265 met-carriers remained in the presence of these additional variants. Despite this significant association, our prediction models showed that genotype moderately improved prediction of 5-HTT in caudate, but effects were not statistically significant after adjustment for multiple comparisons. Our observations provide additional evidence that serotonin-related genetic variants modulate adult human brain serotonin neurotransmission.


Participants
Cross-sectional data were included that was collected previously and were available from the Cimbi database 38 .
We selected healthy participants based on the following inclusion criteria: (1) availability of BDNF val66met (rs6265) and SLC6A4 (5-HTTLPR and SNP rs23351) genotypes, (2) availability of blood samples for additional genotyping, (3) availability of a [ 11 C]DASB PET scan before any intervention, 4) ⩽60 years of age (to avoid age-related biases in brain volumes, as partial volume effects become stronger after 60), 5) self-identification with European ancestry.In addition, we excluded participants who had the following: (1) diagnosis of a severe neurological or systemic disease; (2) diagnosis of a primary psychiatric disease; (3) substance or drug abuse, based on self-reported clinical history and neurological/physical examination.
We identified 140 healthy participants, 84 females and 56 males (mean age: 26.7 ± 7.2; range: 18-51).The sex imbalance is partly because some studies from which data are drawn included only males or females.Demographic data are described in Table 1.
Subsets of the PET and genetic data included in the current study were collected as parts of multiple previous studies and have been included in previous publications.PET scans included in the current analyses were acquired between 2005 and 2015 11,18,[42][43][44][45][46] .All research protocols were approved by the Ethics Committee of Copenhagen and Frederiksberg, Denmark ((KF) 01-124/04, (KF) 01-156/04, (KF) 01 2006-20, H-1-2010-085, H-1-2010-91, H-2-2010-108, amendments included).All participants provided written informed consent after receiving a detailed description of the respective study.All experimental procedures were carried out in compliance with the declaration of Helsinki.
Data included in the current study has been utilized in previous studies 9,10,23,45,47 , some of which focused on the relation between 5-HTTLPR and/or BDNF rs6265 [ 11 C]DASB binding 23,45 .

MRI data acquisition
For each participant, high-resolution T1-weighted structural brain scans were acquired on either a Siemens Magnetom Trio 3 T (N = 81) or a Siemens Verio 3 T (N = 59) (Siemens, Erlangen, Germany magnetic resonance imaging (MRI) scanner.We used structural MRI scans for segmentation and to delineate regions of interest in the PET scans.

C]DASB PET data acquisition
We acquired PET scans for each participant on one of two PET scanners: 1) a Siemens ECAT high-resolution research tomography (HRRT) scanner operating in 3D-acquisition mode with an in-plane resolution of approximately 2 mm (N = 98) or 2) an 18-ring GE-Advance scanner (General Electric, Milwaukee, WI, USA) with a three-dimensional (3D) acquisition mode and an in-plane resolution of approximately 6 mm (N = 42).PET images were acquired on the HRRT scanner using a 6-min transmission scan followed by an intravenous bolus injection of [ 11 C]DASB that was given over 20 s while a dynamic 90-min emission scan was acquired over 36 frames (6 × 10 s, 3 × 20 s, 6 × 30 s, 5 × 60 s, 5 × 120 s, 8 × 300 s, 3 × 600 s).Dynamic PET images were reconstructed with an iterative OP-OSEM3D method using resolution modeling (10 iterations, 16 subsets) 50,51 .Images acquired on the GE-Advance scanner involved a 10-min transmission scan followed by a bolus injection given over 12 s while a 90-min dynamic emission scan was acquired over 36 frames (same time frames as HRRT acquisitions).Dynamic PET images from the GE-Advance scanner were reconstructed using filtered back projection and were corrected for attenuation, dead-time and scatter with a 6-mm Hann filter.
For both scanners, we determined single-subject PET scan motion and realignment using an automatic image registration algorithm 52 .Next, we smoothed the PET images using a 10 mm (HRRT) or a 12 mm (GE-Advance) within-frame Gaussian filter and subsequently aligned the volumes.Using the scaled least squares cost function, we estimated the rigid translation/rotation parameters that align each PET frame to a reference frame with sufficient structural information (frame 26: 20-25 min post injection).We resliced non-filtered PET images and co-registered the high-resolution MR with PET using SPM (HRRT) or automatic image registration (GE-Advance).Co-registration was based on the mean of frames 10-26, i.e., a flow-weighted image; the accuracy of this step was confirmed visually.We automatically delineated brain regions on structural MRI scans using Pvelab 53 .We determined the time-activity curves across gray matter voxels within each region except for the midbrain region, where we used the mean time-activity across all voxels.We determined regional 5-HTT non-displaceable binding potential (BP ND ) using kinetic modeling of regional time-activity curves in PMOD (Zurich, Switzerland) We applied the multilinear reference tissue model (MRTM/MRTM2) with a fixed k2' estimated for each individual in a high binding region (caudate, putamen, and thalamus) and with cerebellum as reference region 54 .We calculated regional BP ND values bilaterally, computing a volume-weighted mean from the right and left hemisphere.

Data analysis
We carried out the statistical analyses in R v4.1.2(https:// cran.r-project.org/).Consistent with previously related analyses, we considered regional 5-HTT BP ND within caudate, putamen, midbrain, thalamus, hippocampus, amygdala and neocortex as our regions of interest 23 .We selected a single neocortical region (including orbitofrontal-, parietal-and occipital cortex and superior frontal-, pre/post central-, superior temporal-and middle/ inferior frontal gyrus) because of previous evidence of high correlation between different cortical areas 55 .For all analyses, we grouped genotypes as follows: SLC6A4 5-HTTLPR and rs23351: L A L A versus S'-carriers (individuals carrying at least one S-or one L G allele); HTR1A rs6295: G-carriers (carrying at least one G-allele) vs CC homozygotes; HTR2A rs7333412: G-carriers vs AA homozygotes; MAOA rs1137070: T-carriers versus CC; BDNF rs6265: met-carriers versus val/val.In addition, we included age, sex, MRI and PET scanner type as covariates, consistent with previous findings 23 .We mean-centered all continuous variables.Although previous evidence showed a seasonal and body mass index (BMI) effect on 5-HTT availability 45,55 , no statistically significant effect of neither daylight minutes nor BMI on 5-HTT BP ND was previously observed on the same cohort 23 so these variables were not included as covariates.We considered additional lifestyle factors, including alcohol consumption (i.e., alcohol units per week), smoking status, total scores of the Pittsburgh Sleep Quality Index 56 and of the Cohen Perceived Stress Scale 57 .None of these measures were significantly associated with regional 5-HTT BP ND (p > 0.05) and were not included as covariates in the main analyses.For all the models, we set the statistical significance threshold to p < 0.05 (two-sided tests).
Association analyses.We fit a linear latent variable model (LVM) to evaluate associations between genotypes and 5-HTT BP ND within our set of pre-defined brain regions (i.e., caudate, amygdala, hippocampus, putamen, thalamus, midbrain and neocortex) as described previously 23 .LVMs are a type of multivariate linear regression that effectively models associations in the presence of inter-correlated variables.In this case, 5-HTT BP ND is highly intercorrelated between the brain regions that we considered 23 .Thus, using the lava package v 1.6.10 58in R, we modeled this shared correlation of regional 5-HTT BP ND values with a latent variable (5-HTT LV ).Next, we modeled all the genotype and covariate effects on 5-HTT LV .We included additional covariance links (caudate-putamen, amygdala-hippocampus and thalamus-midbrain) based on the model framework previously reported 23 .In addition, we used Wald tests of improvement in model fit to find additional paths.To control the www.nature.com/scientificreports/false-discovery rate across all possible paths, we included the paths with a false-discovery rate of p FDR < 0.05 calculated using the Benjamini-Hochberg test across all paths.We used caudate as a reference region, so that the covariate effects reported here can be interpreted as effects on caudate BP ND (corresponding to the reference scale).
In addition, we estimated multiple linear regression models including all genotypes and covariates (age, sex, PET and MR scanner type) for each brain region and we used them to report percent differences in 5-HTT BP ND between genotype groups.We reported all results with parameter estimates and 95% confidence intervals within brackets, including the associated units (Fig. 1).
Finally, we tested for a main effect of each variant using a likelihood ratio test comparing a model including all covariates and genotypes with a model including all covariates but not genotypes.Prediction analyses.To determine whether genotype information predicts regional 5-HTT BP ND , we trained and tested a random forest model 59 .We used the randomForest package v.4.7-1.1 in R with default data sampling and model fitting parameters; this included p/3 sampled features per tree (where p represents the total number of features) and 500 trees per forest.As we were specifically interested in how well genetic information predicted 5-HTT BP ND , above and beyond other covariates, we constrained our feature set to include only the five genotypes and we evaluated the prediction of regional 5-HTT BP ND , adjusted for relevant covariates.First, we fitted a multiple linear regression model regressing each region, e.g., caudate 5-HTT BP ND , against age, sex, PETscanner, and MR-scanner.The residuals from this multiple linear regression model were used as the outcome in our random forest machine learning models with genotype status for rs6265, 5-HTTLPR and rs23351, rs6295, rs7333412, and rs1137070 as model features.Prediction models were estimated using five-fold cross-validation; prediction of 5-HTT BP ND on held out datasets (test data) within each fold was used to determine unbiased, predictive model performance.Model performance was calculated as the root mean-squared error (RMSE) of prediction on test data, across all folds.To account for fold-assignment related variance, model performance was assessed by repeating the five-fold cross-validation 10 times; overall performance was the mean of these 10 resamples.Statistical significance of performance was calculated from an empirical null distribution derived from 10,000 permutations of the resampled residuals.Notably, we performed five-fold cross-validation with 10 resamples within each of the 10,000 permutations so that each permutation more closely reflects our observed model structure.We fitted models for 5-HTT BP ND in each of our seven regions of interest separately, i.e., seven different prediction models.In addition to expressing model performance in terms of RMSE, we express the percent change in RMSE of the prediction models accounting for genotype information (RMSE genotype ) compared to the RMSE computed on residual 5-HTT BP ND values (RMSE residual , i.e., ΔRMSE = 100*(RMSE residual − R MSE genotype )/RMSE residual ).Statistical significance (p-value) of model performance is expressed both uncorrected (p unc ) and corrected (p FWE ) for the seven models estimated, using Bonferroni-Holm, which controls the familywise type-I error rate 60 .

Genotyping
Genotype distribution and allelic frequencies are depicted in Table 1.rs6265, rs1137070, rs7333412, rs6295 and rs1137070 were in Hardy-Weinberg equilibrium (all p > 0.16).MAOA is x-linked so the allele frequency for rs1137070 in men was compared with the frequency in women using a chi-squared test to determine whether there were significant differences between allele frequencies in males and females.We did not find a statistically significant difference in allele frequencies between males and females (P = 0.67).Assessment of Hardy-Weinberg for the 5-HTTLPR is not valid because this polymorphism was an inclusion criterion for some studies from which these data are derived, i.e., participants were not sampled independent of 5-HTTLPR genotype 9,43 .

Association analyses
The likelihood ratio showed that the LVM including genotype information was statistically significantly different from the LVM not including this information (p = 0.002), suggesting that genotype significantly contributed to the model.
The results from the LVM are depicted in Fig. 1.Consistent with previous observations, we observed that all regional 5-HTT BP ND values loaded strongly on to 5-HTT LV (all p < 10 -12 ).After adding genetic variants and covariates (sex, age, PET and MRI scanner) to the model a Wald test did not support the inclusion of additional paths to the model (p FDR > 0.05).

Prediction analyses
Across the seven regions evaluated, we observed that the set of genetic variant features slightly improved prediction of caudate 5-HTT BP ND compared to the model not including genetic information (caudate: RMSE residual = 0.262, RMSE genotype = 0.266, ΔRMSE = 1.6%, p unc = 0.036) (Table 2).However, this effect, and the effect in all other regions was not statistically significant after correction for seven models (Table 2).

Discussion
We observed that MAOA rs1137070 T-carriers had higher 5-HTT availability compared to CC individuals.T-carriers showed higher 5-HTT availability in all the seven brain regions examined, with the highest binding in caudate (~ 11%) and the lowest in amygdala (~ 2%).Conversely, variants in the HTR1A and HTR2A genes were not associated with 5-HTT availability.Despite observing evidence for a statistically significant association, the genetic variants were not significantly informative for predicting brain 5-HTT available above chance.Taken together, our findings support that genetic variation in the MAOA contributes to variation in brain 5-HTT availability in the healthy adult human brain.
MAOA degrades monoamines in the brain, including serotonin, for which it has a preferential affinity compared to its other substrates 61 .MAOA knock-out rodents have increased extracellular serotonin levels and abnormal affective behavior 30,[62][63][64] as well as reduced 5-HTT expression 63,65 , suggesting that genetically altered MAOA signaling can affect regulation of 5-HT levels, which may in turn modulate 5-HTT levels e.g.via downregulation 31 .Similarly, inhibition of MAO activity by monoamine oxidase inhibitors increases extracellular serotonin levels 62 and is associated with reduced 5-HTT BP ND in rhesus monkey and rats 66 .Notably, reductions in 5-HTT BP ND following an acute pharmacologically-induced serotonin increase may not only reflect a downregulation of 5-HTT but also increased serotonin levels competing for the radioligand rather than a change in 5-HTT gene expression 67 .
The rs1137070 T-allele has been previously associated with lower MAOA enzymatic activity compared to the C-allele both in human fibroblasts in vitro and in post-mortem brains 35,68 .Studies on clinical populations have reported mixed findings 41,69,70 .Although some studies suggest a link between the T-allele and increased MAOA mRNA expression in peripheral blood of patients with depression compared to healthy controls 41 , as well as increased vulnerability to depression 41,69 , other studies reported an association between the C-allele and impaired antidepressant treatment outcome in women 70 .
Conversely, we found that healthy human rs1137070 T-carriers, previously associated with low MAOA activity compared to C-carriers, had greater 5-HTT availability.In this case, putatively lower MAOA activity would correspond to greater amounts of intra-and extracellular serotonin.We can speculate that increased 5-HTT availability reflects increased 5-HTT levels, which might be a compensatory mechanism put in place to reuptake the excess serotonin and maintain extracellular serotonin levels constant.Nonetheless, findings from preclinical research point towards an effect opposite to what we observed, whereas the studies in humans provided mixed findings.Thus, the ambiguity provided by previous evidence in humans does not allow us to draw conclusions about the relationship that we detected between MAOA rs1137070 and 5-HTT BP ND .
We did not observe evidence for an effect of HTR2A rs7333412 on 5-HTT availability.A previous study in patients with major depressive disorder, bipolar disorder, and healthy participants reported an effect on thalamus 5-HTT levels 25 .We did not replicate this effect in our study, as indicated by a comparison of Fig. 2 of their manuscript and our observed group differences.This is possibly because our study is based on a larger and more homogeneous cohort of healthy participants, whereas the previous study included patients with major depressive disorder and bipolar disorder as well as individuals with varying ethnic backgrounds.Taken together, our findings do not support that this HTR2A variant is associated with changes in 5-HTT availability in healthy adults.
Similarly, we did not observe evidence that HTR1A rs6295 is associated with 5-HTT availability, suggesting that whatever effects this polymorphism may have directly on the serotonin 1A receptor, those effects do not significantly modulate 5-HTT availability in healthy adults as measured with [ 11 C]DASB PET.
Notably, we previously reported an association between BDNF rs6265 and 5-HTT BP ND , such that met-carriers showed reduced binding in subcortical areas compared to val-homozygotes 23 .Although this effect was marginally above our threshold for statistical significance in the current model, the effect size was very similar, suggesting independent contributions of MAOA rs1137070 and BDNF rs6265 to 5-HTT availability in healthy humans.As we reported previously, 5-HTTLPR was not significantly associated with 5-HTT BP ND in this cohort 23,23 .
Regarding our prediction model, we observed that using genotype information led to a marginal improvement in predicting caudate 5-HTT BP ND vs not using genotype information, but this effect was not significant after correcting for multiple comparisons.In addition, we could not predict 5-HTT BP ND in any other brain regions.This limited performance may be because we evaluated only five variants, whereas genetically induced variation in 5-HTT levels likely stems from many variants.
Previous studies underscore the limited extent to which candidate variants exert main effects on complex behavioral traits or related features of brain activity 20,71 .Direct measures of discrete neurobiological features, e.g., serotonin transporter protein levels, may however be more susceptible to genetic variants that modulate the relevant neurotransmission pathways 72,73 .Nevertheless, alternative genetic analysis strategies such as GWAS would undoubtedly provide a more comprehensive evaluation of genetic contributions to 5-HTT levels in the human brain.However, an exploratory GWAS requires either thousands of datasets or very large effect sizes (i.e., > 20% difference in 5-HTT BP ND , Cohen's d > 1) to establish statistical significance.GWAS-based polygenic risk scores, e.g., for psychiatric disorders or independent-dataset, hypothesis-generation via, e.g., expression quantitative trait loci (eQTL) databases may instead provide informative and statistically viable strategies for resolving genetic contributions to variation in brain serotonin neurotransmission measured with PET.Our cohort of 140 healthy participants stands as the largest single database of 5-HTT PET brain scans in the world.Future studies probing genetic contributions to brain serotonin-related PET scans would likely benefit from pooling data via, e.g., OpenNeuro (https:// openn europ et.github.io/).
Although genetic variation is plausible and partially supported by our findings, environmental factors are also likely to contribute to 5-HTT levels 67 .A previous study based on the same cohort used for the present study reported no association between 5-HTT levels and daylight minutes or body mass index, contrarily to what reported by earlier studies based on smaller cohorts, part of which are included in our analyses 42,55 .In addition, we could not find any effect of variables reflecting lifestyle measures such as smoking, alcohol consumption, sleep and perceived stress on 5-HTT levels, suggesting that such environmental contributions in this cohort did not confound the observed genetic effects.
All participants self-identified with European ancestry.Self-reports of ancestry can be inaccurate 74 and the lack of ethnic diversity in our sample limits the generalizability of our findings.Notably, we exclude that the association between 5-HTT availability and the rs6265 and rs1137070 genotypes is due to a direct effect of genotype on BP ND .BP ND is proportional to the amount of target proteins available for binding (B avail ), i.e. 5-HTT, the affinity constant of the radioligand for its target (K D ) and the free fraction of ligand in the non-displaceable tissue compartment (f ND ) 75 .We infer the observed genetic effects to be primarily related to change in B avail , although we cannot rule out effects on K D or f ND .Nonetheless, this seems unlikely because rs6265 and rs1137070 are proximal to SLC6A4, but a two scan study structure could more directly disentangle effects of B avail and K D 76 .In conclusion, we report evidence for the association between MAOA rs1137070 genotype and brain 5-HTT availability.We did not observe evidence for an effect of HTR1A and HTR2A variants previously associated with brain serotonin markers, suggesting that their contribution may not be relevant to 5-HTT availability in the healthy adult human brain.Future studies considering additional genetic variants as well as environmental factors in larger datasets are critical for improving our understanding of the factors shaping serotonergic neurotransmission in health and disease.

Figure 1 .
Figure 1.(A) Latent variable model (LVM) used to compute the association between each genotype independently from the other genotypes and 5-HTT BP ND in caudate, putamen, midbrain, thalamus, hippocampus, amygdala and neocortex.PET and MRI scanner are not shown but were included as covariates.Yellow hatched boxes to the left represent the genetic variants and the covariates.The genetic variants depicted correspond to the following genes respectively: 5-HTTLPR-> SLC6A4, rs6295-> HTR1A, rs7333412-> HTR2A, rs1137070-> MAOA, rs6265-> BDNF.The orange boxes represent the covariate effects on the latent variable (5-HTT LV ), which is represented by the central blue ellipse.Light blue boxes show the loadings on the latent variable of observed regional 5-HTT BP ND (in the blue solid boxes to the right).β values in light blue and orange boxes indicate the parameter estimates for each model parameter with either its respective p-value (orange boxes) or 95% confidence interval (light blue).Hatched lines between regions indicate interregional shared correlations.Hatched circles on the brain regions represent the included error estimates.Arrows from the yellow to the blue boxes (sex-> caudate, rs6265-> neocortex) represent direct covariate effects on binding.All brain regions significantly loaded on to 5-HTT LV (p < 1 × 10 -12 ).(B) Boxplots showing representative effects of HTR1A rs6295, HTR2A rs7333412 and MAOA rs1137070 on caudate 5-HTT BP ND .The y axis represents 5-HTT BP ND , adjusted for age, sex, MRI scanner, and PET scanner.Gray dots represent 5-HTT BP ND from each participant adjusted for age, sex, MRI and PET scanner.The larger solid dots and lines represent respective group means and ± 1 SD.The boxes represent datapoints from the 25% to the 75% quantile.5-HTT:serotonin transporter; SLC6A4: 5-HTT gene; HTR1A: serotonin receptor 1A gene; HTR2A: serotonin receptor 2 gene; MAOA: monoamine oxidase A gene; BDNF: brain-derived neurotrophic factor gene; MRI: magnetic resonance imaging (MRI); PET: positron emission tomography; BP ND : non-displaceable binding potential.

Figure 2 .
Figure 2. Random forest model performance.The light blue dots represent the individual RMSE values obtained from resampling (for display purposes, the distribution of the RMSE values from resampling is derived from a model run with 100 instead of 10 resamples) of the model including genotype information (RMSE genotype ).The dark blue error plot displays the mean ± standard deviation of the distribution.The red hatched line indicates the 2.5% quantile of the average RMSE value derived from 10,000 permutations in the model that did not include genotype information (RMSE residual ).Dark blue dots below the red hatched line indicate that the model performed significantly better than chance upon adding genotype information.

Table 2 .
Uncorrected (p unc ) and corrected (p FWE ) p-values for each brain region; root mean squared error computed using residual 5-HTT BP ND values (RMSE residual ), using genotype information (RMSE genotype ) and percent change in RMSE between RMSE genotype and RMSE residual (ΔRMSE).* indicates p unc < 0.05.