Large-scale genome-wide analyses of stuttering

Abstract

Developmental stuttering is a highly heritable, common speech condition characterized by prolongations, blocks and repetitions of speech. Although stuttering is highly heritable and enriched within families, the genetic architecture is largely understudied. We reasoned that there are both shared and distinct genetic variants impacting stuttering risk within sex and ancestry groups. To test this idea, we performed eight primary genome-wide association analyses of self-reported stuttering that were stratified by sex and ancestry, as well as secondary meta-analyses of more than one million individuals (99,776 cases and 1,023,243 controls), identifying 57 unique loci. We validated the genetic risk of self-reported stuttering in two independent datasets. We further show genetic similarity of stuttering with autism, depression and impaired musical rhythm across sexes, with follow-up analyses highlighting potentially causal relationships among these traits. Our findings provide well-powered insights into genetic factors underlying stuttering.

Fig. 1. Miami plot of EUR female and EUR male GWAS.

The EUR female association study (top panel) included 570,071 total samples (40,137 self-reported stuttering cases) and 29,449,463 autosomal variants. Nine loci reached genome-wide significance (dotted line, P < 5.00 × 10⁻⁸) through logistic regression (see Methods). The EUR male association study (bottom panel) included 374,279 total samples (38,257 self-reported stuttering cases) and 29,409,446 autosomal variants. Ten loci reached genome-wide significance (dotted line, P < 5.00 × 10⁻⁸) through logistic regression (see Methods). The x axis represents chromosome base pair coordinates in human genome build 37, and the y axis represents observed −log₁₀(P) for each analysis. Annotated genes for each GWAS are the predicted functional gene for each locus (when available) according to the Open Targets Genetics V2G pipeline (see Methods).

Fig. 2. Forest plot showing genetic correlations for stuttering and previously reported comorbid traits.

a, Ancestry-specific and sex-specific genetic correlations performed for each indicated trait with self-reported stuttering in EUR males and EUR females. Each trait is color-coded according to descriptive category (behavioral, circadian rhythm, immune, metabolic, motor, neurological, speech and language). b, Ancestry-specific and sex-specific genetic correlation estimates (EUR male and EUR female) and sex-combined, ancestry-specific (sex-combined EUR) genetic correlation estimates for beat synchronization and autism. Male-specific correlations are designated by triangles, female-specific correlations are designated by circles and sex-combined correlations are designated by squares. Data points represent the correlation coefficients; error bars, SE. Asterisks denote significant genetic correlations by LDSC. See Supplementary Table 15 for full information on traits used for analyses.

Fig. 3. Forest plot showing results of Mendelian randomization analysis for stuttering and previously reported comorbid traits.

a, Debiased inverse-variance weighted Mendelian randomization analysis estimating causal inference (comorbid trait → stuttering). b, Debiased inverse-variance weighted Mendelian randomization analysis estimating causal inference (stuttering → comorbid trait). Data are represented as estimates; error bars, SE. Filled circles are statistically significant (debiased inverse-variance weighted estimator P < 3.33 × 10⁻³ after Bonferroni correction for testing 15 associations that were genetically correlated with stuttering). Full information on traits used for analyses can be found in Supplementary Table 15. Full results along with other Mendelian randomization methods can be found in Supplementary Table 17. ASD, autism spectrum disorder.

Fig. 4. Performance of self-reported stuttering PRS model in independent EUR stuttering datasets.

PRS were developed using EUR male or EUR female GWAS results and applied to clinically validated ISP and self-report Add Health subjects, and demonstrate increased stuttering liability within stuttering cases. The model was developed and trained using the default auto-phi shrinkage parameter through PRScs. LD panels were constructed using 1000 Genomes Project phase 3 EUR reference data. a, AUC model performance in the ISP cohort. b, AUC model performance in the Add Health cohort.

Extended Data Fig. 1. Study design.

a, Overview of primary and secondary GWAS analyses. b, Flow chart of study design and downstream analyses. Primary analyses are indicated in yellow. Secondary analyses are indicated in blue. Green boxes represent analyses that consist of using data from both the primary and secondary analyses. ‘*’ indicate analyses performed in European ancestry only.

Extended Data Fig. 2. Manhattan and Q-Q plots for primary genome-wide association analyses.

Manhattan plots for genome-wide association analyses performed via logistic regression (see Methods). The dotted line represents genome-wide significance (P < 5 × 10⁻⁸). Q-Q plot x-axes represent expected −log₁₀(P-value), and the y-axes represent observed −log₁₀(P-value) with 95% CI indicated by shaded region. a, The AFR female association study included 24,783 total samples (3,008 self-reported stuttering cases) and 29,555,214 autosomal variants. No loci reached genome-wide significance (dotted line P < 5 × 10⁻⁸). λ_gc = 1.0165. b, The AFR male association study included 13,977 total samples (2,480 self-reported stuttering cases) and 29,443,778 autosomal variants. Three loci reached genome-wide significance (dotted line P < 5 × 10⁻⁸). λ_gc = 1.0156. c, The AMR female association study included 74,187 total samples (7,209 self-reported stuttering cases) and 29,814,878 autosomal variants. One locus reached genome-wide significance (dotted line P < 5 × 10⁻⁸). λ_gc = 1.0340. d, The AMR male association study included 45,206 total samples (6,346 self-reported stuttering cases) and 29,759,321 autosomal variants. One locus reached genome-wide significance (dotted line P < 5 × 10⁻⁸). λ_gc = 1.0177. e, The EAS female association study included 13,347 total samples (1,205 self-reported stuttering cases) and 10,595,969 autosomal variants. No loci reached genome-wide significance (dotted line P < 5 × 10⁻⁸). λ_gc = 1.0081. f, The EAS male association study included 7,169 total samples (1,134 self-reported stuttering cases) and 10,586,108 autosomal variants. No loci reached genome-wide significance (dotted line P < 5 × 10⁻⁸). λ_gc = 1.0182. g, Q-Q plot for the EUR female GWAS; λ_gc = 1.0855. h, Q-Q plot for the EUR male GWAS; λ_gc = 1.0834.

Extended Data Fig. 3. Forest plots showing significant differences in effect size heterogeneity.

Wald tests were performed to test for heterogeneity of effect estimates for the 24 sentinel variants that reached genome-wide significance in the eight ancestry- and sex-specific GWAS compared to the other seven primary GWAS. Plot displays the eight sentinel variants that demonstrated significant differences in effect size by ancestry and/or sex. Data points represent effect size +/− standard error. Each genetic ancestry is color coded. Female-specific effect estimates are shown using circles. Male-specific effect estimates are shown using triangles. Effect estimates that are not shown indicate that the variant was not found. The symbol ‘+’ represents the effect of the result in the Discovery GWAS. The symbol ‘*’ denotes significant differences in effect size estimates compared to the Discovery GWAS after Bonferroni correction. Full GWAS samples sizes can be found in Table 1.

Extended Data Fig. 4. Manhattan and Q-Q plots for the sex-combined meta-analyses.

Manhattan plots for sex-combined meta-analyses performed using the inverse variance weighted option in METAL (see Methods). The dotted line represents genome-wide significance (P < 5 × 10⁻⁸). Q-Q plot x-axes represent expected −log₁₀(P-value), and y-axes represent observed −log₁₀(P-value) with 95% CI indicated by shaded region. a, The meta-analysis includes 944,350 total samples (78,394 self-reported stuttering cases) and 29,566,899 autosomal variants observed across EUR male and EUR female GWAS; λ_gc = 1.0586. b, The meta-analysis includes 119,393 total samples (13,555 self-reported stuttering cases) and 29,948,846 autosomal variants observed across AMR male and AMR female GWAS; λ_gc = 1.0070. c, The meta-analysis includes 20,516 total samples (2,339 self-reported stuttering cases) and 10,656,266 autosomal variants observed across EAS male and EAS female GWAS; λ_gc = 1.0023. d, The meta-analysis includes 38,760 total samples (5,488 self-reported stuttering cases) and 29,705,842 autosomal variants observed across AFR male and AFR female GWAS; λ_gc = 1.0019.

Extended Data Fig. 5. Miami plot for female-specific meta-analysis and male-specific meta-analysis.

Ancestry-combined, sex-specific meta-analyses were performed using multi-ancestry meta-regression via MR-MEGA (see Methods). a, The multi-ancestry female meta-analysis (top panel in orange) included 682,388 total samples (51,559 self-reported stuttering cases) and 41,960,804 autosomal variants (see Methods). Five loci reached genome-wide significance (dotted line P < 5 × 10⁻⁸). The multi-ancestry male meta-analysis (bottom panel in blue) included 440,631 total samples (48,217 self-reported stuttering cases) and 41,937,718 autosomal variants (see Methods). Three loci reached genome-wide significance (dotted line P < 5 × 10⁻⁸). The x-axis represents chromosome base pair coordinates in human genome build 37, and the y-axis represents observed −log₁₀(P-value) for each GWAS. Annotated genes for each GWAS includes the predicted functional gene for each loci (when available) according to Open Targets V2G pipeline. b, Q-Q plot for the female-specific meta-analysis. c, Q-Q plot for the male-specific meta-analysis. Q-Q plot x-axes represent expected −log₁₀(P-value), and y-axes represent observed −log₁₀(P-value), with 95% CI indicated by shaded region. Note: P-values were adjusted for genomic control.

Extended Data Fig. 6. Manhattan and Q-Q plot for the ancestry- and sex-combined meta-analysis.

The multi-ancestry meta-regression meta-analysis performed via MR-MEGA (see Methods) included 1,023,243 total samples (99,776 self-reported stuttering cases) and 42,031,853 autosomal variants observed across all eight independent sex and ancestries-specific GWAS. The dotted line represents genome-wide significance (P < 5 × 10⁻⁸). Q-Q plot x-axis represents expected −log₁₀(P-value), and y-axis represents observed −log₁₀(P-value) with 95% CI indicated by shaded region. Note: P-values were adjusted for genomic control.

Extended Data Fig. 7. Partitioned heritability estimates for broad functional categories.

Partitioned heritability was performed using the 52 broad functional annotations defined by Finucane et al. Enrichment is the proportion of h²/proportion of SNPs. The horizontal, dotted line indicates no enrichment. Error bars indicate enrichment estimates +/− standard error. Asterisks indicate enrichment estimates are significant based on a Bonferroni P-value of 9.6 × 10⁻⁴ (52 gene-sets). a, Broad functional categories partitioned heritability estimates in EUR female self-reported stuttering. b, Broad functional categories partitioned heritability estimates in EUR male self-reported stuttering. c, Broad functional categories partitioned heritability estimates in sex-combined EUR self-reported stuttering. Full results are presented in Supplementary Tables 8–10.

Extended Data Fig. 8. Partitioned heritability estimates for brain cell types.

Partitioned heritability was performed using the three cell-type annotations defined by Finucane et al. The −log₁₀(P-value) for enrichment estimates is indicated on the y-axis. The horizontal dotted line indicates significant after Bonferroni correction (P < 1.7 × 10⁻², three tests). a, Brain cell type partitioned heritability estimates in EUR female self-reported stuttering. b, Brain cell type partitioned heritability estimates in EUR male self-reported stuttering. c, Brain cell type partitioned heritability estimates in sex-combined EUR self-reported stuttering. Full results are presented in Supplementary Tables 11–13.

Extended Data Fig. 9. Partitioned heritability estimates for tissue-specific gene expression.

Partitioned heritability was performed using the tissue-specific expression annotations defined by Finucane et al. The −log₁₀(P-value) for enrichment estimates is indicated on the y-axis. The horizontal dotted line indicates significant after Bonferroni correction (P < 6.25 × 10⁻³, eight tests). a, Tissue-specific gene expression partitioned heritability estimates in EUR female self-reported stuttering. b, Tissue-specific gene expression partitioned heritability estimates in EUR male self-reported stuttering. c, Tissue-specific gene expression partitioned heritability estimates in sex-combined EUR self-reported stuttering. Full results are presented in Supplementary Tables 11–13.

Extended Data Fig. 10. Partitioned heritability estimates for tissue-specific chromatin marks.

Partitioned heritability was performed using the tissue-specific chromatin annotations defined by Finucane et al. The −log₁₀(P-value) for enrichment estimates is indicated on the y-axis. The horizontal dotted line indicates significant after Bonferroni correction (P < 2.5 × 10⁻³, 20 tests). a, Tissue-specific chromatin partitioned heritability estimates in EUR female self-reported stuttering. b, Tissue-specific chromatin partitioned heritability estimates in EUR male self-reported stuttering. c, Chromatin partitioned heritability estimates in sex-combined EUR self-reported stuttering. Full results are presented in Supplementary Tables 11–13.