Analyzing longitudinal trait trajectories using GWAS identifies genetic variants for kidney function decline

Abstract

Understanding the genetics of kidney function decline, or trait change in general, is hampered by scarce longitudinal data for GWAS (longGWAS) and uncertainty about how to analyze such data. We use longitudinal UK Biobank data for creatinine-based estimated glomerular filtration rate from 348,275 individuals to search for genetic variants associated with eGFR-decline. This search was performed both among 595 variants previously associated with eGFR in cross-sectional GWAS and genome-wide. We use seven statistical approaches to analyze the UK Biobank data and simulated data, finding that a linear mixed model is a powerful approach with unbiased effect estimates which is viable for longGWAS. The linear mixed model identifies 13 independent genetic variants associated with eGFR-decline, including 6 novel variants, and links them to age-dependent eGFR-genetics. We demonstrate that age-dependent and age-independent eGFR-genetics exhibit a differential pattern regarding clinical progression traits and kidney-specific gene expression regulation. Overall, our results provide insights into kidney aging and linear mixed model-based longGWAS generally.

Fig. 1. Conceptual illustration of genetic variant association with eGFR over time/age and phenotypic models.

a Genetic variant (SNP) associations with eGFR can arise through one allele (risk allele A) that accelerates eGFR-decline over time/age (left) or lowers eGFR in a constant fashion over time/age (right) as compared to the other allele (a). This suggests that genetic variants associated with eGFR-decline are found among genetic variants associated with eGFR cross-sectionally. Shown is a schematic for persons with A/a versus a/a. b Temporal change of eGFR can be modeled in longitudinal data in various ways (phenotypic models): as (i) difference between last and 1st eGFR value of a person (difference model; assessments in-between 1st and last unused and thus depicted as circles); (ii) eGFR over time via linear mixed model (LMM) with person-specific intercepts and slopes (LMM time model RI&RS; time = 0 corresponds to an individual’s 1st eGFR assessment); (iii) eGFR over age (LMM age model RI&RS); or (iv) eGFR over age without random slopes (LMM age model RI-only; time model RI-only possible, but not applied/shown). Shown is a schematic of the phenotypic modeling for two example persons.

Fig. 2. Twelve variants identified for eGFR-decline by focused search among 595 variants.

We selected 595 SNPs previously reported for association with eGFR in cross-sectional data and tested them for association with eGFR-decline using the one-stage LMM age model RI&RS 350K (UKB 350K; n = 348,275, m = 1,520,382). a Shown are P values (P_decline) versus chromosomal position. We identified 12 variants (10 loci) for eGFR-decline at Bonferroni(595)-corrected significance (P_decline < 0.05/595 = 8.4 × 10⁻⁵, brown dashed horizontal line; including 6 with P_decline < 5 × 10⁻⁸, red dashed horizontal line), consisting of 5 novel and 7 known variants for eGFR-decline (blue or green, respectively). Also color-coded are two variants known for eGFR-decline not identified here (orange) and three variants known for not being associated with eGFR-decline (red). Variants with small minor allele frequency (MAF < 5%) are shown as circles. b Shown are genetic effect sizes for eGFR-decline (β_decline from LMM age model RI&RS 350K) versus effect sizes for association with eGFR cross-sectionally (β_{cross-sectional}: eGFR~sex, age, SNP, PCs; eGFR from UKB baseline study-center assessment, n = 341,073). Color and symbol codes are as in (a), additionally highlighting 11 stable-effect variants (black; P_main < 5 × 10⁻⁸, |β_main| > 0.50 mL/min/1.73 m²/allele; P_decline ≥ 0.1; |β_decline| < 0.005 and SE_decline < 0.005 mL/min/1.73 m²/allele and year) that include the CPS1 variant (rs1047891; red in (a)). Effect allele was the cross-sectionally eGFR-lowering allele (unconditioned analyses in EUR). The exact numerical values are provided in Supplementary Data 2.

Fig. 3. Differential pattern between decline-associated versus stable-effect loci regarding age-dependency, clinical progression traits, and tissue-specific gene expression regulation.

We contrasted the 12 decline-associated variants versus 11 stable-effect variants and underlying loci. a Shown are genetic effects on eGFR for 40-, 50-, 60-, 70-year-old individuals using LMM age model RI&RS 350K (beta derived as β_main + (age-50)*β_decline) for decline-associated variants (left; blue: novel, green: known) and stable-effect variants (right; black). Effect allele was the cross-sectionally eGFR-lowering allele (Supplementary Data 3). b We tested the 12 + 11 variants for association with two clinical progression traits using UKB 150K, rapid decline (n_cases = 1211, n_controls = 63,392, logistic regression) and decline in CKD (n_CKD =13,116, m_CKD = 116,944, LMM time model RI&RS; “Methods” section and Supplementary Table 6). Significant enrichment (P_enrich < 0.05) of directionally consistent nominally significant associations was found among the 12 (left; 8/12, 4/12), but not among the 11 SNPs (right; 0/11, 1/11). c We evaluated genes in loci of the 12 + 11 variants regarding tissue-specific enrichment of differentially expressed genes (DEGs): shown are enrichment P values in decline-associated loci (left, among 256 genes) and stable-associated loci (right, among 182 genes; using FUMA, testing 54 tissue types, showing top 25; “Methods” section). Significant enrichment for DEGs (FDR < 0.05, red) was found for decline-associated loci only in kidney cortex (upregulated) and for stable-effect loci in various tissues (mostly downregulated, e.g., in liver, heart, muscle, pancreas, kidney cortex).

Fig. 4. LongGWAS is viable with GMMAT/MAGEE and identifies five loci with genome-wide significance for eGFR-decline.

We conducted a genome-wide search for genetic variant association with eGFR-decline (P_decline, GC-corrected, lambda = 1.06) using the LMM age model RI&RS 350K implemented in GMMAT/MAGEE^, (UKB 350K; n = 348,275, m = 1,520,382; testing 11 million SNPs with MAF ≥ 0.5%, imputation quality INFO ≥ 0.6). a Shown are association P values versus chromosomal position. We identified five loci at genome-wide significance (P_decline < 5 × 10⁻⁸; red dashed horizontal line). Coloring highlights the overall 11 loci identified for eGFR-decline: 10 loci around the 12 variants identified by 595-search (P_decline < 0.05/595 = 8.4 × 10⁻⁵, brown dashed horizontal line; 4 novel and 6 known for eGFR-decline in blue or green, respectively), and one novel locus for eGFR-decline now identified by longGWAS (cyan; lead variant rs2075570 in the 424 loci, but not among the 595 variants). Loci were derived by clumping based on variant position (d > 500kB between loci, “Methods” section). b Shown is the Quantile–Quantile (QQ) plot comparing the distribution of observed P_decline with the distribution of P_decline expected under the null hypothesis of “no association with eGFR-decline” (green: all variants; cyan: excluding the 10 loci around the 12 decline-associated variants; black: excluding the 424 loci around the 595 variants).