### Introduction

### Methods

### Study participants

### Genotype data

^{−6}, resulting in 4,518,929 SNPs. For sample QC, we only preserved participants with consistent sex and calling rates >90%, resulting in 9,331 individuals. The SNP QC and sample QC were performed using PLINK (v1.90) [17]. We then imputed the genotyped SNP data using the Beagle 5.0 software [18] after the SNP and sample QCs. We restricted the analysis to Hapmap3 SNPs with MAF >0.05 and imputation R

^{2}>0.80, which consisted of 3,848,960 SNPs.

### GWAS analysis

^{−8}) that were independent from each other and formed clumps of other SNPs that were located within 1,000 kb from each index SNP and in LD with the index SNP based on r

^{2}> 0.1. Clumps annotated to the same gene were further combined to the clump with the most significant index SNP.

### Meta-analysis using trajectory model

*i*at time

*j*:

*y*

*is a phenotypic value;*

_{ij}*u*

_{0}

*is a random effect allowing variation around the intercept*

_{i}*β*

_{0};;

*u*

_{1}

*is a random effect allowing variation around the slope for age covariate,*

_{i}*β*

_{1};

*x*

*(*

_{hij}*h*= 1, …,

*q*) are

*q*covariates including age (

*h*= 1) with fixed effects

*β*

*; and*

_{h}*e*

*is a random error. We assumed that the random intercept (*

_{ij}*u*

_{0}

*) and the random slope (*

_{i}*u*

_{1}

*) were independent of each other and normally distributed. The models were adjusted for age squared, sex, drinking, and smoking status. We then used the predicted random slope (*

_{i}*u*

_{1}

*) from the model as the phenotype for a GWAS analysis of trajectory outcomes. All GWAS results for average and trajectory outcomes were obtained using BOLT-LMM software. In the GWAS for trajectory outcomes, the models were adjusted for only 10 genotype PCs to avoid over-adjustment.*

_{i}### Meta-analysis using a random-effects model

### Grid-based Bayesian mixed model

*y**= (*

_{i}*y*

_{i}_{1}, …,

*y*

_{in}_{i})

*is an*

^{T}*n*

*× 1 phenotype vector of individual*

_{i}*i*;

*μ**=*

_{i}*μ*

**1**

_{n}_{i}is an

*n*

*× 1 overall mean vector;*

_{i}*q*time/environmental covariates,

*p*SNPs, two-way interactions among

*p*SNPs (resulting in total of

*p(p–1)/2*terms), and

*pq*SNP-time/SNP-environment interactions;

**= (**

*β*

*β*

^{tT}*,*

**β**

^{gT}*,*

**β**

^{ggT}*,*

**β***)*

^{gtT}*is a vector of genetic effects, time/environmental effects, epistasis effects, and SNP-time/environment interactions;*

^{T}

**e***is an*

_{i}*n*

*× 1 vector of random errors with*

_{i}

**e***~*

_{i}*N*(0,

*σ*

^{2}

*I*

_{n}_{i}). For a given trait, we had

*n*individuals, where individual

*i*has phenotypic values measured at

*n*

*time points (*

_{i}*i*= 1, …,

*n*), and the total number of observations was

*k*pre-specified grid points,

*= (*

**t***t*

_{1}, …,

*t*

*)*

_{k}*, and defined*

^{T}

**ν***as a*

_{i}*k*× 1 vector of random effects at the grid time points with

**ν***~*

_{i}*N*(0,

**), where**

*D***is a**

*D**k*×

*k*covariance matrix. We also defined the incidence matrix

*p**as*

_{i}*k*observations measured exactly on the

*k*grid time points, then

*p**became an identity matrix. We applied a linear interpolation procedure to any observation that did not fall on any one of*

_{i}*k*grid time points. When the

*j*th measurement of individual

*i*fell at time

*t*, which was in between the grid points

*t*

*and*

_{r}*t*

_{r}_{+1}(

*t*

*≤*

_{r}*t*≤

*t*

_{r}_{+1}), we set

*t*=

*t*

*, we get*

_{r}*k*×

*k*covariance matrix

**, resulting in the decomposition,**

*D***=**

*D*

*ΔΨΨ*

^{T}**, where**

*Δ***is a nonnegative**

*Δ**k*×

*k*diagonal matrix and

**is a**

*Ψ**k*×

*k*lower triangular with 1’s in the diagonal elements. We re-parameterized model (1) as

*b**= (*

_{i}*b*

_{i}_{1}, …

*b*

*)*

_{ik}*such that*

^{T}*b*

*~*

_{ij}*N*(0,1) and

*j*= 1, …,

*k*.

**and**

*Δ***, we defined two vectors**

*Ψ***= (**

*δ**δ*

*:*

_{l}*l*= 1, …,

*k*)

*and*

^{T}**= (**

*ψ**ψ*

*:*

_{ml}*m*= 2, …,

*k*;

*l*= 1, …,

*m*− 1)

*. The prior distribution for*

^{T}**is**

*δ***was**

*ψ*

*ψ***and**

_{0}

*R***are pre-specified hyperparameters. The joint prior distribution of**

_{0}*a*th genetic effect was a normal distribution,

*χ*

^{2}distribution. The prior for the overall mean

*μ*was given by

*σ*

^{2}was chosen as a scaled inverse

*χ*

^{2}distribution,

_{2}(BF) > 6.8). A critical issue with the proposed Bayesian model was how to choose an optimal number of grid points,

*k*. We achieved the goal by evaluating the goodness of the predictive distributions of our Bayesian models. Spiegelhalter et al. [11] proposed the DIC and Ando [12] developed the BPIC. We select the optimal number of grid points for our model by minimizing the DIC or simplified BPIC scores. The proposed methods were implemented in an R package named Gridbayes [10], which is available for download at https://github.com/wonilchung/GridBayes.

### Results

### GWAS analysis

*R*

^{2}> 0.8 were retained for the GWAS analyses, resulting in 3,848,960 SNPs.

^{−8}), and for HIP, five SNPs on chromosome 11 (baseline), one SNP on chromosome 12 (third follow-up), two SNPs on chromosome 12 (fourth follow-up), and one SNP on chromosome 6 (sixth follow-up) exceeded the threshold for GWAS significance in Supplementary Table 1. For WST, 11 SNPs on chromosome 11 (baseline) and for WHR, 16 SNPs on chromosome 12 (second follow-up), one SNP on chromosome 12 (fifth follow-up), one SNP on chromosome 18 (second follow-up), and one SNP on chromosome 18 (fifth follow-up) reached GWAS significance (see Supplementary Table 1).

### Meta-analysis

_{10}(p-value) in GWAS for four longitudinal traits between the Meta-Random method and Meta-Fixed method. As expected, the p-values with the Meta-Random method were smaller than those obtained with the Meta-Fixed method, but there was reasonable concordance between them. In Fig. 3, it can be seen that p-values of the Meta-Random method were more significant than those of the baseline analysis in general.

### Bayesian analysis

_{2}(BF) > 6.8) for BMI, 20 SNPs for HIP, 10 SNPs for WST, and 6 SNPs for WHR in Supplementary Table 4. We compared the results of our Bayesian analysis with p-values from the baseline, average, trajectory, and Meta-Random methods and discovered various newly detected novel loci and reasonable concordance among different approaches.

*y*

*is a phenotypic value,*

_{ij}*u*

_{0}

*is a random effect allowing variation around the intercept*

_{i}*β*

_{0},

*β*

*denotes genetic and nongenetic effects for eight covariate values*

_{h}*x*

*(including SNP-age interaction), and*

_{hij}*e*

*is a random error. Based on p-values for SNP-age interaction in the LMM, we found 30 significant SNPs interacting with age for BMI, three SNPs for HIP, three SNPs for WST, and one SNP for WHR in Table 1, meaning that almost 25% of the identified SNPs using Gridbayes had significant p-values for SNP-age interaction. Based on the results of SNP-age interaction, as shown in Table 1, we displayed age-specific changes in BMI for three different genotypes (0, 1, and 2) at four SNPs (Chr5:149489242, Chr6:96377858, Chr10:18571215, and Chr20:40534573), age-specific changes in HIP at one SNP (Chr8:98187086) and age-specific changes in WST at one SNP (Chr2:133806432) in Fig. 5. These figures clearly showed different slopes in phenotypic values over six time points for three genotypes, meaning that there were SNP-age interactions at the identified genetic loci by Gridbayes.*

_{ij}