A Multi-Modal MRI Analysis of Cortical Structure in Relation to Gender Dysphoria, Sexual Orientation, and Age in Adolescents

Abstract

Gender dysphoria (GD) is characterized by distress due to an incongruence between experienced gender and sex assigned at birth. Sex-differentiated brain regions are hypothesized to reflect the experienced gender in GD and may play a role in sexual orientation development. Magnetic resonance brain images were acquired from 16 GD adolescents assigned female at birth (AFAB) not receiving hormone therapy, 17 cisgender girls, and 14 cisgender boys (ages 12–17 years) to examine three morphological and microstructural gray matter features in 76 brain regions: surface area (SA), cortical thickness (CT), and T1 relaxation time. Sexual orientation was represented by degree of androphilia-gynephilia and sexual attraction strength. Multivariate analyses found that cisgender boys had larger SA than cisgender girls and GD AFAB. Shorter T1, reflecting denser, macromolecule-rich tissue, correlated with older age and stronger gynephilia in cisgender boys and GD AFAB, and with stronger attractions in cisgender boys. Thus, cortical morphometry (mainly SA) was related to sex assigned at birth, but not experienced gender. Effects of experienced gender were found as similarities in correlation patterns in GD AFAB and cisgender boys in age and sexual orientation (mainly T1), indicating the need to consider developmental trajectories and sexual orientation in brain studies of GD.

1. Introduction

The Diagnostic and Statistical Manual of Mental Disorders, 5th edition (DSM-5) [1] characterizes gender dysphoria (GD) as distress due to an incongruence between a person’s experienced gender and sex assigned at birth. The number of adolescents referred clinically for GD has steadily increased over the past two decades [2,3], with a recent population prevalence estimate of 0.5% in adolescents and adults [4]. Reports indicate a recent reversal in the sex ratio of clinic-referred adolescents, with more adolescents experiencing GD assigned female at birth (GD AFAB) than assigned male at birth (GD AMAB) [2,3]. Despite this increase in GD adolescents, and in particular GD AFAB, the neural basis of GD development remains poorly understood. The current study provides insight by examining three structural brain features—cortical thickness, surface area, and T1 relaxation time—in GD AFAB adolescents not receiving puberty blockers or hormone therapy.

One prior study examined cortical and subcortical gray matter volume in GD individuals ages 10 to 21 years [5] and included GD participants not receiving hormone therapy as well as those receiving puberty blockers or gender-affirming hormone therapy. (Given our focus on GD AFAB not receiving puberty blockers or hormone therapy, we do not report the results of Hoekzema et al. [5] from their sample of GD adolescents receiving puberty blockers or hormone therapy.) GD AFAB (n = 17) and GD AMAB (n = 11) not receiving puberty blockers or hormone therapy generally had cortical and subcortical gray matter volumes similar to individuals of the same sex assigned at birth. In a small number of regions displaying sex differences, gender identity differences were also found: in GD AFAB, a smaller left superior medial frontal cortex volume vs. cisgender (i.e., experienced gender and sex assigned at birth align) girls (n = 52), and in GD AMAB, a smaller volume in the left superior posterior hemisphere of the cerebellum and marginally smaller hypothalamic volumes vs. cisgender boys (n = 44). Hoekzema et al. [5] concluded that their findings supported the idea that some brain regions in GD adolescents develop in the direction of the experienced gender.

The neurohormonal hypothesis is the most commonly invoked theoretical framework for explaining patterns of brain development in GD [6,7,8,9,10,11]. It posits that prenatal androgen exposure organizes the brain, including the organization of neural sex differences that influence sex-differentiated psychological and behavioral characteristics, such as gender identity and sexual orientation. The surge in sex hormones during adolescence influences brain regions to be further expressed in a relatively male- or female-typical manner, based on the earlier prenatal brain organization [12].

Consistent with this framework, adolescent sex differences in cortical structure have been reported. Cisgender adolescent boys have been shown to have larger surface area than cisgender adolescent girls [13,14,15]. Evidence for sex differences in cortical thickness is more mixed, with some studies finding males have thicker cortices [13] and others finding no differences [14,15]. Sex differences in measures sensitive to underlying microstructure, such as T1 relaxation time, have either not been assessed or have not been observed [16,17,18]. However, a study in 8- to 23-year-olds of gray matter density, a value resulting from soft segmentation of T1-weighted (T1w) images and considered to indicate the amount of gray matter in a given voxel, found higher mean values of gray matter density in females than in males in conjunction with smaller regional cortical volumes [19]. One interpretation was that this finding reflected reduced neuronal density and increased neuropil mass in females, consistent with stereologic morphometry findings [20]. Gray matter density calculations rely on T1w contrast, which is predominantly driven by T1 relaxation time (T1), suggesting that maps of T1 (T1 maps) could be sensitive to sex- and gender-related tissue microstructure differences.

In light of the neurohormonal hypothesis and these observations concerning adolescent brain sex differences, adolescents who experience GD might show cortical surface area, thickness, and T1 features consistent with the experienced gender. An important consideration in evaluating this prediction is that puberty blocking and gender-affirming hormone therapies are often prescribed to ameliorate the distress associated with GD and to promote social transition, including development of the desired physical appearance [11,21]. Because hormone therapy may also impact adolescent brain development [12], it is critical to test predictions about brain organization in GD prior to the onset of hormone therapy [6].

The neurohormonal hypothesis has also been invoked to explain the development of sexual orientation [22,23,24,25], but sexual orientation has seldom been investigated in brain studies of GD ([6]; but see [26] for a study in adults). For example, Hoekzema et al. [5] did not examine associations with sexual orientation, likely because all GD AFAB were gynephilic (i.e., attracted to girls/women) and the GD AMAB were mostly androphilic (i.e., attracted to boys/men). As such, sexual orientation could be an important variable to examine in relation to brain features associated with GD [6].

Adolescence is a dynamic time for sexual orientation development. Retrospective reports indicate that sexual attractions first occur around age 10 for cisgender boys and girls [27,28,29], and often before any first sexual behaviors are experienced [30,31]—although some degree of instability in sexual attractions during adolescence exists, suggesting ongoing exploration with sexuality [32,33]. Developmentally, the strength of sexual attractions, independent of the target, is also an important variable to consider. For example, intensity of sexual behaviors increased (i.e., from hugging to having sexual intercourse) with age in adolescents, and the likelihood of engagement in same-sex behaviors increased in later adolescence [34].

The brain undergoes many changes during adolescence. Cortical volume, which can be decomposed into cortical thickness and surface area, is mostly established by approximately age five and is largely reflected by the development of surface area [13]. There is a surge in the growth of surface area in the first year of life followed by more gradual growth until childhood or adolescence, whereas cortical thickness peaks by approximately age four or five and then slowly declines throughout adolescence and into adulthood [13,35,36,37]. Furthermore, surface area is reported to be influenced by different genetic factors than those involved in cortical thickness [38,39] and surface area is affected more by early prenatal environment factors than cortical thickness [40,41]. During adolescence, the commonly observed decreases in cortical thickness are partially attributable to increases in cortical myelination, which continues well into adulthood [42]. T1 can be a sensitive marker of these myelinating processes, including during adolescence [17,18,43,44,45,46].

To investigate associations between cortical structure and GD, the present study examined cortical gray matter in GD AFAB adolescents not receiving hormone therapy relative to cisgender girls and boys. Using the multivariate partial least squares (PLS) approach, we examined three structural brain features (cortical thickness, surface area, and T1) to characterize patterns across 76 brain regions. The primary aim was to identify multivariate brain patterns that distinguished the three groups. A secondary aim was to evaluate group differences in multivariate brain patterns related to sexual orientation and age. Two aspects of sexual orientation were examined: degree of androphilia-gynephilia, and strength of sexual attractions. In line with the prior studies reviewed above, we predicted that there would be sex differences in surface area, cortical thickness, and/or T1 between cisgender boys and girls [13,14,15,19,20]. GD AFAB not receiving puberty blockers or hormone therapy were expected to differ from cisgender girls and to be similar to cisgender boys on these cortical structural metrics [5]. Based on prior findings, it was predicted that this would especially apply to surface area (e.g., [15,16]). Furthermore, the cortical features of gynephilic GD AFAB should especially reflect those of gynephilic cisgender boys [6].

2. Materials and Methods

2.1. Participants

A total of 47 adolescents participated in the study from 2014 to 2018: 16 AFAB adolescents with a DSM-5 diagnosis of GD not receiving hormone therapy, 17 cisgender adolescent girls, and 14 cisgender adolescent boys. Participants were 12 to 17 years old, with a mean age of 15.32 years (SD = 1.64). The majority were “European”/“White” (n = 27, 57.4%) and the remainder were other ethnicities (n = 20, 42.6%) (see Supplementary Materials for more information). Data from four additional participants were removed: one GD AFAB adolescent due to artifacts from a metal orthodontic expander, one cisgender boy due to excessive head motion, and one GD AFAB adolescent and one cisgender boy due to unusable images needed to calculate the T1 maps.

Inclusion criteria for the study were age 12 to 17 years. GD AFAB participants had to have received a diagnosis of GD via clinician assessment. Exclusion criteria were: receiving any form of hormone therapy (apart from oral contraceptive pills), a known disorder of hormone regulation or sex development, any contraindications to magnetic resonance imaging (MRI, e.g., braces, pregnancy), experience of any head trauma, or insufficient English language proficiency to complete study measures. For cisgender individuals, additional exclusion criteria were: a previous mental health diagnosis, being in a special education class in school, involvement with a child protection agency, and feeling uncomfortable with their sex assigned at birth and identifying or wishing to identify as a member of another gender.

GD AFAB participants were recruited from the Gender Identity Service at the Centre for Addiction and Mental Health or were referred by a clinician in private practice specializing in GD. At the Gender Identity Service, all potential participants within the requisite age range were informed of the opportunity to participate in the study. If they were interested in doing so, a member of the study team, who was not involved in their clinical care, shared information about the study, answered questions, performed screening, and obtained informed assent and/or consent to participate. Cisgender participants were recruited from the community via advertisements describing the study as being about brain development and gender. Advertisements were posted online (Kijiji, Facebook) and on bulletin boards in the Greater Toronto area (e.g., City of Toronto message boards). Study information was also distributed to cisgender participants via word of mouth to facilitate snowball sampling.

2.2. Procedures

Prospective participants were screened for eligibility in person or over the phone. If eligible, participants between 12 and 15 years old provided verbal assent and a parent/guardian provided informed consent. Participants 16 to 17 years old provided informed consent. After enrolment into the study, about 20 mL of blood was drawn to measure hormones and other indices of physiological functioning, a 1-h MRI session was conducted, and participants completed a brief intelligence assessment and questionnaire package. If GD participants completed any study measures during a clinical assessment, these data were obtained from their clinical record. For most participants, all study procedures were completed on the same day, and the order was allowed to vary. At the end of the study visit, participants were thanked and received an honorarium of $20 CAD per hour or part thereof for their time and were also compensated for travel expenses. Only measures relevant to the aims of this study are reported.

2.3. Measures

All numerical values for variables reported and most of the code used are available from Scholar’s Portal Dataverse [47]. For measures where a mean is calculated, the mean was computed for those participants who responded to at least 75% of the items.

2.3.1. Age

Age in months was calculated by subtracting the date that the consent form was signed from the date of birth and was rounded to the nearest month. The MRI scan was conducted on the same day as consent for 42 participants, within one month for four participants, and within two months for one participant.

2.3.2. Gender Dysphoria

The Gender Identity/Gender Dysphoria Questionnaire for Adolescents and Adults (GIDYQ-AA) [48,49] is a 27-item questionnaire assessing concurrent gender identity and gender dysphoria with good discriminant validity and clinical utility [47,48,49,50]. Cisgender boys completed the male version and cisgender girls and GD AFAB completed the parallel female version. Each item was rated on a 5-point Likert-type scale ranging from 1 (always) to 5 (never) with respect to the past 12 months. An example of a female version item is: “In the past 12 months, have you felt unhappy being a girl?” A mean score was calculated, and lower scores indicated stronger gender dysphoria.

2.3.3. Sexual Orientation

The Erotic Response and Orientation Scale (EROS) is a 16-item self-report measure assessing sexual orientation with regard to attractions and fantasies over the past 6 months with good discriminant validity [51,52,53]. Half of the questions pertain to attractions/fantasies toward boys/men (i.e., androphilia; e.g., “How often have you noticed you had any sexual feelings [even the slightest] while looking at a boy/man?”). The other half used similar language to assess attractions/fantasies toward girls/women. Frequency of occurrence for each item was rated on a 5-point scale ranging from 1 (not at all) to 5 (almost every day). Mean androphilia and mean gynephilia scores were derived for each participant, where higher scores reflected more attractions/fantasies. Internal consistency on both scales was high (Cronbach’s alpha: androphilia = 0.94, gynephilia = 0.95).

EROS scores measure two aspects of sexual orientation: strength of attractions (i.e., none to many attractions, regardless of orientation) and degree of androphilia-gynephilia. To characterize these two aspects, two metrics were created for each participant. A vector was created by calculating the strength (magnitude) and degree (phase, θ) of sexual attractions from the individual androphilia and gynephilia scores, which were not significantly correlated with each other (see Results). Androphilia was arbitrarily placed on the x-axis and gynephilia on the y-axis. Based on EROS scores (ranging from 1 to 5 on each scale), the strength (magnitude) of sexual attractions was denoted by the length of the vector, calculated for each participant using the following formula:

$$\sqrt{({\left(androphilia score\right)}^2+ {\left(gynephilia score\right)}^2)}$$

Degree of androphilia-gynephilia was calculated for each participant using the following formula:

$$\left(\mathit{arccos}\left(\frac{androphilia score}{magnitude}\right)\right)\ast \left(\frac{180°}{\pi }\right)$$

These values are symmetrically distributed ±34° around 45° such that an 11° phase indicates exclusive androphilia, a 79° phase indicates exclusive gynephilia, and a 45° phase indicates equal scores on androphilia and gynephilia (e.g., asexual or ambiphilic) (Figure 1).

2.4. Magnetic Resonance Imaging (MRI) Methods

2.4.1. Image Acquisition

All participants were scanned at the Centre for Addiction and Mental Health in Toronto, Ontario, Canada on a GE MR750 3T magnetic resonance scanner (General Electric, Milwaukee, WI, USA) with an 8-channel head coil (General Electric, 8HR BRAIN, GE Standard 8-Channel Head Coil). T1w images were acquired with a T1 BRAVO pulse sequence in the sagittal plane. The acquisition parameters were: inversion time = 650 ms, echo time = 3 ms, repetition time = 6.8 ms, flip angle = 8°, field of view = 23 cm, 256 mm × 256 mm matrix, 200 isotropic 0.9-mm thick slices, acquisition time = 4:42 min. Calculation of the T1 maps with B1 correction required four acquisitions using sagittal spoiled-gradient echo (SPGR) pulse sequences. These scans included: two fast-SPGR scans (repetition time = 10.6 ms, echo time = ~4.4 ms, acquisition time = 2:59 min each) with 1 mm isotropic resolution (field of view = 25.6 cm, 160 slices) at two flip angles (14° and 3°). These two images will be referred to as “Flip14” and “Flip3,” respectively, and were used for T1 mapping calculations using the variable flip angle method [54,55]. Also, to account for B1 inhomogeneities, B1 maps were computed using an extrapolation to signal null as per the method of slopes [56] from two SPGR scans (repetition time = 50–60 ms, echo time = 5 ms, acquisition time = 2:24 min each) with low, 4-mm isotropic resolution (field of view = 25.6 cm, 40 slices) and two flip angles (130° and 150°). Total acquisition time for the four scans was 9:46 min.

2.4.2. T1w Image Processing

T1w images were pre-processed using the minc-bpipe-library on SciNet [57,58]. Briefly, images were realigned to roughly align with the MNI template, signal intensity non-uniformity was corrected [59], and the brain was extracted [60] (see also [61,62]). Extensive quality checking of images for intensity non-uniformity, motion, and quality of brain extraction was performed after pre-processing. The extracted brains were edited manually to remove obvious segmentation errors.

Next, the bias-corrected and brain-extracted T1w images were submitted to the CIVET 2.1.0 processing pipeline [63,64] on SciNet [57,58] to estimate cortical thickness and surface area. Briefly, CIVET extracts a white matter surface (the boundary between cortical gray matter and subcortical white matter) and a pial surface (the boundary between cortical gray matter and cerebral spinal fluid) containing 40,962 vertices per hemisphere. Cortical thickness (mm) is calculated as the distance between the white matter and pial surface at each vertex [65], and surface area (mm²) is calculated at the mid-point between the two surfaces using a Voronoi parcellation based on the 6 triangles adjacent to each vertex [66]. A smoothing kernel of 20 mm and 40 mm was applied to statistically regularize the cortical thickness and surface area data, respectively, and the quality of the surface extractions was verified. Mean surface area and mean cortical thickness for 76 regions of interest (ROI) (38 per hemisphere) were computed based on the Anatomical Automatic Labelling (AAL) atlas [67].

2.4.3. T1 Relaxation Time Map Creation

Cortical thickness and surface area are derived from T1w images, but additional, quantitative information regarding subtle cortical microstructural characteristics can be obtained using methods to calculate tissue T1. T1 reflects the rate at which magnetic spins, in a static magnetic field, lose energy and return to their original equilibrium state after absorbing energy by a perturbing time-varying magnetic field. The magnetic spins of interest in brain MRI are the hydrogen nuclei of water molecules and hence, T1 reflects the rate of energy loss, which depends on the interactions of these water molecules with surrounding molecules, including other water molecules, lipids, proteins, and other macromolecules. T1 maps (and extracted T1 values) reflect those interactions at the voxel level and can thus reveal microstructural differences across tissue types and brain regions. Specifically, T1 is shorter in white matter compared with gray matter due to the relative larger numbers of other molecules and fewer free-water molecules (i.e., water that is not in the myelin sheath) [43]. T1 in gray matter becomes shorter with increasing age [44,46], including during adolescence [17,18,45].

All four SPGR images used to create the T1 maps were reoriented into the space defined as halfway between the Flip3 and Flip14 images using FLIRT (FMRIB’s Linear Image Registration Tool) from the FSL (FMRIB Software Library, version 5.0.9) tool library. Specifically, the “halfway flirt” command from FSL’s SIENA pipeline was used. The Flip3 image was brain-extracted using FSL’s BET (Brain Extraction Tool) and the resulting mask was used to extract the brain of the Flip14 image. T1 maps were then calculated and calibrated as described in Chavez [68]. In this procedure, B1 maps were generated using the method of slopes to obtain the true flip angle [56] and T1 maps were computed using the variable flip angle method with a B1 correction incorporated [55].

2.4.4. T1 Map Processing Pipeline

To extract the ROI data from the T1 maps, the reoriented Flip14 image was registered to the brain-extracted T1w image from the bpipe pipeline. The resulting transformation matrix was used to move the T1 map into T1w space using FSL’s FLIRT. Brain-extracted T1w images were registered to the 1-mm MNI template using antsRegistration (version 2.1.0; [69,70]). The resultant matrix was then inverted and used to move the AAL template from MNI space to T1w space using ANTs’ antsApplyTransforms with the generic label option. Each registration step and transformed image was inspected for registration errors.

To better match the cortical thickness data, T1 data were extracted only from gray matter voxels in each ROI. FSL’s Automated Segmentation Tool (FAST) was used to segment the brain-extracted T1w images from the bpipe pipeline into gray matter, white matter, and cerebral spinal fluid, allowing partial volume estimates and disabling the bias field correction. Using an in-house MATLAB (2017b, version 9.3.0.713579 [71]) script, T1 data were extracted for the same 76 AAL atlas ROIs as were used in CIVET. T1 data were included only for voxels where gray matter probability was > 0.40 (i.e., the voxel contained at least 40% gray matter). A histogram of T1 values was plotted for each ROI using T1s 500 ms ≤ T1 ≤ 5000 ms and a bin width of 30 ms. The distributions were generally unimodal (some ROIs included white matter which led to a small peak at lower T1), but were often positively skewed (partly due to partial volume effects between gray matter and cerebral spinal fluid). To avoid overestimating T1 in the ROI, the T1 mode, rather than the mean, was used to represent the T1 value in the ROI. This was accomplished by fitting a spline function to the histogram and finding the T1 value at the peak (see Figure S1 in Supplementary Materials). Histograms of T1 values were spot checked for anomalies.

3. Statistical Analyses

3.1. Group Differences in Demographic and Psychosexual Variables

Using SPSS version 27 (IBM Corp., 2020), group differences for age, the GIDYQ-AA (gender dysphoria), and EROS (sexual orientation) variables were examined with one-way analyses of variance (ANOVA). In the presence of a significant omnibus effect, post hoc comparisons were conducted with least significant difference (LSD) tests. The GIDYQ-AA had a significant Levene’s test, so the Games-Howell post hoc test, robust to heterogeneity of variance, was used for this measure. A two-tailed critical p-value of 0.05 was used.

We also investigated ethnicity, parent education, parent marital status, subtests on the Weschler intelligence scales capturing verbal comprehension and visual spatial indices [72,73,74], pubertal development [75], regularity of menstrual cycle, medication use, externalizing mental health challenges, internalizing mental health challenges, and total mental health challenges [76,77,78]. Tests of group differences for these additional demographic variables are provided in the Supplementary Materials. Also, correlations between these variables and the ROI data are presented, corrected for multiple comparisons. No demographic variable correlated with the ROI data after correction, and age, ethnicity, parent education, verbal comprehension index, visual spatial index, pubertal development, regularity of menstrual cycle, and externalizing mental health challenges did not show significant group differences. Thus, apart from including age as a variable of interest in specific analyses, no other demographic variables were included in the main analyses.

Total brain volume was not included as a covariate because of known sex assigned at birth differences in total brain volume [79], which have been shown to impact detectability of sex differences in cortical thickness [80]. Given that we were interested in examining how GD AFAB adolescents are positioned relative to likely differences between cisgender girls and boys, controlling for total brain volume could overly remove informative sex differences. Furthermore, not all studies have controlled for total brain volume when investigating cortical thickness and surface area (e.g., [13]). Raznahan et al. [13] demonstrated that changes in cortical volume are reflected more in changes in surface area rather than cortical thickness changes, making controlling for total brain volume most relevant for surface area. However, in a recent large-scale analysis of sex differences including surface area, controlling for total brain volume did impact sex differences in cortical thickness and did not impact sex differences in surface area [15]. Given the mixed impact that controlling for total brain volume has had on the detectability of sex differences in surface area and cortical thickness, we chose not to apply this correction.

3.2. Partial Least Squares Analyses

The multivariate statistical technique partial least squares (PLS, [81,82]) was used to examine the relationship between the ROI data and group in two analyses. PLS works well with data that are not independent from each other, as well as in situations where the number of observations is larger than the number of participants, both of which are true for brain data. Furthermore, PLS identifies distributed patterns in brain data, including structural MRI data [83]. PLS uses singular value decomposition of a data matrix to produce latent variables with three components: group and/or task contrasts (v), brain saliences/weights (u), and a measure of the strength of that relationship (s). We were primarily interested in identifying the commonalities and differences among the brain metrics as they related to group differences and to group differences in correlation patterns (rather than identifying differences across metrics). To achieve these goals, the data matrix was created as a single row for each participant containing 76 cortical thickness measures, 76 surface area measures, and 76 T1 measures. The results will be referred to by ROI or region “type” (e.g., cortical thickness/surface area/T1). We noted that the distribution of raw surface area values was not normal, but they were for cortical thickness and T1 data. To minimize the impact of the differences in data distributions as well as differences in metric-specific ranges of values on the results, all brain data were normalized by converting them to z-scores within each ROI (i.e., across all participants; e.g., [84]).

Permutation tests were used to assess the number of times the strength of the permuted latent variables exceeded the observed strength and provided an exact probability. Bootstrap resampling was used to estimate the standard error of each brain salience, to assess the reliability of their contribution to the observed pattern. The ratio of the salience to its standard error approximates a z-score, and for all analyses a threshold of ± 3 was used to identify stable, reliable contributions. Confidence intervals (CI, 95%) around the point estimates were calculated using the bootstrap sampling distribution. For all analyses, 1000 permutations and bootstrap samples were used. For all analyses, “brain scores” were calculated. These scores are the projection (dot-product) of the brain saliences (u) with a participant’s data and provide an indication of how strongly each individual reflects the contrast identified on the latent variable.

The first PLS analysis was a “Task PLS,” which in this case, distinguished the three groups based on an optimal combination of surface area, cortical thickness, and T1. Age was not included in this analysis given no group differences in age and given no significant correlations between age and the ROIs after correcting for multiple tests (see also Results and Supplementary Materials). The second analysis was a “Behavior PLS,” which identified the commonalities and differences in brain-behavior correlations across the three groups. For this Behavior PLS analysis, the behavior variables were the two EROS variables (i.e., strength of attractions and degree of androphilia-gynephilia) and age, and the analysis assessed the degree to which they correlated with the brain measures. Age was included because the strength of attractions scores correlated with age (see Results). For the Behavior PLS analysis, the correlations of brain scores with behavior measures, which reflect the observed contrast (v), were used [82]. Both PLS analyses were conducted using PLS software (v. 6.1311050) using MATLAB (2014a, version 8.3.0.532, [85]). Correction for multiple comparisons is not required because statistical significance is only evaluated in relation to the strength (s) of the latent variable [82].

4. Results

4.1. Demographic and Psychosexual Variables

Descriptive statistics by group are shown in

Table 1. Descriptive statistics for demographic and psychosexual variables.

	Cisgender Boys	GD AFAB	Cisgender Girls	F (df)	p
n	14	16	17
Age (months)
M	184.93	193.38	191.71	0.75 (2, 44)	0.480
SD	25.61	14.20	19.00
Range	147–216	162–216	152–214
GIDYQ-AA
M	4.91	2.21	4.90	693.34 (2, 44) c	<0.001
SD	0.12	0.35	0.15
Range (1–5) a	4.63–5.00	1.74–3.04	4.48–5.00
Strength of attractions
M	3.25	3.00	2.83	0.54 (2, 44) d	0.587
SD	1.28	1.30	0.76
Range (1.41–7.07) b	1.70–5.37	1.41–6.05	1.41–4.44
Degree of androphilia-gynephilia
M	64.60	46.38	30.94	25.93 (2, 44)	<0.001
SD	9.29	15.23	13.18
Range (11–79) b	43.96–75.55	17.35–75.55	14.04–68.20

Note. GIDYQ-AA = Gender Identity/Gender Dysphoria Questionnaire for Adolescents and Adults; GD = gender dysphoria; AFAB = assigned female at birth. ^a Absolute range. ^b Possible range for magnitude and phase (θ) based on EROS scores with absolute range of 1 to 5. ^c A significant Levene’s test (p = 0.009) warrants reporting of the robust tests of equality of means via the Welch test statistic (2, 27.46) = 422.50, p < 0.001, and the Brown–Forsythe test statistic (2, 24.31) = 713.80, p < 0.001. ^d A significant Levene’s test (p = 0.032) warrants reporting of the robust tests of equality of means via the Welch test statistic (2, 25.76) = 0.60, p = 0.557, and the Brown–Forsythe test statistic (2, 35.89) = 0.52, p = 0.599.

. There were no extreme deviations from normality based on skewness and kurtosis values, which were less than |2|.

For age, the main effect of group was not significant (see Table 1). As expected, there was a significant effect for the GIDYQ-AA (gender dysphoria) such that the GD AFAB group scored significantly lower than both the cisgender girls (mean difference (MD) = −2.69, SE = 0.10, p < 0.001) and cisgender boys (MD = −2.70, SE = 0.09, p < 0.001), who did not differ from each other (MD = 0.01, SE = 0.05, p = 0.970). The threshold for a potential diagnosis of GD (“caseness”) on the GIDYQ-AA is ≤3.00 [50]. The lowest score on the GIDYQ-AA from the cisgender participants was 4.48, indicating that none of the cisgender participants met this threshold; the highest score in the GD AFAB group was 3.04, with all other GD AFAB participants scoring below the threshold.

Regarding sexual orientation, a main effect of group was found for degree of androphilia-gynephilia, whereas strength of attractions showed no significant group differences (Table 1). For degree of androphilia-gynephilia, cisgender boys were more gynephilic than both GD AFAB (MD = 18.21, SE = 4.74, p < 0.001) and cisgender girls (MD = 33.66, SE = 4.67, p < 0.001). GD AFAB were more gynephilic than cisgender girls (MD = 15.44, SE = 4.51, p = 0.001). Examination of the scores showed that cisgender boys tended to be gynephilic, cisgender girls tended to be androphilic, and GD AFAB had a range of sexual attractions with a cluster of individuals in the ambiphilic or asexual range (see Figure 2).

Androphilia scores were not related to gynephilia scores across participants (r = −0.09, p = 0.570, n = 47), within cisgender boys (r = 0.43, p = 0.129, n = 14), within GD AFAB (r = 0.25, p = 0.360, n = 16), or within cisgender girls (r = −0.16, p = 0.553, n = 17). Likewise, strength of attractions was not related to degree of androphilia-gynephilia across participants (r = 0.05, p = 0.765, n = 47), within cisgender boys (r = 0.24, p = 0.418, n = 14), within GD AFAB participants (r = −0.15, p = 0.570, n = 16), or within cisgender girls (r = −0.35, p = 0.163, n = 17). Strength of attractions was significantly correlated with age across all participants (r = 0.53, p < 0.001, n = 47), within cisgender boys (r = 0.87, p < 0.001, n = 14) and cisgender girls (r = 0.58, p = 0.014, n = 17), but not within GD AFAB (r = 0.20, p = 0.458, n = 16). Degree of androphilia-gynephilia was not significantly correlated with age across all participants (r = 0.01, p = 0.940, n = 47), nor within any group: cisgender boys (r = 0.39, p = 0.168, n = 14), GD AFAB (r = 0.24, p = 0.370, n = 16), and cisgender girls (r = −0.01, p = 0.958, n = 17).

4.2. Task Partial Least Squares (PLS) Analysis

The Task PLS analysis identified one significant latent variable (p = 0.000) that accounted for 86.59% of the covariance between group and the ROI data. This latent variable differentiated the cisgender girls and GD AFAB from the cisgender boys. Thus, the overall pattern reflected a sex assigned at birth difference (Figure 3, Panel A). The majority of brain saliences and all stable brain saliences were negative, indicating that cisgender boys had generally thicker cortices and larger surface areas than GD AFAB and cisgender girls. T1 did not reliably contribute to the sex difference. The sex difference was most stably expressed in 25 regions, 5 reflecting cortical thickness and 20 reflecting surface area differences (Figure 3, Panel B). These stable regions were distributed throughout the lobes of the brain, and included primarily frontal, temporal, and occipital ROIs (

Table 2. Stable ROIs contributing to group differences from Task PLS and Behavior PLS.

Task PLS			Behavior PLS
SA/CT	Hem.	ROI	SA/CT/T1	Hem.
Frontal Lobe
CT	Left	Superior Frontal Gyrus: Orbital Part	SA	Left
CT	Right	Superior Frontal Gyrus: Orbital Part
SA	Left	Middle Frontal Gyrus	T1	Left
CT	Left	Rolandic Operculum	T1	Left
SA	Left	Superior Frontal Gyrus: Medial
SA	Left	Precentral Gyrus	T1	Left
SA	Right	Precentral Gyrus
CT	Right	Gyrus Rectus
		Gyrus Rectus	SA	Left
SA	Left	Supplementary Motor Area
		Middle Frontal Gyrus: Orbital Part	T1	Left
		Middle Frontal Gyrus: Orbital Part	SA	Right
		Superior Frontal Gyrus: Dorsolateral	T1	Left
		Inferior Frontal Gyrus: Opercular Part	T1	Left
		Inferior Frontal Gyrus: Triangular Part	T1	Left
		Inferior Frontal Gyrus: Orbital Part	T1	Left
Parietal Lobe
SA	Left	Angular Gyrus	T1	Left
		Angular Gyrus	CT; T1	Right
		Supramarginal Gyrus	SA	Left
		Supramarginal Gyrus	CT	Right
		Postcentral Gyrus	T1	Left
		Precuneus	T1	Left
		Precuneus	T1	Right
Temporal Lobe
CT	Left	Parahippocampal Gyrus
SA	Left	Temporal Pole: Superior Temporal Gyrus
SA	Left	Superior Temporal Gyrus	T1	Left
		Superior Temporal Gyrus	T1	Right
SA	Left	Middle Temporal Gyrus	T1	Left
		Middle Temporal Gyrus	T1	Right
SA	Right	Inferior Temporal Gyrus	T1	Right
		Temporal Pole: Middle Temporal Gyrus	CT	Left
		Fusiform Gyrus	T1	Left
		Heschl Gyrus	T1	Left
		Heschl Gyrus	T1	Right
Occipital Lobe
		Middle Occipital Gyrus	T1	Left
SA	Right	Middle Occipital Gyrus	T1	Right
		Inferior Occipital Gyrus	T1	Left
SA	Right	Inferior Occipital Gyrus	T1	Right
SA	Left	Lingual Gyrus	SA; T1	Left
		Lingual Gyrus	T1	Right
SA	Left	Calcarine Fissure and Surrounding Cortex
SA	Left	Cuneus	T1	Left
		Cuneus	T1	Right
		Superior Occipital Gyrus	T1	Left
		Superior Occipital Gyrus	T1	Left
Insula and Cingulate Gyri
SA	Left	Insula
SA	Right	Insula
SA	Right	Anterior Cingulate and Paracingulate Gyri
SA	Left	Posterior Cingulate Gyrus
SA	Right	Posterior Cingulate Gyrus	T1	Right
		Median Cingulate and Paracingulate Gyri	T1	Right

Note. PLS = partial least squares, SA = surface area, CT = cortical thickness, T1 = T1 relaxation time, Hem. = hemisphere, ROI = region of interest.

; Figure 3, Panel C).

4.3. Behavior PLS Analysis

The Behavior PLS analysis identified one significant latent variable (p = 0.000) that accounted for 45.87% of the cross-block covariance between the ROI data and the two EROS variables and age across the groups. This latent variable indicated that cisgender boys and GD AFAB shared similar brain–behavior correlations, with no stable correlations in the cisgender girls. For the cisgender boys, EROS strength of attractions, EROS degree of androphilia-gynephilia, and age were stably related to the ROI data; for the GD AFAB group, similar correlations were found for EROS degree of androphilia-gynephilia and age (Figure 4, Panel A). The majority of brain saliences were positive, indicating that the brain metrics were inversely related to the behavioral measures. The stable regions were dominated by those reflecting correlations of T1 with the behavioral measures, across the brain. Thus, in both cisgender boys and GD AFAB, shorter regional T1s correlated with older age, stronger gynephilia, and—for cisgender boys only—stronger attractions. Across all measures, the brain-behavior correlations in cisgender boys and GD AFAB were most stably expressed in 41 regions, 33 of which reflected regional correlations with T1. The regions contributing to the observed brain-behavior patterns in cortical thickness and surface area regions were fewer, with only 3 stable contributions in cortical thickness regions, and 5 stable contributions in surface area regions (Figure 4, Panel B). With the exception of one cortical thickness region, the brain saliences were positive, reflecting negative correlations with the behavioral measures for both cortical thickness and surface area, similar to what was found for T1. As was found in the Task PLS results, the stable regions were distributed throughout the lobes of the brain (Table 2; Figure 4, Panel C).

5. Discussion

Using multivariate PLS analyses, we examined group differences (GD AFAB, cisgender boys, cisgender girls) and associations between sexual attraction, age, and three structural brain metrics: cortical thickness, surface area, and T1 relaxation time. The Task PLS results indicated the structural metrics in GD AFAB adolescents were similar to cisgender girls, and both groups differed from cisgender boys. This sex assigned at birth difference was driven mostly by surface area, with some cortical thickness contributions, such that cisgender boys had larger surface areas and thicker cortices than GD AFAB and cisgender girls. Thus, when examining group differences, a sex difference was revealed, rather than a gender identity difference, which did not support our prediction that the brains of GD AFAB would be similar to their experienced gender. However, our prediction that a sex difference would be driven mostly by surface area was supported.

Many studies, including those with cisgender adolescents, show a sex assigned at birth difference in surface area such that AMAB individuals have larger surface areas than AFAB individuals [13,14,15,36]. Therefore, it is not surprising that surface area contributed the most to the sex difference in the current study. The stable surface area ROIs contributing to this sex difference were distributed across the brain, with most in the frontal, temporal, occipital, and insula/cingulate gyri regions, again aligning with earlier studies of widespread sex differences in surface area [15]. The more limited contribution of cortical thickness to the sex difference was predominant in fronto-temporal ROIs. This aligns with at least one large-scale study that included adolescents and found regions in the frontal and temporal lobes where males had thicker cortices than females, although the regions did not overlap completely with those found in the current study [15]. Not all regions in Wierenga et al. [15] showed thicker cortex in males (vs. females) and there have been mixed results regarding adolescent sex differences in cortical thickness [13,14]. T1 did not stably contribute to the sex difference, consistent with early studies that did not find sex differences in regional T1 [16,17,18].

Our finding that group differences in the Task PLS were related to sex assigned at birth partly aligns with the only other study on the cortical structure of the GD adolescent brain. Hoekzema et al. [5] found that cortical and subcortical gray matter volume in GD AFAB and GD AMAB adolescents not receiving hormone therapy mostly reflected sex assigned at birth and not gender identity, with a few exceptions. Specifically, in a small number of regions displaying sex differences, gender identity differences were also found: in GD AFAB, a smaller left superior medial frontal cortex volume vs. cisgender girls, and in GD AMAB, a smaller volume in the left superior posterior hemisphere of the cerebellum and marginally smaller hypothalamic volumes vs. cisgender boys. We did not find any exceptions to our pattern of sex differences in the Task PLS and in that regard, our study does not fully align with Hoekzema et al. [5].

Some insights regarding mechanisms underlying the development of surface area and cortical thickness have been advanced. The radial unit hypothesis partly explains the prenatal cellular contributions to regional variations in surface area and cortical thickness [87,88]. Specifically, surface area is represented by the number of neuronal columns in an area, determined during in utero brain development by neurons migrating along radial glial cells from the ventricular zone to the cortical plate. Cortical thickness is determined by the number of neurons within these columns. There is evidence from imaging studies that somatic and maturational variables (i.e., age, sex, and physical growth) are more strongly associated with surface area and cortical thickness in the adolescent brain than variables reflecting relatively more environmental influences (e.g., general intelligence, prenatal parental smoking; [89]). If so, it might explain why group differences in the present study were reflected in birth-assigned sex differences in cortical structure, particularly surface area. Imaging studies have further demonstrated that surface area is influenced more by early prenatal environment factors than is cortical thickness [40,41]. The exact mechanisms have not been elucidated, but given the likely prenatal basis for the development of cortical thickness, and particularly surface area, prenatal organizational mechanisms may play a role.

The neurohormonal hypothesis argues that prenatal androgens shape the brain to be organized in a particular manner, including organization of sex differences in the brain. Then, during adolescence, the pubertal surge in sex hormones further influences brain regions to be expressed in a relatively male- or female-typical manner, based on the earlier brain organization [12]. Our finding of a sex assigned at birth difference within the cisgender adolescents aligns with the neurohormonal hypothesis. Thus, variation in prenatal androgen exposure could be one mechanism that produces larger surface areas and somewhat thicker cortices in individual AMAB versus AFAB. However, despite a likely prenatal basis for the development of these two structural brain features, they continue to develop throughout childhood and adolescence [13,35,36,37]. Thus, additional mechanisms, which are not necessarily mutually exclusive, may play a role. For example, circulating pubertal hormones could influence neural reorganization [90]. The pubertal increase in steroid sex hormones and accompanying physical maturation have been related to cortical development in cisgender adolescents [91]. Furthermore, the effects of these physical changes on the brain are speculated to be influenced further by indirect mechanisms such as interactions with the social environment [92].

In contrast to the group difference associated with sex assigned at birth, when brain structure was assessed in relation to sexual orientation and age, the findings indicated similarities between GD AFAB and cisgender boys that were not shared with cisgender girls. In the Behavior PLS analysis, we found that cisgender boys and GD AFAB adolescents had similar correlations between brain measures and age and sexual orientation. Specifically, gynephilia, older age, and stronger attractions (in cisgender boys only) were related predominantly to shorter T1 in regions distributed throughout the frontal, temporal, parietal, and occipital lobes. Thus, in GD AFAB and cisgender boys, brain patterns showed similarities in how they varied in relation to sexual orientation as well as developmentally (i.e., in relation to age). Given these patterns were not seen in cisgender girls, the similarities between GD AFAB and cisgender boys in the context of sexual attractions and age provide some support for our prediction that the brains of GD AFAB would be similar to their experienced gender.

For T1, 33 regions stably contributed to the similarities between GD AFAB and cisgender boys such that T1 was shorter in older GD AFAB and cisgender boys. Shorter T1s are associated with more myelination and other macromolecules, more iron, less free water, and increased gray matter tissue density [43]. T1 typically shortens in gray matter with age, including during adolescence [17,18,43,44,45,46]. The adolescent brain also experiences changes in surface area and cortical thickness, with some minor growth in surface area and decreases in cortical thickness to allow for the increase in white matter (without much growth in brain volume) during this time [13,35,36,37,42]. In our sample, only five surface area regions and three cortical thickness regions combined with the globally shorter T1s to contribute to the observed similarities in correlations between brain metrics and sexual orientation measures and age in GD AFAB and cisgender boys. Thus, one way in which the adolescent brains of GD AFAB are similar to those of cisgender boys could be through a similarity in timing of age-related changes in mostly cortical tissue microstructure related to sexual attraction.

The results of the Behavior PLS suggest that only focusing on cortical thickness and surface area may overlook insights provided by mechanisms that influence variation in cortical tissue microstructure. These mechanisms are not fully understood, but findings related to the development of dendritic arborization and myelin could provide insights (e.g., [93]). Given associations of the brain metrics with sexual orientation and age, mechanisms that influence these variables may also provide a good starting point for speculation. There is some support for the neurohormonal hypothesis in underlying the development of sexual orientation, particularly in same-sex attracted cisgender women (i.e., gynephilia) [23,30,94]. Our results show that cortical structure and microstructure of GD AFAB is similar to cisgender boys in the context of gynephilia, which is consistent with the neurohormonal hypothesis and the role of prenatal androgens. However, alternative explanations are also plausible. For instance, if the cortical microstructure similarities are attributable to myelination [95], it has been demonstrated in rodent models that estradiol decreases myelination during puberty [96], whereas androgen exposure leads to thicker myelin sheaths [97]. However, these cellular properties have not yet been related to any psychological traits in humans (e.g., gynephilia).

Shorter T1 was related to stronger attractions in cisgender boys only, highlighting the likely developmental importance of this aspect of sexual attractions. Development of the strength of attractions, and mechanisms underlying the development of this aspect of sexual orientation, are understudied in the sexual orientation literature. Xu et al. [34] recently found that the intensity of sexual behaviors increased (i.e., from hugging to having sexual intercourse) with age in cisgender adolescents, and the likelihood of engagement in same-sex behaviors increased in later adolescence, but these relationships were found in both cisgender boys and cisgender girls, strength of sexual attractions was not examined, and the findings were not discussed further. Given that cisgender boys showed strong relationships between tissue microstructure, degree of androphilia-gynephilia, strength of sexual attractions, and age, incorporating similar measures in future studies of gender and sexual orientation in adolescence is warranted.

Limitations

The present results should be interpreted with caution in light of criticisms of neuroimaging research conducted with small sample sizes [98]. We ameliorated some concerns of our small sample size by utilizing a multivariate statistical approach and by interpreting distributed patterns rather than focusing on individual brain regions. Furthermore, our results largely align with the only other study of adolescent GD brain structure, which had a similar sample size of GD AFAB adolescents not on puberty blockers or gender-affirming hormone therapy, but did have larger samples of cisgender adolescents [5].

An additional limitation is that the cross-sectional non-experimental design does not allow us to determine causal relationships. With our study design, we cannot infer whether changes in the measured brain features result in GD and/or sexual orientation, or vice versa. Moreover, we cannot know whether changes in the measured brain features and GD and/or sexual orientation are the result of a third variable (e.g., a shared underlying cause such as variation in prenatal hormone exposure). We also used a convenience sample, which may affect representativeness. We did not directly examine other non-brain mechanisms (e.g., peripheral hormone levels) to interpret the results. Finally, we did not measure consumption of various substances (e.g., water, caffeine, smoking) prior to the MRI scan, or quantify differences in vasculature and metabolism (e.g., cerebral blood flow) that have been shown to impact measurements of T1 and cortical thickness (e.g., [99,100,101]).

Our results can only address GD AFAB adolescents given GD AMAB adolescents were not included. Also, our results cannot generalize to the broader sexual orientation literature given the small number of cisgender same-sex attracted adolescents in this sample. The restricted range in degree of androphilia-gynephilia in cisgender boys and cisgender girls may have also impacted our findings in the Behavior PLS. Furthermore, our results cannot generalize to ages beyond 12- to 17-years-old. Studies on adults have examined cortical structure related to GD [9,26,102,103,104,105,106,107,108], but we know of no studies examining cortical structure of GD children. Lastly, we examined three cortical structural brain features; other aspects of adolescent GD brain structure (e.g., white matter microstructure) await investigation.

6. Conclusions

In our sample of GD AFAB adolescents not receiving hormone therapy, cisgender girls, and cisgender boys, a multivariate sex assigned at birth difference in cortical structure was identified, driven mostly by variability in surface area. When assessing correlations of cortical structure with sexual orientation and age, the brains of GD AFAB adolescents were similar to cisgender boys, driven mostly by T1 relaxation time. Our findings suggest that the cortical structure of GD AFAB adolescents aligns with the experienced gender in the context of age-related changes in sexual attractions during adolescence reflected mostly in tissue microstructure.