Sex-biased gene expression during neural differentiation of human embryonic stem cells

Sex differences in the developing human brain are primarily attributed to hormonal influence. Recently however, genetic differences and their impact on the developing nervous system have attracted increased attention. To understand genetically driven sexual dimorphisms in neurodevelopment, we investigated genome-wide gene expression in an in vitro differentiation model of male and female human embryonic stem cell lines (hESC), independent of the effects of human sex hormones. Four male and four female-derived hESC lines were differentiated into a population of mixed neurons over 37 days. Differential gene expression and gene set enrichment analyses were conducted on bulk RNA sequencing data. While similar differentiation tendencies in all cell lines demonstrated the robustness and reproducibility of our differentiation protocol, we found sex-biased gene expression already in undifferentiated ESCs at day 0, but most profoundly after 37 days of differentiation. Male and female cell lines exhibited sex-biased expression of genes involved in neurodevelopment, suggesting that sex influences the differentiation trajectory. Interestingly, the highest contribution to sex differences was found to arise from the male transcriptome, involving both Y chromosome and autosomal genes. We propose 13 sex-biased candidate genes (10 upregulated in male cell lines and 3 in female lines) that are likely to affect neuronal development. Additionally, we confirmed gene dosage compensation of X/Y homologs escaping X chromosome inactivation through their Y homologs and identified a significant overexpression of the Y-linked demethylase UTY and KDM5D in male hESC during neuron development, confirming previous results in neural stem cells. Our results suggest that genetic sex differences affect neuronal differentiation trajectories, which could ultimately contribute to sex biases during human brain development.


Introduction
The importance of studying sex differences is today well recognized in many research fields, including neuroscience (Jazin and Cahill, 2010;Zagni et al., 2016;Golden and Voskuhl, 2017;Choleris et al., 2018).However, no studies have yet specifically addressed sex differences in multiple human embryonic stem cell lines (hESC).Models of stem cell differentiation are invaluable for this kind of research because they allow the study of genetic sex differences in the absence of sex hormones and their effects.Sex differences can have critical implications for the tissues they affect.In the developing nervous system, they have the potential to influence, e.g., the number of neurons, glial cells, synapses, and their organization.As a consequence, neural networks may develop that respond to stimuli in sex-specific ways (Arnold and McCarthy, 2016).Sex differences appear to be particularly relevant for several neurological diseases.Autism spectrum disorder, Tourette syndrome, and attention-deficit/hyperactivity disorder, for example, exhibit sex biases in symptoms, onset, and prevalence (Boyle et al., 2011;Loke et al., 2015;Tesic et al., 2019).These disorders are believed to have a neurodevelopmental origin, and it is likely that sex differences in early brain development act as a major contributor to the sex-biased development and manifestation of these diseases.Currently, sex differences in the human brain are largely attributed to the effects of sex hormones during critical periods of fetal development and adolescence (McCarthy, 2009;Arnold and McCarthy, 2016).However, an increasing body of evidence demonstrates that not only hormones alone play a role but also genetic factors, such as sex chromosome dosage compensation (Raznahan et al., 2018;Sidorenko et al., 2019;Warling et al., 2021), and differences in gene expression contribute to sex differences in the developing human brain (Carruth et al., 2002;Corre et al., 2016;Arnold, 2017;Gegenhuber and Tollkuhn, 2020;Cabrera Zapata et al., 2022).Despite the growing attention that research in sex differences receives, only a few studies specifically aim to investigate sex differences in embryonic stem cells, and to our knowledge, none of them include more than 1 cell line of each sex.Particularly, the role of sex-biased gene expression in early neuron development remains unclear.To address this gap, we collected four male and four female human embryonic stem cell lines, differentiated them into a mixed neuronal population over 37 days, and analyzed their gene expression profiles for sex differences using both total RNA sequencing and qPCR analysis.Our results reveal the presence of gene expression sex differences in all investigated stages of neuronal differentiation, but most pronounced after 37 days.We identified a number of candidate genes with potential contributions to sexbiased neuronal development.Finally, we demonstrate how overexpression of X-linked gametologous genes that escape X chromosome inactivation (XCI) is compensated by the expression of corresponding Y-linked homologs.
In summary, our research addresses a critical gap in understanding how sex differences in gene expression may influence early neuronal development.These findings could have implications for understanding the basis of sex-specific responses in the nervous system and their relation to neurological disorders.

RNA sampling and sequencing
RNA sampling for qPCR as well as for RNA sequencing was performed in the following way: cells were washed with 2 mL PBS, thereafter 1 mL TriZol (ThermoFischer 15596026) was used for lysation.The lysate was passed six times through a syringe (BD Microlance 3 (27G ¾ (Nr.20) 0.4 × 19 mm) Ref: 302200) to promote cell lysis.The TriZol phases were separated using chloroform and the RNA phase was purified with the PureLink RNA Mini Kit (ThermoFischer 12183018A).Subsequently, the RNA was depleted of DNA using the DNA-free DNA Removal Kit (ThermoFischer AM 1906).All steps were executed according to manufacturer's instructions.
Cell lysates were obtained from one 6-well-plate well before differentiation (D0) and after start of differentiation (D4, D9, D17, D27, D37) using 1 mL TriZol (ThermoFischer 15596026).TriZol phases were separated using chloroform and the RNA phase was purified with the PureLink RNA Mini Kit (ThermoFischer 12183018A) according to manufacturer's instructions.Subsequently, the RNA was depleted of DNA using the DNAfree DNA Removal Kit (ThermoFischer AM 1906).
RNA samples were sent in triplicates for sequencing to the SNP&SEQ Technology Platform in Uppsala where the RNA quality was verified using a Bioanalyzer (RIN ≥8).Samples were then treated with the TruSeq Stranded Total RNA and Ribo-zero Gold kit (Illumina).The treatment preserves strandedness of transcripts but depletes the samples of cytoplasmatic and mitochondrial ribosomal RNA.Sequencing was performed on the NovaSeq instrument with an S4 lane (Illumina) resulting in approximately 27.7 million 150 bp pair-read read-pairs per sample.The data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus (Edgar et al., 2002) and are accessible through GEO Series accession number GSE249035 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE249035).

Differential gene expression and gene set enrichment analysis
The samples were aligned to GRCh38_r96 (Schneider et al., 2017) with Hisat2 (Kim et al., 2019) and sequencing reads were assigned to genomic features using featureCounts (Liao et al., 2014).Quality control was performed with FastQC and MultiQC (Ewels et al., 2016).The differential gene expression analysis was performed with DEseq2 (Love et al., 2014).The two factors "Sex" and "Time" were joined and the contrast function was used to compare the differential gene expression between and within the biological sexes (XX to XY) at each timepoint (D0, D4, D9, D37), Figure 3 A-D.Differentially expressed genes were then filtered for a minimal average of 50 counts in any of the cell lines of one sex, and for a fold-change bigger than 0.5 with an adjusted p-value lower than 0.05.The count limit of 50 accounts to approximately 2 CPM in our data and is considered stringent.The value was selected after inspection of the histogram of average logCPM values and an appropriate threshold was selected between the low-CPM peak representing non-expressed genes and the high-CPM peak representing expressed genes.
To acquire genes that are truly sex different, a stringent filter was implemented.It consisted of the standard parameters: fold-change (≥0.5) and adjusted p-value (≤0.05), but also includes the coefficient of variation for genes within cell lines of each sex (≤0.5) and a minimal difference between the expression of cell lines of the opposite sex (≥20% difference in gene expression).The filter kept only those genes that are expressed at a similar level within cell lines of the same biological sex and if none of the cell lines of one sex expresses the gene of question as high as any cell line of the other sex.For this, the coefficient of variation for genes within the cell lines of each sex was required to be lower than 0.5, to remove genes that are highly variable in expression and hard to predict.Additionally, the difference in gene expression between the cell lines of the opposite sex had to be higher than 20%.The latter required all 9 comparisons (3 male vs. 3 female lines) to display a difference of at least 20% in DEseq2 normalized counts.To investigate genes that contribute significantly to the neuronal differentiation, XX and XY counts were combined and the timepoints were contrasted (D0 to D4, D0 to D9 and D0 to D37), Figures 3E, F.
For the gene set enrichment analysis, the GSEA/MSigDB software was used, software version 4.2.3 and molecular signatures database version 7.5.1 (Mootha et al., 2003;Subramanian et al., 2005).Preranked GSEAs were performed using the stat parameter of the DEseq2 analysis, which represent the difference between the selected comparison including the magnitude of fold-change and p-adjusted.The C5 ontology gene sets were used consisting of gene sets for biological processes (BP), cellular components (CC) and molecular functions (MF).The resulting data was used to create enrichment maps using the Cytoscape software version 3.91 (Shannon et al., 2003) and the auto annotate app for cluster annotation based on gene set description.Enrichment maps were created using an FDR q-value cut-off of below 0.1 and a perfused force directed layout based on the similarity coefficient of gene sets.

qPCR
For the gene expression analysis by qPCR, cDNA was synthesized from the purified and DNA free RNA using Maxima H Minus Reverse Transcriptase (ThermoFischer EP0752) with anchored Oligo (dT)20 Primer (ThermoFischer 12577011).The following protocol was used for cDNA synthesis in a 15 µL reaction: 500 ng rDNAse-treated RNA, 1.25 μg/μL Oligo (dT) 20 Primer, 1 µL dNTP mix (10 mM each, Thermo Fisher R0191), 4 µL RT buffer, 0.5 µL RiboLock RNAse inhibitor (ThermoFisher EO0381), 0.5 µL Maxima RT enzyme and RNAse free water was used to top up to the desired volume.Incubate at 50 °C for 30 min, terminate at 85 °C for 5 min.A NanoDrop spectrophotometer was used to assess the RNA amount and quality.RNA integrity was measured using gel electrophoresis and comparison of 28S/18S ribosomal RNA ratios.The qPCR analysis was performed on the Applied Biosystems 7500 Real-Time PCR System using PowerUp SYBR Master Mix (Applied Biosystem A25742) and standard cycling mode.For every reaction 1 ng of cDNA was used with 7.5 μL of SYBR master mix and 1 μL of 10 mM forward and reverse primers.Molecular grade H 2 O was added for a total volume of 15 μL.qPCR measurements were taken in triplicates.The thermal cycling profile was selected according to master mix recommendations (1 cycle UDG activation, 50 °C, 2 min; 1 cycle activation of Dual-Lock DNA polymerase, 95 °C, 2 min; 40 cycles of denaturation, 95 °C, 15 s, and annealing/extension, 62 °C, 1 min), all runs were performed at an annealing temperature of 62 °C.Amplification cycles were repeated 40 times followed by a melting stage, serving as a diagnostic tool to assess amplicon length.Primer sequences can be found in Supplementary Table S9.

Similar neuronal differentiation of male and female human embryonic stem cells
To identify potential sex differences in the development of human embryonic stem cells (hESC), we differentiated four male and four female-derived hESC lines towards a population of mixed neurons.The differentiation protocol utilizes adherent cell culture combined with small molecules for neural induction and maturation, as described in detail in the Methods section.The cells were differentiated for a total of 37 days, and samples were collected at six timepoints: D0, D4, D9, D17, D27, and D37 (Figure 1A).Six cell lines were analyzed by total RNA sequencing at four timepoints (D0, D4, D9, D37), and all 8 cell lines were analyzed via qPCR at all six timepoints (Table 1).
A principal component analysis of expressed genes revealed that male and female cell lines show similar gene expression at D0 but display cell line-specific differences at D4 and D9 (Figure 1B).After 37 days of differentiation, the differences decreased, and male as well as female cell lines show similar gene expression.However, we noticed a clear sex-specific separation between them according to principal component 2 (PC2).
In general, the cell lines of both sexes have differentiated similarly towards neurons, displayed by very similar expression dynamics of marker genes for pluripotency (OCT4, NANOG), Similarities between male and female hESC in neuronal differentiation: (A) Overview of the experimental procedure: Four hESC lines of each sex were differentiated for 37 days to a population of mixed neurons induced by small molecules in 2D culture.Bulk RNA sequencing and qPCR.(B) A PCA analysis of expressed genes displays similar gene expression patterns at D0, some inter-cell line differences during early differentiation, and a sex-biased separation at D37 (XY: WA14, HS401, HS983a, XX: WA09, HS975, HS980).(C) Marker gene expression for pluripotency (OCT4, NANOG), ectoderm (PAX6, SOX1), immature neurons (DCX, NCAM), mature neurons (MAP2, DLG4, SYP), post-mitotic neurons (ELAVL2, TUBA1A, NDN), as well as for proliferation (MKI67, PCNA), and apoptosis (CASP3, BCL2) show a similar expression in male and female cell lines.(D) Markers for specialized neurons (GABAergic, Dopaminergic, Glutamatergic) display a premature state of the neurons after 37 days of differentiation.(E) A gene set enrichment analysis at D4, D9, and D37 reveals an upregulation of biological processes implicated in neuronal differentiation in all cell lines.The GSEA was summarized in an enrichment map.Each dot represents a gene set; the size of the dot corresponds to the size of the gene set.If the dots are connected or in close proximity, the more the gene sets are overlapping, and the closer related are their associated biological functions.Gene sets were grouped by a common GO-term classification and named after an overarching biological function.Biological functions implicated in neuronal development are highlighted in blue.Gene sets are filtered by FDR q-values <0.1.Figure 1A created with BioRender.com.
ectoderm (PAX6, SOX1), immature neurons (DCX, NCAM1), mature neurons (MAP2, DLG4, SYP), post-mitotic neurons (ELAVL2, TUBA1A, NDN), as well as for proliferation (MKI67, PCNA) and apoptosis (CASP3, BCL2), Figure 1C.In addition, the pluripotency state of all cell lines was investigated and identified as primed by the overexpression of naive over primed pluripotency marker genes at T0, Supplementary Figure S1.To estimate the replicability of the differentiation protocol, the experiment was repeated two more times, and the expression of NES, PAX6, and DCX was measured using qPCR.We observed consistent gene expression levels in each repetition, indicating low variability and a robust differentiation protocol.Additionally, we identified sex differences which were not detected in the RNAseq data for NES at D9 (1.8 FC, p < 0.0001), in PAX6 at D37 (1.3 FC, p < 0.05), and in DCX at D17 (2.2 FC, p < 0.0001), Supplementary Figure S2.
An investigation of marker genes for specialized neurons has revealed that male as well as female cells similarly differentiated to premature neuron-like cells after 37 days (Figure 1D).For example, the premature state of GABAergic neurons is suggested by the high expression of GABBR1 (GABA type B receptor unit) combined with high expression levels of MEIS2 and ABAT (both GABAergic precursor markers), while other GABAergic markers, such as the type A receptor GABRA1/2 or the neurotransmitter transporters SLC6A1 (GAT) and SLC32A1 (VGAT), as well as the rate-limiting enzyme for GABA synthesis GAD1 are expressed very low or not at all.Similarly, dopaminergic markers NR4A2 and LMX1A, which are transcriptional regulators for the development of DA neurons, indicate the presence of neurons that are not fully mature.However, the expression of the catecholamine precursor L-DOPA-synthesizing-enzyme tyrosine hydroxylase (TH) and the L-DOPA to dopamine converting enzyme dopa decarboxylase (DDC), as well as the dopamine receptor D2 (DRD2), point to an, at least, early post-mitotic state of some dopaminergic neurons.The dopamine transporters SLC6A3 (DAT1) and SLC18A2 (VMAT2), as well as the transcription factor FOXA2 that supports maintenance of mature DA neurons, are expressed at low levels or not at all.As for the glutamatergic neurons, we can see high expression of the glutamate transporters GLT1 and VGLUT2, as well as for the glutamate AMPA type receptors subunits GRIA1 and GRIA2.Glutamate NMDA type receptors were also expressed, albeit at a moderate to low level (data not shown).Interestingly, the VGLUT2 gene shows a sex difference and is higher expressed in male cells than in female cells after 37 days of differentiation (p = 0.009 **).
Gene set enrichment analyses (GSEA) also revealed a similar differentiation of male and female cell lines.We performed GSEAs for each sex separately (e.g., XY D0 to D4 and XX D0 to D4) and compared them.We found that the upregulated gene sets representing biological processes in the cells during the differentiation experiment are almost identical with only minor variations between the sexes (data not shown).A simple way to visualize this is by performing a GSEA with pooled data from all cell lines.This has resulted in a single enrichment map in which gene sets will only appear if they are upregulated in all cell lines, male and female (Figure 1E).Each dot in the map represents an upregulated gene set.Upregulated processes involved in neuronal development are, e.g., neural tube development at D4, forebrain development at D9, and neurogenesis, nervous system, and neuron projection development at D37.The processes involved in nervous system development are highlighted in blue.A full list of the enriched gene sets can be found in Supplementary Table S1.

Differences in XCI erosion in female cell lines displays little effect on autosomal and X-linked gene expression
Since X chromosome inactivation (XCI) erosion can affect the gene expression of both X-linked and autosomal genes, we investigated the XCI erosion status in our female cell lines.A key indicator for the presence of XCI erosion is the absence of XIST expression.In our study, one of the three primed female cell lines used for RNAseq analysis expresses XIST at a relatively low level (HS980 with 1.5 TPM at D0) without significant changes during differentiation.The other 2 cell lines (WA09/LT2e and HS975) do not express XIST at any time point.To assess the extent of XCI erosion in the female cell lines, we calculated the expression ratio of X-linked genes between male and female cell lines.Since XCI erosion leads to an increased expression of X-linked genes, we can use the XX:XY expression ratio to estimate the degree of XCI erosion.In line with Bar et al., 2019, we consider a ratio of >1.1 as a state in which X-linked genes in female cell lines are increasingly expressed due to XCI erosion.XCI escape genes display an average XX:XY ratio of 1.29 (data not shown) and were removed before the analysis.
On average, the majority of X-linked genes (67%) in the three female lines demonstrate an expression ratio indicative of one X chromosome being inactivated (XaXi), while the remaining 33% of X-linked genes reflect either a state of XCI erosion (XaXe) or two active X chromosomes (XaXa), Supplementary Figure S3A.When analyzed individually, the cell lines exhibit different states of XCI erosion.The lines HS975 and WA09/LT2e display a higher degree of XCI erosion with 31% and 48% of X-linked genes, respectively, with an XX:XY ratio above 1.1 and thus within the XaXe or XaXa range.While in the line HS980 only 16% of X-linked genes display a ratio above 1.1 (Supplementary Figure S3B).To estimate the number of affected genes, we used the average XX:XY gene expression ratio of chromosome 5 and 16 as a baseline, given their similarity in the number of protein-coding genes compared to chromosome X.After subtracting the baseline, the final fraction of X-linked genes upregulated due to the potential effect of XCI erosion was estimated to be on average 6%, corresponding to 38 out of 634 genes (4% or 25 genes for HS980, 4.5% or 29 genes for HS975% and 10% or 63 genes for WA09/LT2e).The complete list of X-linked genes and their corresponding XX:XY expression ratios can be found in Supplementary Table S2.Based on this, only two of the gametologous genes (CASK and SHROOM2) were identified as possibly affected by XCI erosion since they exhibit an XX:XY ratio above 1.1 and are not XCI escape genes.Furthermore, no significant difference was found between the average X-linked gene expression of the cell lines using a oneway ANOVA with Tukey's multiple comparison test, (Supplementary Figure S3C).Moreover, the variance of the average gene expression from the X-chromosome among the male cell lines (F = 1.33) was twice as high than among the female cell lines (F = 0.65).This suggests that the variation in XIST expression does not critically impact X chromosome gene expression levels in the female lines, rendering them comparable in their X-linked gene expression despite their difference in XIST expression and XCI erosion.
To evaluate whether XIST expression and XCI erosion status affect the autosomal gene expression in the female cell lines, we investigated the relation between X chromosomal and autosomal gene expression.For this, we used the TPM gene expression values to calculate the X:A (Supplementary Figure S3D).We did not observe significant difference in autosomal gene expression in the Xist-positive cell line HS980.A slight trend for increased expression can be noted on chromosomes 8, 9, 11, and 22 with their ratios increased by 0.15-0.2.However, the average X:A ratio of all autosomes was only slightly increased by 0.09 points in HS980, compared to the average of the Xist-negative lines (HS980: 1.16, LT2e: 1.11, HS975: 1.03).

X and Y-linked genes contribute to neuronal differentiation
To identify how many Y-linked, X-linked, and autosomal genes contribute to neuronal differentiation, we investigated their gene expression levels and changes during the differentiation experiment (Figures 2A-C; Supplementary Table S3).We used the average gene expression value of all cell lines grouped by sex and implemented a rigorous threshold by excluding genes with a count number below 50.This stringent criterion was applied to effectively remove a large number of lowly expressed genes.Subsequently, we categorized resulting genes based on the timepoint at which they are expressed highest.We noticed a total of 15,222 autosomal genes, 568 X-linked genes and 23 Y-linked genes to be expressed in the embryonic stem cell lines (Supplementary Table S3).For autosomal and X-linked genes, a similar proportion of genes were upregulated at certain timepoints during differentiation.The largest share of genes was upregulated at D37, accounting for up to 34% of autosomal genes and 39% of the X-linked genes.This is followed by genes that were downregulated during differentiation (26% autosomal, 28% X-linked) and, again in a similar proportion, genes with unchanged expression during differentiation (25% autosomal, 26% X-linked).Only a small percentage of X-linked and autosomal genes seemed to be upregulated at the stages of neuronal induction at D4 (4%-5%) and at D9 (6%-7%).As for the Y-linked genes, the largest number of genes was upregulated early during neuronal induction (35%, 8 genes) at D4.Only one gene is upregulated at D9 in the late stage of neuronal induction (4%, 1 gene).
An equal number of Y-linked genes are upregulated at D37 (21%, 5 genes) or decreased during neuronal differentiation (26%, 6 genes), while 14% (3 genes) of Y-linked genes do not change their expression during differentiation.We have labeled the four genes with the highest expression during differentiation at their corresponding timepoints in the panels D4, D9, and D37 in Figure 2. The Y-linked genes with the highest expression were TXLNGY, NLGN4Y, and ZFY at D4, RPS4Y1 at D9, and USP9Y, UTY, KDM5D, and PCDH11Y at D37 (Figure 2A).The X-linked genes with the highest expression were ZIC3, MBNL3, and KLHL4 at D4, EFNB1, TMSB15A, and FHL1 at D9, and DCX, MAGED2, GPM6B, and MAGED1 at D37 (Figure 2B).The autosomal genes with the highest expression were BCAT1, PDOXL, GJA1, and FAT1 at D4, RN7SL1, RN7SL2, SNORD3A, and ENO1 at D9, and MALAT1, TUBA1A, ACTG1, and TXNIP at D37 (Figure 2C).The expression data of these and additional genes, as well as genes that significantly decrease during differentiation or remain stable in expression, can be found in Supplementary Table S3 and Supplementary Figure S4.

Identifying genes with sex-biased expression and impact on neuronal differentiation
To identify sex-biased gene expression during neuronal differentiation, we performed a differential gene expression analysis with the total RNA sequencing (RNAseq) data.In the assay, we compared 3 male and 3 female cell lines at 4 time points during differentiation, this resulted in the following four comparisons XX to XY at D0, D4, D9, and D37.After applying a filter to exclude genes with low expression (<50 counts), and an adj.p-value >0.05, we obtained 1,152, 2,330, 1,497, and 2,915 differentially expressed genes (DEG) for D0, D4, D9, and D37, respectively (Table 2).However, many of those genes did not meet our requirements for sex-based differential expression.Some DEG displayed a too high variation in expression between the cell lines of one sex to be a reliable candidate for a gene with a sex-biased expression.In other cases, individual cell lines of one sex showed a too similar expression level to cell lines of the opposite sex, so that also these candidates were not a good match for a gene with a sexbias in expression.To ensure that we select true sex differences, we applied another set of stringent post hoc filtering rules that only assign genes as sex-biased if they are expressed at a similar level within cell lines of the same biological sex (coefficient of variation <0.5), and that required all cell lines of one sex to express a gene at least 20% higher (based on DEseq2 normalized counts) than all cell lines of the other sex.Additional details can be found in the methods section.After filtering, the numbers of differentially expressed genes have reduced substantially to 65, 218, 114, and 148 for D0, D4, D9, and D37, respectively (Table 2).The complete list of DEG can be found in Supplementary Table S4.
As expected, genes of the sex chromosomes (X and Y) make up a considerable part of the sex-biased DEG at all time points (Table 2).The Y chromosome genes display the highest fold-changes among the DEG in male cells.Similarly, the X chromosome genes make up a large part of the overexpressed genes in female cells.Not surprisingly, the X-linked genes that escape X-inactivation (XCIescapees) show the highest fold-change among the DEG in female cells.Apart from the genes on the sex chromosomes, there is also a significant share of DEG expressed from the autosomes, 43% at D0, 78% at D4, 73% at D9, and 76% at D37.Interestingly, the number of sex-biased DEG increases with the start of differentiation at D4 (Figures 3A-D).While at D0 there are 65 sex-biased genes expressed (27 XY,38 XX), at D4 it is 218 genes (137 XY, 81 XX), at D9 it is 114 genes (28 XY, 86 XX), and at D37 it is 148 genes (104 XY, 44 XX).Next, we investigated the sex differences at D37 in more detail.Most of the sex-biased genes at D37 are higher expressed in the male cell lines, 104 compared to 44 genes in the female lines (Figure 3D; Supplementary Table S4).From the upregulated genes in the female lines, 21 are X-linked genes, and of these, 9 are known to escape XCI, the remaining 22 genes are autosomal.In contrast, the largest number of upregulated genes in males are expressed from the autosomes (89 genes).Only a small number of the upregulated genes at D37 in the male lines are encoded on the Y chromosome (15 genes).We have annotated all sex-biased genes at D37 for their appearance in neuronal, synaptic, as well as developmental gene sets using the gene ontology database (GO) and collected evidence from available scientific literature regarding their implication in neuronal processes (Supplementary Table S5).Of the 104 DEG in male cell lines, 32 genes were implicated in neuronal processes (Top 5: NLGN4Y, UTY, PCDH11Y, NHLH2, EBF1), while of the 43 DEG Gene expression of sex chromosome and autosomal genes during differentiation: Gene expression of (A) Y-linked, (B) X-linked, and (C) autosomal gene expression during neural differentiation.The pie-chart shows the number of genes peaking in expression at different timepoints (D4, D9, or D37), as well as the number of genes decreasing during differentiation and genes that maintained a stable expression level.Additional panels display the expression dynamics of genes peaking at D4, D9, and D37.The genes with the highest expression levels at the corresponding timepoint are labeled.According to their expression peaks, Y-linked genes play a role at different timepoints during neural differentiation in male stem cell lines.*XCI-escaping genes selected according to Katsir and Linali (2019), Garieri et al. (2018), Zhang et al. (2013).
in female cell lines it was 7 genes (Top 5: SYAP1, TACR3, GRM1, NDNF, AMOT).Genes that change their expression pattern during neuronal differentiation and do that in a sex-dependent manner are likely to be important candidates that cause sex differences in brain development.To identify these genes, we looked for sex-biased genes among the highly up-and downregulated genes during differentiation (D0 to D4, D0 to D9, and D0 to D37) and found between 2,218 and 5,518 genes that are up-or downregulated, respectively (Figures 3E-G).Approximately 1% of the genes show a sex-bias in gene expression (Supplementary Table S6).
We have added this sex-biased expression during neural differentiation as a parameter in Supplementary Table S5 and used it for the selection of candidate genes later.

Sex-biased gene expression in neural processes after 37 days of differentiation
To identify potential differences in biological processes between cell lines of different sexes, we conducted additional GSEAs comparing male and female gene expression at the timepoints D0, D4, D9, and D37.Enrichment maps were generated from the GSEA results to identify clusters of upregulated gene sets with similar functions.Clusters of gene sets are more likely to have a biological relevance than isolated gene sets.Using this approach, identified various biological processes before differentiation and at the early (D4) and late stages of differentiation (D37).Interestingly, at D9, the analysis did not reveal any differences between the cell Sex-specific differences at key timepoints and in genes that change significantly during differentiation: (A-D) Sexually differentially expressed genes comparing male (XY) against female (XX) cell lines at different time points of neuronal differentiation.Genes overexpressed in male cell lines are marked in blue; green if they are located on the Y chromosome.Similarly, genes overexpressed in female cell lines are marked in red; purple if they are expressed from the X chromosome and pink if they escape XCI.The number of sex-biased genes overexpressed in male and female cell lines is as follows: (A) D0, 65 genes in total, 25 in XY and 38 in XX; (B) 218 genes in total, 137 in XY and 81 in XX; (C) 114 genes in total, 28 in XY and 86 in XX; (D) 148 genes in total, 104 in XY and 44 in XX.The number of genes with sex-biased expression from sex chromosomes remains relatively constant with an average of 15 Y chromosome and 22 X chromosome genes expressed at each time point.XCI-escapee genes are among the sex-biased genes with the highest p-value in the XX cell lines.(E, F) display the DEG that are up-and downregulated over the course of neuronal differentiation.Genes that are additionally expressed in a sex-biased manner are highlighted in red and blue according to their sex-bias in female and male cell lines, respectively.On average, 3,500 genes are up-and downregulated at each time point, of which 1% display sex-biased expression.(E) From D0 to D4: 2,218 genes are upregulated (26 in XY and 24 in XX with a sex-bias) and 3,135 are downregulated (23 in XY and 52 in XX with a sex-bias); (F) From D0 to D9: 3,566 genes are upregulated (0 in XY and 34 in XX with a sex-bias) and 2,961 are downregulated (30 in XY and 9 in XX with a sex-bias); (G) From D0 to D37: 5,518 are upregulated (58 in XY and 17 in XX with a sex-bias) and 5,456 are downregulated (12 in XY and 21 in XX with a sex-bias).
lines.Before differentiation (D0), upregulated gene set clusters in male cell lines were associated with ribosomal and mitochondrial functions (Supplementary Figure S5).At D4, a substantial number of actin-based processes were increased in male cells, while in female cells, cilium-associated processes, as well as DNA packaging and chromatin organization, were upregulated (Figures 4A, C).At D37, male cells exhibited upregulation in cilium and microtubule-based processes, as well as in various neural processes such as synaptic signaling, synaptic vesicle transport, axon structure, axon guidance, neuron recognition, and action potential.In contrast, female cell lines were upregulated in cell division and ribosomal processes (Figures 4B, D).Details for all enriched gene sets within a cluster can be found in Supplementary Table S7.
At D37, the GSEA identified three clusters of sex-biased neural processes in male cell lines.To pinpoint the specific genes contributing to these clusters, consisting of 22 individual gene sets, we extracted the genes with the highest contribution to the gene set enrichment, known as core-enriched genes (Figure 5A).Enriched gene sets at D4 and D37: Gene set enrichment analysis was based on differentially expressed genes at D4 and D37.Enriched gene sets in male and female cell lines are colored in blue and red, respectively.Panels (A, B) display an enrichment map in which a dot represents an upregulated gene set.If gene sets are connected, and the closer they are in proximity, the more they are related in function.Clusters of gene sets are encircled and named after their common biological function.(A) Enrichment map for D4: male cell lines display upregulated gene sets that are associated with actinbased processes, such as "Cell-substrate adhesion", "Anchoring junction" and "Actin cytoskeleton and actin filament organization".Female cell lines at D4 display larger gene set clusters of cilium-based processes as well as "DNA packaging and chromatin organization".(B) Enrichment map for D37: male cell lines show three clusters of enriched gene sets that are associated with neuronal processes (cluster names highlighted in blue) and a larger cluster that combines cilium and microtubule-based processes.Female cell lines at D37 show an upregulation of gene sets involved in cell division and ribosomal processes.Subfigures (C) and (D) summarize upregulated gene sets in male and female cell lines (at D4 and D37) and display the number of gene sets within a cluster.Enriched neuronal gene sets and corresponding core-enriched genes in male cell lines at D37: Core-enriched genes from enriched neuronal gene sets in male lines at D37 were contributing to candidate gene selection.(A) Gene sets upregulated within the neuronal clusters.Numbers in brackets show the total number of genes within the gene set.Circle size illustrates the percentage of core-enriched (highly contributing) genes within the gene sets.Color indicates the false-discovery rate (q-value) of the gene set at comparison of male and female cell lines.(B) Gene expression heatmap of coreenriched genes from all gene sets within the sex-biased neural process clusters.Colors indicate the level of expression comparing male and female cell lines based on a Z-score.Red indicates an overexpression, black indicates an equal expression in both cell lines, green indicates a reduced expression.(C) Using stringent filtering, we identified 17 core-enriched and significantly sex-biased genes within the clusters.These genes later contributed to the selection of candidate genes.
This resulted in 362 unique core-enriched genes (Figure 5B), which we further filtered using previously mentioned stringent filtering rules, to identify 17 sex-biased genes implicated in neural processes (Figure 5C).The full list of core-enriched genes is available in Supplementary Table S8, with the 17 sex-biased genes highlighted.
Combined information returns candidate genes that show robust sex-biased overexpression in the cell lines and are implicated in neural processes Next, we aimed to identify candidate genes demonstrating robust sex differences in expression and implication in the development and function of neurons.For this, we combined all the information gathered in this study to annotate the sex-biased genes that resulted from the DEG analysis at D37 (104 upregulated genes in XY cells and 44 in XX cells).The following information was used for the annotation: (a) GO-terms associated with neurodevelopment, (b) upregulation at D37 of differentiation, (c) appearance in clusters of neuronal processes in the GSEA, (d) indications from the literature for involvement in neurodevelopmental processes, and (e) the 17 sex-biased genes resulting from the GSEA analysis of D37.Based on these annotations more than 60 of the 148 genes display implications in neuronal development.Among these genes, 13 displayed an expression level above the 70th percentile, which we decided to select as candidate genes.Ten of these genes display an upregulated gene expression in male cell lines (NHLH2, EBF1, SLC17A6, RUNX1T1, KIF5A, AKAP12, MDGA1, ONECUT2, P2RX3, LMX1B), and 3 genes are upregulated in female cell lines (SYAP1, AMOT, PAK3), as detailed in Table 3.The full list of annotated genes can be found in Supplementary Table S5.To validate the expression tendencies of the candidate genes, we conducted qPCR analyses on material from a replicate differentiation experiment, including two additional cell lines for each sex (KARO1-XX and HS1001-XY) and two extra timepoints, D17 and D27.By including timepoints D17 and D27, we observed that most candidate genes display an initial expression earliest at D17 (Figure 6).Considering this and the onset of DCX expression at D17 (Supplementary Figure S2), we regard D17 as the earliest timepoint where neurons begin to emerge in our differentiation experiment.The sex-biased expression of candidate genes was affirmed at D37 in all cell lines, with some genes even exhibiting a sex difference as early as D27.The difference between the male and female samples was tested using a repeated measurements two-way ANOVA with Šidák's multiple comparisons test.
Gene dosage compensation of X/Y homologs escaping XCI through Y homolog and significant Y homolog overexpression in TXLNG/Y and KDM6A/UTY Gametologous genes are X/Y homologous genes located on opposite sex chromosomes in non-recombining regions.Gene dosage of gametologous genes can contribute to sex-biased gene expression.The analysis of these genes in male and female cell lines has revealed several sex differences (Table 4).Since the impact of XCI erosion on the gene expression of X-linked genes cannot be entirely ruled out, the interpretation of gametologous gene expression should be done with caution, particularly when there is a significantly increase of expression in the X homologs in this section.We categorized the gametologous gene pairs into four groups based on their expression pattern during the neural differentiation period from D0 to D37.Group A includes genes that are significantly downregulated at D37 (GYG2/P1, ARSL/P1, ANOS1/P2, GPR143/P, RBMX/Y1A1; Figure 7A).In contrast, Group B comprises genes upregulated at D37 (ARDS/P1, MXRA5/Y, SHROOM2/P1, TSPYL2/TSPY, CASK/P1; Figure 7B), while Group C shows no change in gene expression from D0 to D37 (STS/P1, TMSB4X/Y, BCOR/1, SOX3/SRY, OFD1X/P1Y; Supplementary Figure S6).Group D includes all gene pairs that show substantial Y chromosome homolog expression contributing to gene dosage compensation (ZFX/Y, USP9X/Y, DDX3X/Y, KDM5C/D, RPS4X/Y1, NLGN4X/Y, PRKX/Y, EIF1AX/Y, PCDH11X/Y, TLX1X/Y; Figure 8A).This group includes two gene pairs (TXLNG/Y, KDM6A/UTY) that display a significant overexpression of the Y chromosome homolog leading to sex biases in gene dosage, Figure 8B.More than half of the X-linked genes of the X/Y homologous gene pairs (19 out of 27) are escaping XCI according to Katsir and Linial (2019), Garieri et al. (2018) and Zhang et al. (2013).We observe that for many of the XCI-escapees (12 out of 19), the Y chromosome compensates the increased X chromosome expression, and therefore, these gene pairs do not contribute to sex differences.These genes were clustered in Group D, with the exception of one gene (ANOS1), which sharply drops in gene expression at D37 and thus is categorized into Group A.
We detected a female-biased overexpression of gametologous genes in Group A (GYG2, ANOS1, BCOR1) during early differentiation (D0, D4, D9).This sex bias could contribute to different prerequisites that affect subsequent differentiation steps.We also noticed a female-biased overexpression in genes of Group B: ARSD at D0, D4, and D37 and MXRA5 at D37.A sex bias in Group B at D37 suggests relevance in late neural differentiation or maintenance.Of the Group C genes, the gene RPS4X is significantly overexpressed in females at D37 due to insufficient Y chromosome gene dosage compensation.Moreover, the Y homolog expression of the Group D gene pairs TXLNG/TXLNGY and KDM6A/UTY is very high at nearly all timepoints, leading to a significant increased gene dosage in male cell lines.

Discussion
Understanding sex differences in the field of neuroscience is important to predict, monitor and treat neurological disorders, especially since a majority of early-onset neurodevelopmental disorders show a clear sex bias (Boyle et al., 2011;Loke et al., 2015;Tesic et al., 2019).In this study, we sought to determine sexspecific gene expression during early neural differentiation.We used four male and four female-derived hESC lines and applied a differentiation protocol resulting in a population of mixed neurons.Using hESCs as a research model has several advantages.By using epiblast-derived embryonic stem cells, we aimed to keep the epigenetic variability low, that can otherwise arise as a result of "the epigenetic memory" of the cell-of-origin when cells are reprogrammed to an induced pluripotent state (Narsinh et al., 2011).Furthermore, we can exclude the effects of sex hormones by using a cell culture model in a controlled in vitro environment.Thus, we consider our model a valuable system to study genetic sex differences that affect early nervous system development.It proved to be a robust system as it resulted in reproducible and comparable differentiation in all cell lines (Figures 1B-E).While our focus was primarily on neurons, it is important to acknowledge the complexity of neural development, which encompasses various cell types such as astrocytes, oligodendrocytes and microglia.Our findings demonstrate that under the culture conditions employed, neuronal development was favoured, with no discernible markers for other cell types within the neuronal lineage detected.While it's important to interpret these findings thoughtfully, it's worth noting that the culture conditions utilized, though not exact replicas of organotypic environments, still provide valuable insights contributing significantly to our understanding of sex-specific mechanisms in neuronal differentiation.While X-linked genes have already been extensively linked to the development and functions of the nervous system (Skuse, 2005;Nguyen and Disteche, 2006;Mallard et al., 2021;Raznahan and Disteche, 2021;Cabrera Zapata et al., 2022), evidence for the involvement of Y-linked genes in male neuronal development is scarce (Vakilian et al., 2015;Meyfour et al., 2017;Johansson et al., 2019;Pottmeier et al., 2020;Cabrera Zapata et al., 2022).In this study, we confirm the involvement of X-and Y-linked genes in the neuronal differentiation of male and female human embryonic stem cells.Furthermore, we provide evidence for sex-biased gene expression during neural differentiation, most prominent after 37 days of differentiation, and highlight 13 candidate genes that show sex-biased expression, thus possibly impacting neuronal development.Finally, we display gene dosage compensation of X homologs that escape X-inactivation through their Y homolog in male cell lines suggesting a more important role for the Y chromosome in neurodevelopment than generally believed.
The comparison of gene expression in male and female stem cell lines presents a number of challenges.One of them, which specifically affects X-linked gene expression is the state of XCI in Candidate gene expression at all time points: Gene expression was measured by qPCR in replicates of the neural differentiation experiment (n = 3).All candidate genes (MDGA1, NHLH2, LMX1B, ONECUT2, EBF1, KIF5A, P2RX3, AKAP12, SLC17A6, RUNX1T1) show significant overexpression in all four male cell lines at D37.Many genes start expressing at D17. RUNX1T1 shows a significant sex difference already at D9. ONECUT2 shows a female-biased sex difference at D27, which then changes to an overexpression in male cells at D37. Differences between the male and female samples were measured using repeated measurements two-way ANOVA with Šidák's multiple comparisons test.
female stem cell lines.It is well-documented that commonly used cell culture conditions lead to the erosion of XCI in female cell lines (Hanna et al., 2010;Cloutier et al., 2022), which in turn can result in increased expression of X-linked and autosomal genes.The absence of XIST expression in two of the female cell lines and the number of genes with an XX:XY expression ratio above 1.1 indicates different XCI erosion statuses in our cell lines.However, despite this erosion, the very similar expression of X-linked genes and the low number of X chromosome genes with an XX:XY expression ratio above 1.1 indicates that the XCI erosion has only a minor effect on their X-linked gene expression.The same is true for the autosomal gene expression as the comparison of the average X:A ratio revealed no significant increase in the XIST positive cell line.In addition, our female cell lines were able to differentiate equally well to neurons, whereas under XCI erosion cell lines have been reported to be impaired in neural differentiation (Motosugi et al., 2022).Moreover, we must consider whether our analyses to detect sex differences in the study might be influenced by genes exhibiting high variance, possibly attributed to the erosion of XCI.Here, it's worth noting that such genes were excluded based on the stringent filtering criteria we implemented for identifying sex-biased genes.Consequently, we assert that the findings of our study remain unaffected by the varying XCI erosion statuses observed in the XX cell lines.The use of the XX:XY ratio to estimate the number of genes affected by XCI erosion is an approximation and may not cover all cases.The observed disparities for chromosomes X, 5, and 16 reflect a combination of cell-line-specific variability and authentic effects of XCI erosion, which are inherently intertwined and indistinguishable.Certain scenarios may remain unaccounted for, such as instances where physiologically low X-linked gene expression in a female cell line, yielding an XX:XY ratio of 0.5, is subsequently doubled due to the reactivation of the second X chromosome through XCI erosion.In this scenario, although the gene is affected, it does not meet the criteria for being considered "affected by XCI erosion" (ratio>1.1).This methodology assumes uniform X-linked gene expression between female and male cells under identical experimental conditions.We consider this approach suitable to estimate the general progression of XCI erosion and its impact on X-linked gene expression in female cell lines, given that average X-linked expression levels are commonly reported to be elevated in XIST -hPSCs compared to XIST + and male hPSCs (Bar et al., 2019;Motosugi et al., 2022;Sarel-Gallily and Benvenisty, 2022).Nevertheless, results including X-linked genes should be interpreted with caution as a possible effect from XCI erosion cannot be entirely excluded.
Considering that the majority of gametologous genes are XCI escape genes (19 out of 27), transcribed from two X chromosomes, it is improbable that XCI erosion significantly affects their gene expression.Only two gametologous genes (CASK and SHROOM2) exhibit signs of XCI erosion, as indicated by an XX: XY ratio above 1.1 and are not XCI escape genes.Therefore, we consider the analysis of gametologous genes presented in this study to be largely unaffected by XCI erosion.

Genetic sex differences are not routinely investigated
Even though 'sex differences' is a subject that receives increasing awareness among researchers and is regularly debated among neuroscientists (DeCasien et al., 2022), the investigation of sexbiased gene expression is not pursued frequently.Most studies that investigate sex-biased gene expression during neurodevelopment are transcriptomic studies that utilize brain tissues of the developing embryo (Shi et al., 2016;O'Brien et al., 2018;Johansson et al., 2019; # marks genes escaping XCI, XX (red) and XY (blue) displays increased gene dosage in female or male hESC, lines, respectively.
Frontiers in Cell and Developmental Biology frontiersin.orgGametologous genes with low Y homolog expression: X/Y homolog expression was analyzed in total RNA sequencing data (3 male and 3 female cell lines).A hashtag after the X homolog indicates that the gene is escaping XCI.Panel (A) displays the 5 gene pairs that decrease (GYG2/P1, ARSL/P1, ANOS1/ P2, GPR143/P, RBMX/Y1A1), and (B) the 5 gene pairs that increase (ARDS/P1, MXRA5/Y, SHROOM2/P1, TSPYL2/TSPY, CASK/P1) in expression during neural differentiation, comparing the period D0 to D37.An increased gene dosage in females was found in GYG2/GYG2P1 at D0, D4, and D9 (p < 0.0001, p < 0.001, p < 0.001, respectively) and ANOS1 at D0 and D4 (p < 0.001, p < 0.05).Differences between the male and female samples were measured using repeated measurements two-way ANOVA with Šidák's multiple comparisons test.For TXLNG/Y, this led to a male bias in gene dosage at all time points D0, D4, D9, and D37 (p < 0.05, p < 0.0001, p < 0.0001, p < 0.001, respectively).For the gene pair KDM6A/UTY, a male-biased gene dosage was found in D4, D9, and D37 (p < 0.01, p < 0.01, p < 0.001, respectively).Differences between the male and female samples were measured using repeated measurements two-way ANOVA with Šidák's multiple comparisons test.Liu et al., 2020;Oliva et al., 2020;de Toledo et al., 2022).We have only found one recent study that, like ours, uses human stem cell lines to investigate gene expression during neural differentiation (Waldhorn et al., 2022).Unfortunately, they only used 1 cell line per sex.Many studies that use multiple human stem cells for neuronal differentiation investigate differentiation characteristics and mechanisms but do not include or consider the sex of the used cell lines (Maroof et al., 2013;Tsankov et al., 2015;Bardy et al., 2016;Hoffman et al., 2017;Schwartzentruber et al., 2018;Burke et al., 2020;Kathuria et al., 2020;Lindhout et al., 2020;Strano et al., 2020;Kuruş et al., 2021).Even large-scale cell bank initiatives that meticulously characterize their stem cells do not investigate differences or similarities between the cell lines of different sexes (Salomonis et al., 2016;Kilpinen et al., 2017;Panopoulos et al., 2017;D'Antonio et al., 2017).This underlines that sex differences are still not routinely investigated even though the transcriptomic data is at hand.A large number of studies use only a single cell line for the investigation of neuronal differentiation properties, neglecting to account for interindividual differences of cell lines and missing the chance to investigate sex differences (Wu et al., 2010;van de Leemput et al., 2014;Chu et al., 2016;Close et al., 2017;Song et al., 2017;Meisig et al., 2020;Allison et al., 2021;Solomon et al., 2021).
Y-linked genes in neuronal development: Insights from the differentiation of male embryonal stem cells It has been reported that a major contribution to sex differences arises from the Y chromosome genes or the dosage of X chromosome genes (Jazin and Cahill, 2010).The Y-linked genes detected during our differentiation experiment seem to be robustly expressed since most of them were also found in undifferentiated male human iPSCs (13 out of 16, compared with D0) as well as after 14 days of neuronal differentiation (14 out of 17, compared with D9 and D37) (Waldhorn et al., 2022).Based on their gene expression patterns, a number of Y-linked genes (Figure 2A; Supplementary Table S4) are likely to be involved in the development of neurons in male stem cell lines.The peaking expression of the genes TXLNGY, NLGN4Y and ZFY at D4 indicates a role in early neuronal induction.Interestingly, TXLNGY is considered a pseudogene but shows expression in multiple tissues, including the brain.To date, the function of TXLNGY is still unknown.Expression evidence from our present and a previous study (Pottmeier et al., 2020) both suggest an involvement of TXLNGY in processes of early neuronal induction.The gene NLGN4Y is a member of the neuroligin family, which is important for the formation and function of synapses.Due to the high sequence similarity of NLGN4X and Y in humans (97%), the two homologs are assumed to have the same function.A recent study, however, showed that a single amino acid difference in NLNG4Y results in a deficit in trafficking that leads to decreased surface expression and thus possibly to decreased synaptic function of NLGN4Y (Nguyen et al., 2020).Nevertheless, NLGN4Y has been reliably detected in the human brain (Trabzuni et al., 2013;Tukiainen et al., 2017;Johansson et al., 2019) and during neuronal differentiation of stem cells (Pottmeier et al., 2020).Our results support the claims for an involvement of NLGN4Y in neuronal development and suggest a role in early neuronal induction.The gene ZFY encodes a zinc finger transcription factor that is able to regulate gene expression and is expressed in the developing and mature human brain (Mayer et al., 1998;Weickert et al., 2009;Shi et al., 2016).The role of ZFY in neuronal differentiation is unknown.In our experiment, ZFY is peaking at D4 of neuronal differentiation, indicating a possible contribution to the early neuronal induction of male stem cells.The gene RPS4Y1 is the only Y-linked gene peaking in expression at D9. RPS4Y1 is a ribosomal protein and thus involved in the synthesis of proteins.Its function in neuronal development is not fully understood; however, it is expressed in the developing and mature brain (Weickert et al., 2009;Kang et al., 2011;Shi et al., 2016;Khani et al., 2022).Moreover, it is reported to play a role in neuronal differentiation (Vakilian et al., 2015), which is supported by our results, specifically in an advanced stage of neuronal induction.The Y-linked genes with the highest expression at D37 are USP9Y, UTY, KDM5D and PCDH11Y.Due to their peaking expression at this timepoint, these genes are likely to play an important role at this stage of differentiation.The gene USP9Y encodes a deubiquitinase and is expressed in multiple tissues, including the developing and mature brain (Weickert et al., 2009;Kang et al., 2011;Shi et al., 2016;Johansson et al., 2019).Ubiquitylation orchestrates core signaling networks essential for stem cell maintenance and differentiation.Consistent with our findings, there is evidence suggesting that USP9Y plays a role in neuronal differentiation (Vakilian et al., 2015).PCDH11Y encodes a member of the protocadherin protein family.Due to their high sequence similarity (98%), PCDH11X and its Y chromosome homolog are believed to perform a similar function.Many protocadherins are predominantly expressed in the brain and play a role in cell-cell adhesion and signaling, specifically in developing neurons (Pancho et al., 2020).PCDH11Y is expressed in developing and adult brain tissue (Mayer et al., 1998;Weickert et al., 2009;Shi et al., 2016) and is thought to play a role in neuronal development and function through its contribution to the formation and maintenance of synapses (Durand et al., 2006).The large increase of USP9Y and PCDH11Y expression at D37 of the differentiation experiment suggests a role in neuronal maturation or maintenance.

Histone demethylases of the Y chromosome in neurodevelopment
Histone modifications play critical roles in a plethora of biological processes, including neuronal development.We found that two Y-linked histone demethylases, UTY (KDM6C), and KDM5D (JARID1D), are highly expressed in male human embryonic stem cells after 37 days of neuronal differentiation (Figure 2A), indicating a contribution to neuronal development.Histone demethylases are particularly interesting for gene expression because they can regulate transcription initiation, elongation, and termination through post-translational histone modification.The gene UTY (KDM6C) and its X chromosome homolog UTX (KDM6A) are members of the histone lysine demethylase subfamily 6 (KDM6).It has not been fully clarified if UTY possesses catalytical activity.There are studies confirming (Walport et al., 2014) but also denying (Hong et al., 2007;Lan et al., 2007;Shpargel et al., 2012) UTY its catalytical ability.UTY seems to have reduced catalytic activity compared to UTX, but has the same protein binding capacity (Walport et al., 2014).It is possible that UTY has access to additional cofactors in vivo that support its histone demethylase activity, explaining its reduced enzymatic activity in vitro.Apart from their demethylating function, it has also been suggested that UTY and UTX possess a histone demethylase-independent function, as the presence of UTY allows UTX-deficient mouse embryos to survive until birth (Shpargel et al., 2012).UTY and UTX share 88% sequence identity, and both are H3K27-specific demethylases that can supports the activation of gene transcription (Hong et al., 2007;Lee et al., 2007;Swigut and Wysocka, 2007).It has been demonstrated that the regulation of H3K27 methylation states is essential for the development and function of the mammalian nervous system and that enzymes involved, such as UTX and likely also UTY, are key regulators of human neural differentiation in stem cells (Lei and Jiao, 2018;Yang et al., 2019;Shan et al., 2020;Tang et al., 2020).Our study supports the previous claims suggesting a contribution of UTY to neural differentiation.
The gene KDM5D (JARID1D or SMCY) and its X chromosome homolog KDM5C (JARID1C or SMCX) are members of the histone lysine demethylase subfamily 5 (KDM5).KDM5D and KDM5C show an 85% sequence similarity.They exhibit conserved protein domains for histone demethylase activity (JmjC), but experimentally determined H3K4 demethylase was lower for KDM5D than KDM5C (Iwase et al., 2007).A demethylation at H3K4 is considered as a sign of repressed gene transcription (Barski et al., 2007).H3K4 methylation states are regulated by several SET-domaincontaining enzymes, including the KDM5/JARID1 family of histone demethylases (Cloos et al., 2008).A dysregulation of H3K4-specific methylation states has been linked to various developmental disorders and impaired cognition (Bale et al., 2010;Millan, 2013;Ronan et al., 2013;Kim et al., 2017).There is evidence that H3K4 methylation landscapes in the human cerebral cortex are dynamically regulated during prenatal development and throughout early childhood years until adolescence (Shulha et al., 2013).Furthermore, it has been shown that neuron-specific H3K4me3 peaks exist in humans, and that they are enriched at promoters of genes that control synaptic function (Cheung et al., 2010).Mutations in the X chromosome demethylase KDM5C are associated with X-linked intellectual disability and epilepsy (Strømme et al., 2002;Najmabadi et al., 2011;Poeta et al., 2013).Further, KDM5C is implicated in the regulation of neural fate commitment in stem cells (Dey et al., 2008;Kim et al., 2018).Y-linked demethylases are among the least studied enzymes involved in histone modifications, most likely because they are thought to have a lower catalytic activity than their X-linked counterparts.In our present study, we show that the Y-linked demethylases KDM5D and UTY are increasingly expressed during neuronal differentiation of human embryonic stem cells.These results are in line with a previous study, in which we have observed an increased expression of KDM5D and UTY during differentiation of neural stem cells (Pottmeier et al., 2020).Our studies suggest a role of Y-linked demethylases in male-specific neural development.Due to their conserved protein domains, it is likely that the X and Y chromosome homologous demethylases perform a similarly important function in human neural development.

Hormone-independent genome-wide gene expression changes during neuronal differentiation
At all timepoints investigated, including the ESC stage at D0, we found the number of sex-biased genes to be increasing as differentiation progressed (Table 2).This indicates that sex differences in gene expression become more pronounced as the nervous system develops.Together with the increasing influence of sex hormones during development, this could be a major factor contributing to the emergence of sex-specific neural circuits in human brain development.Although the number of differentially expressed genes at D0 was small (65), robust sex differences were already present, and interestingly, they originated to a large part from the Y chromosome in male ESCs (16 DEG) and from the X chromosome in female ESCs (21 DEG).A gene set enrichment analysis found that male cells at this early stage are enriched for ribosomal and mitochondrial functions (Supplementary Figure S5).This could point towards increased cell proliferation and differentiation in the male cell lines.Although we have not measured the growth rate in the used hESC lines, this is in line with previous observations where we have seen increased cell proliferation in male neural stem cells (Pottmeier et al., 2020).While the number of sex chromosome genes with sex-biased expression remained constant during differentiation, the number of autosomal genes with sex-biased expression increased in both sexes during differentiation, but to a greater extent and with larger variations in male cells (Table 2).
A sex-bias in autosomal genes has been reported in a number of transcriptomic studies of neurodevelopment (Trabzuni et al., 2013;Ronen and Benvenisty, 2014;de Toledo et al., 2022).In male cells at D37, the number of autosomal genes contributes the most to the sex bias in expression (86%, 89 of 104 genes), while in the female cells, a similar number of autosomal (23) and X-linked genes (47%, 21 of 44) are displaying a sex bias in expression.It has been demonstrated that the sex chromosome complements and especially the Y-linked gene SRY has an effect on global gene expression.The effect is believed to be achieved by the direct activation of key autosomal genes and the autocrine or paracrine effects of sex steroids (Wijchers et al., 2010;Ronen and Benvenisty, 2014).Based on this, it is likely that a gene expression amplifying effect by the Y chromosome is also present during differentiation and could therefore explain the high number of sex-biased autosomal genes during differentiation in our male cell lines.From D0 to D4, the number of upregulated autosomal genes increased by the factor 10 in males and by the factor 3 in females, affecting biological processes related to actin dynamics in male cells and cilia, chromatin and DNA packaging in female cells, according to GSEA (Figure 4).This displays that sex differences in gene expression are already present as early as 4 days after differentiation, and that they are able to affect biological processes relevant for neuronal development (Figure 4).At D9, no gene sets were enriched in male or female cell; however, females displayed 5 times more autosomal DEG than males.At D37, the most advanced stage in our neurodevelopment experiment, we see a total of 148 sex-biased DEG from which a majority (70%) was upregulated in males.Of these DEG, 39 were implicated in neuronal processes, with the majority involved in neurodevelopment.In total, four times more DEG were upregulated in male than in female cell lines.These results are supported by Waldhorn et al. (2022) who also found the majority of DEG to arise from male cells in neuronal differentiation of hiPSC, as well as by Shi et al. (2016) providing evidence for a male-driven sex bias in the transcriptome of the human fetal brain.Our data illustrates how early sex chromosome activity can set the stage for a later manifestation of differences and suggests a larger effect on the male cell lines.While in our study the male DEG at D37 were specifically associated with neurogenesis, cell-cell signaling, synaptic transmission, trans-synaptic signaling and neuronal differentiation, transcriptome investigations of human embryos revealed contrasting results.In two other studies, female samples were enriched for neuronal differentiation and synapse specification, as well as synaptic processes and synaptic signaling, respectively (O'Brien et al., 2018;de Toledo et al., 2022).These discrepancies can be explained by the nature of the samples and the difference in stages of neuronal development, as data from the above-mentioned studies is gathered from mid-gestation embryos, that represent a more advanced stage than our hESC cells after 37 days of differentiation.A limitation of our study is the lack of RNA sequencing data from additional timepoints during differentiation.Even though data was collected over a number of timepoints until the endpoint at D37, for the interpretation of results from the latest stage it has to be considered that possible delays or shifts in timing in neuron development between the sexes cannot be resolved.Moreover, the results are retrieved from pooled material from a mixture of cells and thereby limiting conclusions on the single-cell or cell-type basis.
Among the sex-biased DEG, we found a considerable number of X-linked genes throughout all differentiation stages (Supplementary Table S1).Most notably are the 12 genes ARSD, CTPS2, GYG2, LRCH2, NAP1L3, RBBP7, SYAP1, TXLNG, ZFX, ZRSR2, PLS3, and PABPC5 which show a significant overexpression in female cell lines in every investigated timepoint (D0, D4, D9, D37).All of these genes, except PLS3, NAPL1L3, PABPC5, and LRCH2, have been shown to stably escape X inactivation (Katsir and Linial, 2019), explaining the consistent sex bias in expression during differentiation.PLS3, NAP1L3, and PABPC5 have been described as variable XCI, dependent on the investigated tissue type (Garieri et al., 2018).The genes TXLNG and ZFX are gametologous genes that are significantly overexpressed genes in females.However, we have shown that their Y-linked homolog compensates for the overexpression (Figure 8) leading to relatively balanced gene dosage.Of the above-mentioned X-linked genes with a significant overexpression in females PLS3, SYAP1, TXLNG, ZFX show an implication in neurodevelopment according to GO-terms and the genes ARDS, LRCH2, and ZRSR2 are peaking in expression at D37, also suggesting an involvement in neuronal differentiation.Interestingly, the X-linked genes GABRE and AC244197.3show a significant overexpression in male cell lines at D37. AC244197.3encodes a protein of unknown function.Since AC244197.3significantly increases in gene expression at D37, it suggests an involvement in male-specific neuronal differentiation or maintenance.GABRE encodes the epsilon subunit of the gammaaminobutyric acid A receptor.During neurodevelopment, GABA plays a critical role in regulating the proliferation, differentiation, and migration of neural cells.GABA also plays a role in the formation and maturation of synapses.GABA A receptor is constructed as a pentameric structure of multiple subunits (α, β, γ, δ, ε, θ, π).Sex differences in GABA receptor subtypes were already reported in adult human brain tissue, displaying a higher level of GABA A receptor subunits α1, α2, α5, and β3 in male brain tissue (Pandya et al., 2019).

Compensation of sex differences by gametologous genes
Closer investigation of our RNAseq data for gene dosage in gametologous genes (Figures 7, 8; Supplementary Figure S6) revealed that seven gene pairs showed a higher gene dosage in female cells, with GYG2/P1 and ARSD/P1 higher expressed over several stages, ANOS1/P2 and BCOR/1 at early, and MXRA5/Y, RPS4X/Y1, and SOX3/SRY at the last stage, D37.Two gametologous pairs, TXLNG/Y and KDM6A/UTY showed a higher gene dosage in male cells at all differentiation timepoints due to the contribution of the Y homolog.TXLNG/Y was also overexpressed in undifferentiated ESCs (D0).In a previous study, we found UTY and also KDM5D expression increased upon differentiation in male NSCs (Pottmeier et al., 2020).Interestingly, the gene dosage of the X homolog KDM5C is compensated by the Y homolog even in the undifferentiated ESC stage.The lysine demethylases KDM6A target the Hox gene family and thereby control pluripotency and lineagespecific genes and can be partly compensated by the complement UTY.Despite or rather because of the critical function in development, KDM6A/UTY might be an important contributor to the genetically driven sexual dimorphisms.

Candidate sex-biased genes implicated in neurodevelopment
We annotated the sex-biased DEG at D37 with abundant expression using the following data: GO-terms related to neurodevelopment, upregulation at 37 days of differentiation, enrichment in clusters of neuronal processes in GSEA, a minimum expression level of 500 DESeq2 normalized counts, and presence in literature associated with neurodevelopmental processes (Supplementary Table S5).This has resulted in 10 candidate genes that show upregulated expression in male cell lines (NHLH2, EBF1, SLC17A6, RUNX1T1, KIF5A, AKAP12, MDGA1, ONECUT2, P2RX3, LMX1B) and 3 candidate genes with upregulated expression in female cell lines (SYAP1, AMOT, PAK3), Table 3.The observed sex differences in candidate genes have proved to be highly replicable even across gene expression analysis methods, demonstrated by the consistent results of RNA sequencing (Supplementary Figure S7) and qPCR analysis of replicated differentiation experiments (Figure 6).Noticeably, all of the genes with higher expression in male cell lines are autosomal genes and half of them are brain-specific transcription factors (NHLH2, EBF1, RUNX1T1, ONECUT2, LMX1B) with critical functions in the development and differentiation of neurons.In addition, all of the candidate genes in males are implicated in neurodevelopmental and/or neurological diseases (Huynh et al., 2012;Villaescusa et al., 2016;Afshar et al., 2020;Sherman et al., 2022;Flex et al., 2023;Kimura et al., 2023;Verma et al., 2023).The three candidate genes that are upregulated in female cells are X-linked genes, highlighting the importance of the X-chromosome in neurodevelopment.All of these genes are known to be involved in neurological diseases; AMOT is implicated in Alzheimer's disease and PAK3 is associated with inherited X-linked intellectual disability (Sherman et al., 2022).The human gene encoding SYAP1 is located within chromosomal band Xp22.2, a region associated with mental retardation, developmental delay, and autism spectrum disorder (Sismani et al., 2011;Prasad et al., 2012).The candidate genes display varying degrees of conservation; all are conserved among mammals, with confirmed expression and function in the mouse nervous system.Additionally, some genes show high conservation with orthologues not only in mice but also in other model organisms.For instance, LMX1B, KIF5A, and P2RX3 have orthologues in zebrafish, while ONECUT2, RUNX1T, KIF5A, SLC17A6, SYAP1, and PAK3 have orthologues in drosophila.Moreover, LMX1B, SLC17A6, and PAK3 have orthologues in C. elegans.All these genes exhibit confirmed expression and a role in the development and/or function of the nervous system across multiple species (Wigerius et al., 2020;Bult and Sternberg, 2023), likely suggesting that they have roles fundamental in pathways conserved across a wide range of organisms.
Strikingly, not only the number of DEG in male cells is higher, there are also a considerable number of transcription factors among the candidate genes.In contrast, the candidate genes in female cells are associated with cellular processes and signaling while transcription factors are not overrepresented.This could point towards an active process in male cells based on intrinsic genetic factors that contribute to early programming or priming of the neural cells.Notably, as mentioned earlier, the direct influence of SRY on global gene expression emerges as a significant genetic contributor to sex-specific expression patterns.Epigenetic modifiers, such as histone demethylases like KDM6A, KDM5D, and UTY, have been identified throughout the differentiation phase, with UTY exhibiting markedly elevated levels in male cells compared to female cells.This discrepancy potentially impacts the chromatin landscape surrounding the promoters and enhancers of candidate genes.Using databases at The Biological General Repository for Interaction Datasets (BioGRID4.4https:// thebiogrid.org/)we uncovered several documented direct or mediated protein-protein interactions between epigenetic modifiers and DEG: UTY displayed interactions with KDM6B, ANOS, KIF-paralogues, and USP9X, while KDM5D exhibited interactions with AMOT, Runx1, KDM6A, and NHLHparalogues.TXLNGY was found to interact with AMOT and KIF-paralogues via TRIM14.Additionally, evidence of interactions between KDM6A and AMOT, NHLH-and KIFparalogues, USP9, and EIF1AY was observed.Genetically, interactions were documented between KDM6A and KDM5A and 5C.Moreover, KDM6A was observed to directly interact with the pioneer transcription factor FOXA1, as well as with PAX, LHX, GATA, and SOX family factors.These findings suggest a complex involvement of epigenetic regulatory networks, intertwined with fundamental neurodevelopmental pathways, influencing the expression of sex-biased genes.However, it remains to be determined whether UTY, sharing 82% sequence homology with its gametologue KDM6A, shares the same interaction partners.The significant difference in the number of known interactors between X-gametologues, particularly KDM6A (with 158 interactors), and the Y-gametologue UTY (with 21 interactors), warrants further investigation into the underlying mechanisms governing sexspecific gene expression.The fewer interactors found for Y-linked gene products may stem from a lack of comprehensive information, potentially due to the overlooked roles of Y-chromosome genes.These sex differences in early programming could set the stage for later differentiation events and contribute to the establishment of a neural network distinct from that in females.Additionally, the genetic differences observed could prepare for later responsiveness to hormonal cues during further development or in response to external stimuli.The observed expression bias may contribute not only to physiological sex-specific differences but also to varying susceptibility to neurological diseases as deviations in proliferation and differentiation have been associated with neurodevelopmental disorders (Ernst, 2016) that have a higher prevalence in males such as autism spectrum disorder, schizophrenia, attention-deficit hyperactivity disorder have a higher prevalence in males (Saha et al., 2005;May et al., 2019;Maenner et al., 2020).
Further research is needed to elucidate the specific molecular mechanisms underlying these observations and to understand how they contribute to the development of sexually dimorphic characteristics in the nervous system.

Conclusion
In conclusion, our study identified a very early genetically driven and hormone-independent trait separation according to sex in the neurodifferentiation of human embryonic stem cells.We provide a catalog of genes that display a large degree of sex differences in expression during neuronal differentiation, which holds significant implications for understanding the origins of sex-biased neurodevelopmental diseases.The most substantial contribution was found to arise from the male transcriptome, involving both Y chromosome and autosomal genes, of which the majority have crucial functions in neurodevelopment and differentiation.It is important to note that mRNA and protein levels do not always correlate, and additional layers of post-transcriptional regulation are likely influenced by sex as well.Therefore, we advocate for further investigation into the direct impact of sex on protein expression.To fully validate these findings, additional research is warranted to confirm the observed sex biases and explore the potential functional roles of candidate genes in neuronal differentiation.Ultimately, we hope that our study contributes to the understanding of genetic sex differences in neurodevelopment and sheds light on their potential implications in the development of sex-biased neurological disorders.

Data availability statement
The data presented in the study are deposited in the Gene Expression Omnibus (GEO) repository, accession number GSE249035.

FIGURE 8
FIGURE 8 Gametologous genes with high Y homolog expression: Gene expression was analyzed based on total RNA sequencing data.The hashtag symbol after the X homolog indicates that the gene is escaping XCI.Panel (A) shows gene expression of 10 gametolog pairs that show substantial Y chromosome homolog contribution to gene dosage balance (ZFX/Y, USP9X/Y, DDX3X/Y, KDM5C/D, RPS4X/Y1, NLGN4X/Y, PRKX/Y, EIF1AX/Y, PDCDH11X/Y, TLX1X/Y).A sex difference was noticed in RPS4X/Y, with an increased gene dosage in female cell lines (p < 0.01).Panel (B) displays two gene pairs (TXLNG/Y, KDM6A/UTY) that display excessive Y chromosome homolog expression leading to an unbalanced gene dosage.For TXLNG/Y, this led to a male bias in gene dosage at all time points D0, D4, D9, and D37 (p < 0.05, p < 0.0001, p < 0.0001, p < 0.001, respectively).For the gene pair KDM6A/UTY, a male-biased gene dosage was found in D4, D9, and D37 (p < 0.01, p < 0.01, p < 0.001, respectively).Differences between the male and female samples were measured using repeated measurements two-way ANOVA with Šidák's multiple comparisons test.

TABLE 2
Number of differentially expressed genes between biological sexes during neuronal differentiation.

TABLE 3
Description of candidate genes: Sex-biased genes overexpressed in male (pos.Log2FC) or female cell lines (neg.Log2FC) with implications in neuron development.

TABLE 4
Grouping of X/Y homologs according their expression during neuronal differentiation and sex-bias in gene dosage.