Scientific attention allocated to topics concerning women enabled researchers to connect endocrine activity to anxiety. Indeed, paradigms using naturally-fluctuating hormones (e.g., estradiol [E2] and progesterone [P]) or synthetic hormones (e.g., oral contraceptives [OCs]) have revealed that sex hormones modulate anxiety behaviors (for reviews see 16, 21, 22). Endogenous sex hormones are the end result of the endocrine cascade of the hypothalamic-pituitary-gonadal (HPG) axis and can fluctuate according to the menstrual cycle. That said, used by more than 150 million women worldwide (23), OCs act as a powerful disruptor of the menstrual cycle. Combined OCs (COCs) are the most common type and are made up of synthetic estrogen (e.g., ethinyl estradiol [EE]) with a progestin (e.g., levonogestrel). These exogenous hormones bind to estrogen and P receptors in brain regions involved in the HPG axis, thereby suppressing sex hormone secretion via negative feedback (24–26). Additionally, EE induces hepatic synthesis of the sex hormone binding globulin (SHBG), reducing not only free (unbound) E2 and P levels but also free testosterone (T) (27).

Notably, lowered endogenous E2 levels (either naturally occurring during the early follicular phase of the menstrual cycle or induced by OC use) is a key contributor to fear maintenance as evidenced using fear conditioning and extinction protocols (28–30). At the neural level, sex hormones are known to modulate the brain network implicated in fear processes (28, 31–34). Major brain regions of the fear circuitry include the amygdala, hippocampus, hypothalamus, anterior insular cortex (AIC), dorsal anterior cingulate cortex (dACC), rostral anterior cingulate cortex (rACC), and ventromedial prefrontal cortex (vmPFC) (32, 35–40). Particularly, the amygdala, hippocampus, and vmPFC express high levels of sex hormone receptors and as such, provide insight into how sex hormones exert their influence on fear (6, 28, 41, 42).

Fluctuations in sex hormones have been shown to influence brain activity of the fear circuitry (for reviews see 28, 31, 32). Human neuroimaging studies have associated high E2 with increased activations within this network (33, 43, 44), though fewer studies have investigated the role of sex hormones on its anatomy. The most consistent observation of the latter is the trophic anatomical effects of E2 in the hippocampus (45–48). Although less consistently reported, morphologic changes of fear-related regions have been reported throughout the menstrual cycle (e.g., changes in the dACC (49), AIC (50), and amygdala (51)). OC use has also been associated with a reduction of the hippocampus (52–54), amygdala (55), hypothalamus (56), and prefrontal cortical thickness (57, 58). Additionally, ‘acute’ sex hormone concentrations have been shown to correlate with brain structures (59). For instance, fear-related brain regions have been associated with salivary or serum E2 and T in various hormonal profiles (i.e., OC users, naturally cycling [NC] women, men) (49, 53, 57). Given the scarce literature on the topic and lack of replicability, these associations call for further empirical support. Nevertheless, these findings collectively support the involvement of endogenous and exogenous sex hormones on key regions implicated in fear processing and regulation, as well as their ability to engender significant structural alterations within this network.

The impact of OCs on brain morphology has mainly been studied by comparing current users to NC women (according to specific phases of the menstrual cycle or not) and men (for reviews see 59–61). One of the most important methodological limitations in this type of design is the sample selection bias. Indeed, a ‘survivor effect’ may linger in samples of current OC users, knowing that women for whom the pill does not cause adverse effects are more likely to continue using it (62). Discontinuation of OCs has also been associated with psychiatric symptoms or disorders (63, 64), neuroticism (65), and mood-related side effects (66). Thereby, NC women that are previous OC users could constitute their own group that is distinct from never and current users. However, a recent investigation reported negligible differences between current, past, and never users in sociodemographic variables and personality in a matched Eastern European sample (67). Nonetheless, given that the study of OCs is normally done using cross-sectional designs between users and non-users, careful attention must be paid to the choice of comparison groups.

The vast majority of research on OCs has focused on acute (i.e., ‘here and now’, activational) effects underlying its use. Yet, hormonal events have been shown to exert long-lasting neural changes, as reported in the context of puberty (68, 69), pregnancy (70–72), and menopausal hormonal therapy (73, 74). Regarding OC use, a paucity of studies distinguishing current, past, and never users suggests that OC intake could lead to long-term effects on endogenous sex hormone concentrations (75), brain structures (76, 77), cognition (78, 79), and anxiety symptoms (as experienced during the COVID-19 pandemic, 80). It is still unclear whether such lasting effects exist and if so, by which mechanisms OCs exert their influence after cessation. Still, it appears obvious that the simple comparison of OC users and NC women would generate its share of heterogeneity, as grouping NC women would most likely comprise past and never users.

Generally, studies have homogenized their subset of OC users by selecting COC users only (i.e., excluding progestin-only OC users). If there has been a keen interest over time to refine methodologies to the specific use of COCs, very few studies have considered women’s lifetime history of hormonal contraceptives. Yet, as previously stated, an increasing body of research points to differential effects between current users, past users, and never users (77, 78, 80–82). Knowing that women commonly try several types of hormonal contraceptives in their lifetime, this factor may introduce unwanted noise. Therefore, when examining current and past COC users, rigorous methodological standards can be achieved by considering the history of hormonal contraceptive use. To our knowledge, only two studies have considered contraception history and restricted their sample to women who had only used OCs or COCs throughout their lifetime (76, 83).

Relatedly, COCs is a broad category composed of more than 30 formulations with varying doses and/or synthetic compound compositions (24, 84). Doses of EE range from 10 to 50μg, although COCs containing 35μg of EE or less are the most commonly used. While EE is typically the estrogenic constituent in COCs, progestins (and their dosage) are heterogenous and differ according to their chemical structure and pharmacodynamics (85). Androgenic activity is commonly used to classify progestins and as such, can be classified as high (e.g., levonogestrel), low (e.g., norgestimate, desogestrel), or anti-androgenic (e.g., drospirenone, cyproterone acetate) progestins (24, 86). Of note, androgenicity effects have been examined on cognitive functions (87–91) and neural correlates (79, 82, 83, 92). Anatomically, a larger volume of the bilateral fusiform face area was found in anti-androgenic users compared to androgenic users, while EE dose did not correlate with this region (82). Furthermore, the effects of exogenous hormone concentrations on brain structures have not been reported to date. Synthetic compounds (e.g., EE) are not captured by traditional immunoassays of circulating E2 in serum or saliva, meaning that the current literature has only partly depicted the endocrine activity of OC users (84).

Finally, many metrics to quantify brain morphology have emerged over the years. While this represents a fundamental strength for the field of neuroscience, it has also inevitably led to heterogeneous findings. Structural imaging studies have used different measures to assess brain tissues such as gray or white matter volume, cortical thickness (CT), cortical surface area, and folding/gyrification (61, 93, 94). Voxel-based morphometry (VBM) has traditionally been the most popular method to quantify gray matter volume (GMV), though it can also be measured using surfaced-based morphometry (SBM; as the product of CT and surface area for cortical GMV). However, it is well known that changes in GMV are unspecific and can result in differing interpretations (e.g., changes in CT, folding, or surface area of the cortex) (61, 95). Interestingly, GMV is the most frequent metric used to quantify OC changes in the brain (61) and only one group has investigated the impact of OCs on other gray matter (GM) phenotypes, namely CT (57, 58).

Despite the growing interest in studying women’s health, much remains to be discovered about the impact of sex hormones and OCs on the brain, cognition, behavior, and mental health. In psychoneuroendocrinology, there is a lack of consensus with regard to the effects of sex hormones and OCs, mostly due to the relatively small number of studies and heterogeneous methodologies used (96–102). This study aimed to investigate the current and potentially lasting effects of COC use, as well as the role of endogenous and exogenous sex hormones on structural correlates of the fear circuitry. We recruited women currently using or having used solely COCs in their lifetime, women who never used any hormonal contraceptives, and men. Using structural magnetic resonance imaging (MRI), we extracted GMV and CT of key regions of the fear circuitry. Comparing current users, past users, never users, and men allowed us to examine whether COC use was associated with acute (i.e., current users vs. never users = past users) or long-term (i.e., current users = past users vs. never users) morphologic alterations, as well as to detect sex differences. Additionally, we conducted analyses to explore the associations between endogenous sex hormones and brain morphology in past users, never users, and men, as well as with EE (salivary levels and prescribed doses) and androgenicity in current COC users.

This study stems from a larger research project investigating the neural, cognitive, and endocrine correlates of COC use, which involved three laboratory sessions. Participants were all healthy adults aged between 23 and 35 years old. This minimal age criterion was set to limit brain variability due to late-adolescence maturation (103, 104) without compromising recruitment feasibility (i.e., undergraduates were still allowed to apply). According to their hormonal profile and COC history, participants were divided into groups of men, never users, past users, and current users. Eligibility was confirmed after a thorough phone screening, ensuring that participants from all groups met the following inclusion/exclusion criteria: French speakers, no past or current medical/psychological diagnosis, no current or past pregnancy, no usage of medication affecting the endocrine system (other than COCs), and no regular usage of drugs and alcohol. Participants also completed an MRI screening questionnaire for safety purposes (e.g., no metallic implant, claustrophobia). For COC criteria, past COC users were required to have used their COC and have stopped it for over a year (M_duration = 6.74 years ± 3.21; M_cessation = 3.37 years ± 2.15). Current COC users had to be using it for at least 3 months (M = 8.97 years ± 3.83). Regarding contraception history, past and current users must not have used other types of hormonal contraceptives than COCs in their lifetime (i.e., progestin-only pill, intra-uterine device, patch, vaginal ring, injection, or implant). Never users must not have used any hormonal contraceptives in their lifetime. As for NC women, never users and past users were expected to have a regular menstrual cycle (M = 29.26 days ± 3.01; normal cycle length of 21-35 days, 105). As these women monitored their menstrual cycles, laboratory sessions were scheduled during either the early follicular phase (day 1-5 of the menstrual cycle, where day 1 corresponds to the onset of menses) or the pre-ovulatory phase (1-4 days before ovulation, where ovulation day was calculated as average cycle length minus 14). One never user had a cycle length longer than 35 days and two women (one past user, one never user) were not able to quantify their cycle length, they were therefore scheduled during the early follicular phase. Given that E2 has been found to have a marked influence on the fear circuitry (see section 1.1.), these two cycle phases were selected due to their distinct E2 concentrations. This methodological detail allowed us to investigate the role of cycle phases in NC women, as well as to control for hormone levels within this group of women (i.e., equivalent distributions of cycle phases between both groups).

Data for the present manuscript were collected during the second session of the larger study. Before entering the MRI scanner, saliva samples were collected for later quantification of sex hormone levels. In the MRI scanner, participants first underwent a structural sequence before going through a task for functional sequences. All sessions were scheduled between 8:15AM and 7PM. This study was reviewed and approved by the research ethics board of the Centre intégré universitaire de santé et de services sociaux de l’Est-de-l’Île-de-Montréal. Individuals provided written consent before participating in this study, had the right to terminate the experiment at any time, and were offered monetary compensation for their time. A total of 181 participants completed the procedure for the present study.

The following questionnaires were completed by participants to evaluate potential confounds between our groups and on structural outcomes. Importantly, these constructs have all been associated with COC use and/or brain morphology (65, 106–111). Participants completed the 1) Beck Depression Inventory-II (BDI), a 21-item questionnaire assessing physiological and psychological symptoms related to depression (range 0-63, 112, 113), 2) Brief Trauma Questionnaire (BTQ), a 10-item questionnaire derived from the Brief Trauma Interview (114) that evaluates previous exposure to stressful/traumatic events (115), 3) Neuroticism scale from the NEO-Five Factor Inventory-3, a 12-item scale measuring the tendency towards negative affect and emotional instability (range 12-60, 116), 4) Santa Clara Strength of Religions Faith Questionnaire (SCSRFQ), 10 questions regarding religiosity (range 10-40, 117), and 5) trait form of the State and Trait Anxiety Questionnaire (STAI-T), which measures anxiety as a personality trait through 20 items (range 20-80, 118, 119). All questionnaires (except the BDI) were answered online on Qualtrics XM between the first and second sessions of the larger study. The BDI was completed at the end of the third session and participants scoring ≥ 14 (i.e., mild depression symptomatology) were offered psychological resources.

Salivary sex hormone concentrations were collected using the passive drool method. A 2mL sample was provided in a salivette by all participants and stored temporarily in a -20°C freezer at the MRI facility. Within a period of 2 months, samples were transferred to a -80°C freezer at our research center for long-term conservation. Upon the time of analysis, samples were shipped with dry ice to ZRT Laboratory (Beaverton, Oregon, USA), a lab specialized in salivary assays of sex hormones using liquid chromatography – tandem mass spectrometry (LC-MS/MS). Assays were run with an AB Sciex Triple Quad 5500 system. Saliva was mixed with internal standards, then extraction of steroids was performed by Cl 8 column chromatography. Steroids were eluted from the solid phase extraction and dried under nitrogen. Derivatization was carried out on dried samples (120), after which they were diluted and injected for LC-MS/MS analysis with analytical separation performed on an Agilent Poroshell 120 EC-C8 column and ionization by atmospheric pression chemical ionization (121). The lower limit of quantification (LLOQ) was 0.30 pg/mL for E2, 5 pg/mL for P, 3 pg/mL for T, and 0.40 pg/mL for EE. Given that we had data on which participants were expected to truly have circulating EE levels (i.e., COC users in the active phase of their regimen), we used the lower limit of detection (LOD) of 0.11 pg/mL for EE values below the LLOQ.

Structural images were acquired on a Siemens Magnetom Prisma 3T scan using a 64-channel coil. A T1-weighted Multi-Echo Magnetization-Prepared Rapid Gradient-Echo (MEMPRAGE) sequence was used to collect 176 high-resolution images (TR = 2530ms, TE = 1.62ms, FOV = 176mm, flip angle = 7°, voxel size = 1mm³). Visual inspection for motion artifacts was performed on raw data for each participant based on a 4-point scale protocol (122; https://github.com/CoBrALab/documentation/wiki/Motion-Quality-Control-(QC)-Manual). This procedure resulted in the exclusion of one participant (past COC user).

Raw anatomical images were processed using FreeSurfer (http://surfer.nmr.mgh.harvard.edu). The recon-all pipeline of Freesurfer 6.0 was performed on CBRAIN, an online open-source software (123) and allowed for the automatic generation of a 3D cortical surface model. Segmentation of white matter (WM) and GM (pial) surfaces was visually inspected for every participant. No major flaws were detected by visual inspection (e.g., WM/pial surface segmentation into the skull or ventricles).

We extracted both GMV and CT phenotypes for regions of interest (ROI)-based analyses. The Destrieux atlas (124, 125) was used for the parcellation of cortical ROIs. To select precise parcels of the fear circuitry, we inspected the brain map for the term ‘fear’ on Neurosynth (neurosynth.org), an online platform for the automated synthesis of functional MRI data. Using peak coordinates, we identified which Destrieux’s labels overlapped with fear-related clusters. The dACC was best delimited by the anterior part of the middle cingulate gyrus and sulcus (aMCC, label#7), which was further validated by Vogt et al. (126) as being involved in fear expression. We defined the rACC as the anterior part of the cingulate gyrus and sulcus (ACC, label #6). For the AIC, we combined the short insular gyrus (label #18) and anterior segment of the circular sulcus of the insula (label #47). Finally, the vmPFC was defined as the sum of the straight gyrus (gyrus rectus, label #31) and suborbital sulcus (label #70). At the subcortical level, GMV of the amygdala and hippocampus were obtained from FreeSurfer’s segmentation. The hypothalamus was segmented using the automated tool of Billot et al. (127). GMV of the whole structure was used. All ROIs were extracted by hemisphere. Exploratory vertex-by-vertex whole-brain analyses (WBA) were conducted for cortical GMV and CT using a smoothing full-width half maximum (FWHM) kernel of 10mm and 20mm, respectively.

Additionally, we performed VBM as a secondary analytic method. Data were analyzed using Matlab (R2018b, Mathworks Inc.) and SPM12 (http://www.fil.ion.ucl.ac.uk/spm). Images were reoriented along the anterior-posterior commissure line and the anterior commissure was set as the coordinate origin. Images were segmented into GM, WM, and cerebrospinal fluid (Ashburner and Friston, 2005). The Diffeomorphic Anatomical Registration Through Lie Algebra (DARTEL) algorithm was used for estimating deformations and increasing the accuracy of between-subjects alignment (Ashburner, 2007). Images were modulated and spatially normalized to the Montreal Neurological Institute (MNI) space. Normalized GM images were smoothed with a FWHM kernel of 8 mm. GM voxels with a value < 0.2 were excluded from statistical analyses. Global normalization was set at proportional scaling to account for TIV.

In the present study, we examined two phenotypes of the GM, namely volume and thickness. As volume is an unspecific feature of the GM (95), we also measured CT to refine the interpretation of the results. Further, we used two volumetric methods to generate surface-based and voxel-based GMV. Given the well-documented difficulty of low reproducibility in MRI studies (128–130), this motivated our team to conduct analyses not only in FreeSurfer but in SPM using VBM to evaluate the robustness of analytical variability of GMV results. As such, between-software convergence would strengthen confidence in the obtained results (130). Moreover, given that FreeSurfer’s volumetric WBA is restricted to the cortex, the use of SPM’s VBM allowed us to examine cortical, subcortical, and cerebellar voxels.

Regarding endocrine quantification, the use of LC-MS/MS resulted in a considerable rate of levels falling under the LLOQ. Left-censored undetected values were high for E2 (< 0.3 pg/mL: 110/180, 61.1%), P (< 5 pg/mL: 125/180, 69.4%), and moderate for T (< 3 pg/mL: 51/180, 28.3%). This small proportion of detected values was also reported with LC-MS/MS for salivary E2 (131). Nonetheless, undetected values followed a coherent pattern with nearly all of them being distributed in expected low hormonal profiles. Frequencies of undetected E2 were distributed as follows: 57 COC users, 27 NC women in the early follicular phase, and 23 men (107/110, 97.27%). Given that the luteal phase of the menstrual cycle was not targeted in NC women, all our participants were recruited in a low state of P. This provides insight into the high undetected P rate in every group (≥ 64.5%). For T, 39 COC users, 11 early follicular NC women, and 1 pre-ovulatory NC woman accounted for the undetected values. For COC users in the active phase of their regimen (n = 55), EE levels were below the LOD for 8 women (14.5%). Given that our missing values were not missing at random (NMAR) and as recommended by Herbers et al. (132), we addressed this issue with multiple imputation from a fitted lognormal distribution using their same procedure implemented in R with the fitdistrplus and EnvStats packages (133, 134; the code used by Hebers et al. is accessible at https://osf.io/spgtv/). This imputation method was found to be valid, despite the high proportion of missingness in left-censored datasets (135). Using the SPSS toolbox, we generated 10 datasets with simulated values between 0.01 pg/mL and the LLOQ for each endogenous hormone (e.g., 0.01-0.29 pg/mL for E2). The LOD was used as the upper bound for EE. For statistics, we reported the mean value of each metric (e.g., F_M = average of the 10 Fs). Due to highly unbalanced distributions and to be coherent with the imputation procedure performed, we then log-transformed endocrine variables. Finally, outliers were examined. Z-scores were derived for men, early follicular NC women, pre-ovulatory NC women, and COC users on log-transformed variables. Using a criterion of ≥ ± 3.29 (136, 137), we identified two outliers for E2 (COC users), one outlier for P (COC user), two outliers for T (pre-ovulatory NC woman and man), and one outlier for EE. All outliers were winsorized with respect to groups, using the next highest value of each group (138) on the imputed log-transformed variables. Additionally, we removed endocrine values for one participant (man) due to an invalid assessment (technical issue/low internal standard). For descriptive purposes only, we winsorized imputed raw (untransformed) values and reported endocrine concentrations in pg/mL.

All ROI analyses were performed using SPSS 27 (IBM). We conducted general linear models (GLMs) for each ROI with Group (men, never users, past users, current users) as the between-subjects factor and Hemisphere (left, right) as the within-subjects factor. We controlled for age for both GMV and CT analyses, as well as for TIV in GMV analyses (for a justification, see 139, 140). Given our between-group design, we examined further potential confounders thereafter. Covariates were selected from analyses performed on various sociodemographic and psychological outcomes with a threshold set at p <.10. First, study groups were compared using ANOVAs and chi-square tests for continuous and categorical variables, respectively. Second, correlations were performed between potential confounders and TIV/global CT. Using the ‘disjunctive cause criterion’ approach for covariate selection (141), we selected variables that were either associated with the study groups, brain morphology (TIV/global CT), or both. Thus, we decided to initially run the minimally adjusted models (e.g., age and TIV for GMV) and re-run the fully adjusted models for additional confounders (e.g., ethnicity, relationship status, religiosity). To investigate the interplay between current dynamics of the menstrual cycle and history of COC use, we focused on both groups of NC women (namely, never users and past users, as well as the cycle phase they were in during data collection). Accordingly, we performed Group (past users, never users) x Cycle phase (early follicular, pre-ovulatory) GLMs. Similar to the previous model, we implemented the same covariate selection plan.

To investigate the impact of endogenous sex hormones, we re-ran previous Group x Hemisphere GLMs for each ROI with added hormonal concentrations. Given that these specific analyses relied on the continuous nature of hormonal data, we excluded variables showing a proportion of imputation higher than 60%. This procedure led to the exclusion of COC users in this set of analyses and P levels for the three remaining groups (i.e., men, never users, past users). Conceptually, this decision is supported by the fact that COC users should be in an endogenous ‘shut-down’ state and that P variability was expected to be minimal across participants (due to the nature of our study design). As such, models were built to explore the effects of Group (3) x Hemisphere (2) x E2 and Group (3) x Hemisphere (2) x T, all while controlling for age (and TIV for GMV). In a second step, significant results were re-examined by running a sensitivity analysis with subcategorized hormonal data (e.g., ‘low levels’ = undetected imputed values, ‘medium levels’ = detected but low values, ‘high’ = detected and higher values). This was performed to assess the robustness of the results obtained with multiple imputations. No WBA was performed for this analysis subset given that multiple imputation could not be implemented in FreeSurfer and SPM.

In COC users, we explored the role of exogenous hormones as assessed by salivary EE levels or EE prescribed dose (contained in each pill, in μg) and androgenicity. As EE levels were accounted for, EE models were restricted to women tested in the active phase of their COC (n = 55). Among them, 48 used monophasic regimens, while seven used triphasic regimens with a constant dose of EE across all active pills. Given the few and unevenly distributed data points for EE dose (as reported in Table 1 ), we grouped users into low (10-25μg) or high (30-35μg) dosages. In our sample, no COC users used 50μg of EE. For each ROI, EE GLMs were tested with Hemisphere (left, right) x EE dose (low, high) and Hemisphere x EE levels as variables of interest, where age and TIV (for GMV only) were used as covariates. Of note, we did not investigate the impact of lifetime EE dose in current users or previous EE dose within past users given that these calculations led to insufficient sample sizes (i.e., 33/62 current users and 14/37 past users having used the same dosage over lifetime use). As for progestin androgenicity, we classified users into high (n = 27), low (n = 18), and anti-androgenic (n = 17) COCs categories (24) (see Table 1 ). The influence of androgenicity was explored with Hemisphere (right, left) x Androgenicity (anti, low, high) GLMs, while controlling for age and TIV (for GMV only). Considering the small sample size, we first ran separate models for EE and progestin androgenicity. However, EE dosage and progestin androgenicity can vary within formulations and potentially interact with each other. Indeed, the chi-square test showed that a higher proportion of low EE doses were classified as high androgenic COCs, while high EE doses were more likely to be classified as low androgenic COCs (p = .004). Therefore, in a second step, we adjusted for one or the other in their respective set of analyses. Although statistically underpowered, this allowed us to verify the stability of our results when covarying out an important confounder.

For ROI-based analyses, significance levels were set at p <.05. Trend-level results were defined as.05 ≤ p ≤.08. Bonferroni corrections were applied to p-values in post hoc analyses. Based on the novelty of the study, uncorrected results were presented for multiple comparisons; however, adjustment for multiple comparisons was applied for the seven GMV ROIs and four CT ROIs using the False Discovery Rate (FDR) correction (142). For exploratory volumetric and thickness WBA in FreeSurfer, the Monte Carlo Null-Z Simulation was carried out to control for multiple comparisons (10 000 iterations, cluster-forming two-tailed p <.05, cluster-wise p <.05). A Bonferroni correction was applied as we tested both hemispheres. For VBM, we used a primary uncorrected threshold of p <.001 with a cluster extent of k > 10 voxels and secondary cluster-level FDR-corrected threshold of q <.05. Groups were compared using a F-test orthogonal contrast. Effect sizes were reported as partial eta squared (η_p ² ). Error bars in figures denote ± standard error of the mean.

To account for potential diurnal variations in sex hormone levels (143–145), we examined whether the timing of study sessions, which ranged from 8:15AM to 7PM as a matter of feasibility, influenced endocrine concentrations. Our analysis revealed no significant influence of session timing (AM/PM) on hormone levels within each group (p ≥.130). Groups were also equally distributed for AM/PM testing (p = .190). To confirm hormonal status, GLMs were conducted on groups based on their current hormonal state. Levels of E2, P, and T were compared between men (n = 40), early follicular NC women (n = 45), pre-ovulatory NC women (n = 32), and COC users (n = 62). Analyses revealed a group difference for E2 [F_{M (3, 175)} = 34.80, p_M <.001, η_p ² _M = .369] and T [F_{M (3, 175)} = 213.98, p_M <.001, η_p ² _M = .784]. Games-Howell post hoc tests showed higher E2 concentrations in pre-ovulatory NC women compared to the three other groups (p_M s <.001). Early follicular NC women and men had similar E2 levels (p_M = .898), which were both comparable to COC users (p_M = .267 and.117 respectively). For T, men showed higher concentrations than the three groups of women (p_M s <.001). Pre-ovulatory NC women exhibited greater T levels than early follicular NC women (p_M = .033) and both NC groups had more T levels than COC users (p_M s ≤.006). Further, a Group (2) x Cycle phase (2) GLM was conducted between past users and never users. For the three endogenous hormones, no differences were found between the two groups nor for the interaction term (p_M s ≥.185). For cycle phase, irrespective of groups, E2 and T were significant and followed the same pattern as in the first set of hormonal analyses. Untransformed hormonal means are reported in Table 1 , while both sets of group analyses are summarized in Figure 1 . In sum, the salivary endocrine assessment revealed expected variations in our study groups. Never and past users were equivalent in terms of their hormonal profile.

Sample characteristics are presented in Table 1 . As expected, men had greater TIV than the three other groups of women (p <.001). No differences were observed for global CT. Regarding sociodemographic variables, men, never users, past users, and current COC users did not differ on age, mother tongue, body mass index, and physical activity practice (ps ≥.101). However, relationship status was disproportionately distributed across study groups (p <.001). Based on Bonferroni-adjusted pairwise Z-tests, current users were significantly more in relationships than the three other groups and the proportion of individuals in a relationship was higher in both current and past users than never users and men. Moreover, groups tended to differ on ethnicity (p = .059), where never users were less of Caucasian origin than expected. Although non-significant, a greater proportion of past users were of Caucasian origin (78.4%) than never users (52.5%) in our sample. Years of education were also marginally different between the groups (p = .072), although post hocs did not yield significant group differences (ps ≥.156). Concerning psychological constructs, groups did not differ regarding religiosity, neuroticism, trait anxiety, depressive symptoms, and exposure to adverse events (ps ≥.151). For endocrine-related variables ( Table 2 ), no differences were found for the age at menarche between the three groups of women (p = .498), nor for cycle phase between past and never users (p = .524), and reasons for COC use between current and past users (p = .329).

Correlations between TIV/global CT and potential confounders across the whole sample are reported in Table 3 . According to the disjunctive cause criterion (141), fully adjusted models for GMV included age (conceptually), TIV, ethnicity, years of education, relationship status, and religiosity, whereas age, ethnicity, years of education, relationship status, and practicing physical activity were controlled for in fully adjusted CT models. When restricted to past users and never users, fully adjusted models included age and TIV (conceptually), mother tongue, body mass index, ethnicity, relationship status, religiosity, and adverse events for GMV analyses and age, mother tongue, body mass index, ethnicity, relationship status, and neuroticism for CT analyses (data not shown).