Dimensional gender diversity is associated with greater polygenic propensity for cognitive performance and interacts with other genetic factors in predicting health outcomes

2.1 Primary study: SPARK Research Match

SPARK [16] is a U.S.-based nationwide autism study of over 270,000 participants, in which genetic data is available for many of the participants. For the primary study we conducted a Research Match, in which independent adults, with or without autism, were invited to participate. Those who agreed and consented to participate were asked to complete the Gender Self Report (GSR), the Adult Self Report (ASR) [?], and additional questions regarding their sexual orientation, gender identity, and gender expression, with the final sample size N = 818. This study was approved by the University of Iowa Institutional Review Board (IRB #201611784). SPARK is approved by the Western IRB (#20151664).

2.2 Secondary analysis: Larger SPARK sample

The secondary analysis used the SPARK adult background history and child background history data from the Version 7 phenotype data release. The child background history (parent-report) was filtered to individuals age 12 or older and with cognitive age level not significantly below age. This was merged with the adult background history (self-report). SPARK is approved by the Western IRB (#20151664).

2.3 Genotype quality control and imputation

Version 3 Freeze (2019) and Version 4 (2020) genotypes were first merged using PLINK [20]. The merged genotypes were then lifted from hg38 to hg19 using the LiftOver tool [21]. The genotypes included 43,209 individuals and 616,321 variants that were then quality controlled using the BIGwas quality control pipeline [22], which performed the genotype quality control, sample quality control, and identification of population stratification and sample filtering due to genetic ancestry. The BIGwas default parameters were used, except for skipping Hardy-Weinberg tests and keeping related samples due to the SPARK cohort being family-based and not a general population sample (we later removed related individuals with GCTA [23] within our SPARK sample of interest). The pre-QC annotation step removed 21 variants (N = 616,299 variants remaining). The SNP QC step removed 101,600 variants due to missingness at a threshold of 0.02 (N = 514,699 variants remaining). The sample QC step removed 1,114 individuals due to missingness, 67 individuals due to heterozygosity, and 176 due to duplicates (or monozygotic twins). The population stratification step projects the remaining individuals onto the principal components (PCs) from the combined HapMap and 1000 Genomes PCs and removed individuals who are not within median +/- five times the interquartile range for PC1 and PC2. This removed an additional 9,533 individuals (N = 32,422 individuals remaining). The quality controlled set of N = 514,699 variants and N = 32,422 individuals were then imputed to the TopMed [24] reference panel using the Michigan Imputation Server [25] with the phasing and quality control steps included and to output variants with imputation quality r2 > 0.3. After the genotype imputation, the variants were filtered to only the HapMap SNPs (N = 1,054,330 variants) with imputation quality r2 ≥ 0.8 using bcftools [26]. Next, they were lifted over from hg38 to hg19 using the VCF-liftover tool (https://github.com/hmgu-itg/VCF-liftover) and the alleles normalized to the hg19 reference genome. Finally, the files were converted to PLINK files with N = 1,018,200 final variants.

2.4 Polygenic score calculations

Polygenic scores were calculated using LDpred2 [27] and the bigsnpr tools [28] in R [29]. Because SPARK is family-based, an external LD reference based on 362,320 European individuals of the UK Biobank (provided by the developers of LDpred2) was used to calculate the genetic correlation matrix, estimate heritability, and calculate the infinitesimal beta weights. Polygenic scores were calculated from the following genome-wide association studies performed by the Psychiatric Genomics Consortium: ADHD (2019) [30], anorexia nervosa (2019) [31], autism (2019) [32], bipolar disorder (2021) [33], major depression (2019) [34], OCD (2018) [35], and schizophrenia (2020) [36]. Polygenic scores were calculated from genome-wide association studies performed by the Social Science Genetic Association Consortium for cognitive performance (2018) and educational attainment (2018) [37] and from the UK Biobank for same-sex sexual behavior [13]. The public LDpred2 beta weights from the Polygenic Index Repository [38] were used to calculate polygenic scores for depressive symptoms [39], extraversion [40], loneliness [41], neuroticism [42], openness [43], risky behavior [44], subjective well-being [45], age at first birth [46], number of children ever born (men) [41], and number of children ever born (women) [41].

From the N = 818 Research Match participants, N = 813 also had mental health outcomes data, and N = 457 also had genetic data that passed quality control. This subset was pruned to remove related individuals using GCTA [23] with a relatedness threshold of 0.05 (N = 7 individuals removed, N = 450 remaining). Genetic principal components were re-calculated for these N = 450 participants using the snp autoSVD function from the bigsnpr package [28]. The polygenic scores were residualized using linear regression for the first five genetic principal components, autism diagnostic status, sex, and age and then centered to have a mean of 0 and scaled to have a standard deviation of 1. For the N = 3,466 participants from the larger SPARK sample with sex/gender concordance data available, the same polygenic score processing was performed.

2.5 Polygenic score analyses

In the N = 450 Research Match participants, polygenic scores were correlated with the two GSR scores (Binary and Nonbinary) and the two mental health outcomes from the ASR (Externalizing and Internalizing) using Spearman correlations. For the genetic liability analysis, we residualized the effects of all 20 polygenic scores from both the GSR scores and ASR scores by linear regression. In the correlations stratified by polygenic risk, we grouped people into one of three groups for each polygenic score: upper 75th quartile (N = 113), middle quartile (N = 224), and the lower 25th quartile (N = 113). We then compared Spearman correleation coefficients between the upper versus lower quartiles. We tested for GSR score by polygenic group interaction effects in association with the ASR scores with linear models: ASR score ∼ GSR score + polygenic group + GSR score:polygenic group.

For the N = 3,466 participants from the larger SPARK sample with sex/gender concordance, t-tests were performed.

3.1 Distribution of the GSR scores

The demographic characteristics of the SPARK cohort are shown in Table 2. Approximately one-quarter of the cohort is transgender or gender nonbinary (TGNB). The distribution of the two GSR scores are shown in Figure 2. The two GSR measures are Binary gender diversity, the extent to which one identifies with the binary gender opposite that of their sex at birth, and Nonbinary gender diversity, the extent one identifies with a neutral gender and/or both binary genders. The scores range from 0 (no gender diversity) to 1 (high gender diversity), with the mode being for both near 0 (Figure 2A). The overall trends show higher gender diversity in females (Figure 2B) and autistic individuals (Figure 2C). TGNB and non-heterosexual participants trended towards greater Binary and Nonbinary (Figure 2D and E).

• Download figure
• Open in new tab

Figure 2. Distribution of the two Gender Self Report (GSR) scores.

A Distribution of the GSR scores in N = 450, stratified by B designated sex at birth, C autism diagnostic status, D gender identity: cisgender or TGNB (transgender and gender nonbinary), and E, sexual orientation.

3.2 Polygenic score correlations with GSR scores

The two GSR scores were correlated with twenty human behavior polygenic scores using Spearman correlations, with the correlation coefficients ρ shown in Figure 3A. As expected, higher polygenic scores for same-sex sexual behavior were significantly positively correlated with both Binary (ρ = 0.11, p < 0.05) and Nonbinary (ρ = 0.11, p < 0.05). Interestingly, higher polygenic scores for cognitive performance were positively correlated with Binary (ρ = 0.09, p = 0.07) and Nonbinary (ρ = 0.13, p < 0.05), meaning that polygenic propensity for greater cognitive performance is associated with elevated Binary and Nonbinary gender diversity. To a lesser extent, educational attainment polygenic scores were positively correlated Nonbinary (ρ = 0.08, p = 0.09), but only moderately correlated with Binary (ρ = 0.05, p = 0.26). No neuropsychiatric polygenic scores were significantly correlated with either GSR score. The strongest correlations were Binary being positively correlated with schizophrenia polygenic scores (ρ = 0.07, p = 0.13) and OCD (ρ = 0.06, p = 0.22), and Nonbinary being negatively correlated with bipolar disorder (ρ = -0.06, p = 0.22).

• Download figure
• Open in new tab

Figure 3. Polygenic score correlations with the GSR scores.

A Polygenic score correlations with the GSR scores. The dotted lines represents the minimal correlation coefficient for statistical significance (nominal p < 0.05), which is ρ = 0.09. Polygenic score abbreviations: cog perf = cognitive performance, edu attain = educational attainment, dep = depressive symptoms, extra = extraversion, lonely = loneliness, neurot = neuroticism, open = openness, risky = risky behavior, SWB = subjective well-being, AFB = age at first birth, NEB = number of children ever born, SSSB = same-sex sexual behavior, ADHD = attention-deficit hyperactivity disorder, MDD = major depressive disorder, OCD = obsessive compulsive disorder, SCZ = schizophrenia.

3.3 ASR score correlations with GSR scores

We next wanted to disentangle the relationship between mental health outcomes, GSR scores, and polygenic scores in the N = 450 participants. The two mental health outcomes were the two major subscales from the Adult Self Report (ASR) [15], Internalizing problems and Externalizing problems. Internalizing problems is an overall composite score of anxious depressed, withdrawn depressed, and somatic complaints. Externalizing problems is an overall composite score of aggressive behavior, rule-breaking behavior, and intrusive behavior. We first correlated the ASR scores with the GSR scores, in which all four were significantly positively correlated (Figure 4A). Externalizing problems was positively correlated with Binary (ρ = 0.12, p < 0.05) and Nonbinary (ρ = 0.18, p < 0.05). Internalizing problems was also positively correlated with Binary (ρ = 0.13, p = 0.06) and Nonbinary (ρ = 0.20, p < 0.05).

• Download figure
• Open in new tab

Figure 4. Correlations and interactions between GSR scores, ASR scores, and polygenic scores

A Correlations of the GSR scores (Binary and Nonbinary) with the two ASR scores (Externalizing and Internalizing). B Correlations of the ASR scores with the polygenic scores. The dotted lines represents the minimal correlation coefficient for statistical significance (nominal p < 0.05), which is ρ = 0.09. C Correlations of the GSR scores with the ASR scores stratified by polygenic scores: the upper 75th quartile and the lower 25th quartile groups (N = 113 for each group). The dotted lines represent the minimal correlation coefficient for statistical significance for N = 113, which is ρ = 0.18. E Four examples of the GSR-ASR score correlations stratified by polygenic group. The β in text is the interaction β.

3.4 ASR score correlations with polygenic scores

We next correlated the ASR scores with polygenic scores (Figure 4B). Externalizing was not significantly correlated with any polygenic scores, but it approached significance for loneliness (ρ = 0.08, p = 0.09) and ADHD (ρ = 0.07, p = 0.15). Internalizing was positively correlated with neuroticism (ρ = 0.16, p < 0.05), major depressive disorder (ρ = 0.10, p < 0.05), and depressive symptoms (ρ = 0.10, p < 0.05).

3.5 ASR score correlations with GSR scores stratified by polygenic risk

We first tested if controlling for genetic liability by residualizing out the effects of all polygenic scores from the ASR scores would attenuate the positive correlations. The correlations were hardly changed, with the strongest attentuation being Externalizing problems and Nonbinary, from ρ = 0.18 to ρ = 0.17.

Lastly, we wanted to test if there were differences in the positive relationships between mental health and gender diversity between high versus low polygenic risk groups. We grouped participants into one of three groups for each polygenic score: high risk (upper 75th quartile, N = 113), middle quartile (N = 224), and the low risk (lower 25th quartile, N = 113), and then performed ASR-GSR correlations between the upper 75th quartile and the lower 25th quartile and also tested for interaction effects with linear models: ASR score ∼ GSR score + polygenic group + GSR score:polygenic group. Figure 4C shows the ρ for GSR-ASR correlations for the upper quartile versus lower quartile polygenic risk groups, and Table 5 has the ρ difference between the two quartiles as well as the p-value for the GSR-polygenic group interaction term. Figure 4D shows four examples of the strongest interaction effects.

We observed dramatic differences in correlations when stratified by polygenic risk, most prominently with neuropsychiatric polygenic scores. The upper quartiles for anorexia, bipolar, and schizophrenia polygenic groups had the strongest correlations, whereas the lower quartiles had correlations of roughly zero. The Nonbinary-Internalizing correlation in the schizophrenia upper group was ρ = 0.32, p < 0.05 and in the lower group was ρ = 0, with the interaction being β = 0.32, p < 0.05 (Figure 4D). A similar effect was observed for the Nonbinary-Externalizing correlation–– the schizophrenia upper group was ρ = 0.27, p < 0.05, and the lower group was ρ = 0.01, with the interaction being β = 0.21, p = 0.12. The Nonbinary-Internalizing correlation in the bipolar upper group was ρ = 0.34, p < 0.05, but in the lower group was ρ = 0.07, p = 0.27, with the interaction being β = 0.31, p < 0.05 (Figure 4D). The Binary-Internalizing correlation in the anorexia upper group was ρ = 0.21, p < 0.05 and in the lower group was ρ = 0.01, with the interaction being β = 0.22, p = 0.09.

We also observed interaction effects with cognitive, personality, and sexuality polygenic groups. Age at first birth (AFB) is positively genetically correlated with educational attainment [11]). The Binary-Internalizing correlation in the AFB upper group was ρ = 0.03, but the lower group (i.e. highest polygenic quartile for younger age at first birth) was ρ = 0.25, p < 0.05, with the interaction being β = -0.29, p < 0.05 (Figure 4D). The Binary-Externalizing correlation in the extraversion upper group was ρ = 0.04, but in the lower group (i.e. highest polygenic quartile for introversion) was ρ = 0.22, p < 0.05, with the interaction being β = -0.24, p < 0.06 (Figure 4D).

4.1 Limitations

Most common genetic analyses (genome-wide association studies, heritabilities, polygenic scores) require moderate sample sizes due to the small effects of individuals common variants. Our major limitation is the small sample size, and we therefore were only powered to detect strong polygenic score correlations. Future replication work in larger samples is crucial. Additionally, age, sex, at birth, and autism diagnostic status can be confounding with the results. Autism is confounding at the genetic level due to previous work showing autism and educational attainment to be positively genetically correlated [32] and higher autism polygenic scores positively correlated with cognitive performance in a general population sample [50]. Because of our small sample size, we could only residualize out the effects of age, sex at birth, and autism diagnostic status from both the GSR scores and the polygenic scores. However, future work should analyze the interaction between autism, sex, and polygenic scores in their associations with gender diversity by performing stratified analyses on autism and sex.