Sex differences in gene regulatory networks during mid-gestational brain development

Neurodevelopmental disorders differ considerably between males and females, and fetal brain development is one of the most critical periods to determine risk for these disorders. Transcriptomic studies comparing male and female fetal brain have demonstrated that the highest difference in gene expression occurs in sex chromosomes, but several autossomal genes also demonstrate a slight difference that has not been yet explored. In order to investigate biological pathways underlying fetal brain sex differences, we applied medicine network principles using integrative methods such as co-expression networks (CEMiTool) and regulatory networks (netZoo). The pattern of gene expression from genes in the same pathway tend to reflect biologically relevant phenomena. In this study, network analysis of fetal brain expression reveals regulatory differences between males and females. Integrating two different bioinformatics tools, our results suggest that biological processes such as cell cycle, cell differentiation, energy metabolism and extracellular matrix organization are consistently sex-biased. MSET analysis demonstrates that these differences are relevant to neurodevelopmental disorders, including autism.


. Introduction
Sexual dimorphism in the context of brain structure, function and chemistry has been consistently reported. Total and regional volume, for example, differs between males and females, including during fetal development (Ritchie et al., 2018;Galjaard et al., 2019). Moreover, males present higher connectivity between different regions, while females demonstrate higher connectivity within regions (Gur and Gur, 2016). Several ASD, and biological processes such as cell cycle and cell adhesion (Reinius and Jazin, 2009;Shi et al., 2016). Curiously, there are evidences suggesting that genes on interacting pathways of sexbiased genes in fetal brains are more related to NDs risk than the differentially expressed genes themselves (Werling et al., 2016).
The enormous quantity of data obtained from the multiomics technologies required researchers to analyze the complex interactions between multiple molecules, giving rise to the systems biology field of study (Veenstra, 2021). Complex Network analysis is a systems biology approach to model interactions between biological entities such as co-regulated genes and its regulatory elements (Mähler et al., 2017). Coexpression networks are built based on correlation between the gene expression across samples, and the clusterization of co-expressed genes could indicate that they have related biological functions or participate in the same regulatory pathways (Langfelder and Horvath, 2008). There is evidence of differential patterns of regulation even when genes are not differentially expressed (Gaiteri et al., 2013;Lopes-Ramos et al., 2020). Within this framework, integrative methods that combine co-expression networks with other data types such as proteinprotein interactions (PPI) and regulatory mechanisms are not only more informative (Simões et al., 2015), but they also provide more specific results about the imputed biological data (Michoel et al., 2009). These integrative network approaches have been successfully applied to understand gene regulatory differences between males and females in different scenarios (Michoel et al., 2009;Lopes-Ramos et al., 2018. We hypothesized that the application of systems biology approaches to the understanding of sex differences in fetal brain development could unravel novel interactions and biological mechanisms underlying sex bias in neurodevelopmental disorders. Our goal was to find differentially regulated biological processes between males and females occurring during the midgestational period of brain development, when corticogenesis is the major histo-and neurogenic event. Moreover, we aimed to find signaling pathways and their effectors that were pertinent to these differences. In this context, we applied Complex Network tools and integrated their results to investigate autosomal gene co-expression and regulation in the brain from females and males human fetuses. We removed the genes on the sex chromosomes in order to increase the autosomal gene expression signal (Glass et al., 2014). Then, we conducted a primary analysis using the R package CEMiTool (Russo et al., 2018), which identifies co-expression modules (CM) based on gene expression correlation and reports an first overview of sex differences in module activity. Later, we applied the netZoo package to integrate PPI and motif data, creating regulatory networks for each sex and then comparing their community structure (Glass et al., 2013;Padi and Quackenbush, 2018). Next, we examine how TFs and genes are clustered based on their importance for the modularity of each created module (M), identifying the drivers of regulatory networks sex differences.

FIGURE
Workflow of the analyses that were carried out using CEMiTool and netZoo with male and female RNAseq data from fetal brain samples.
Finally, we searched for enrichment of these modules for biological functions, cell types, and neurodevelopmental and psychiatric disorders. When combined, these results show differences between sexes in fetal brain development mechanisms and hint at some of its possible effectors.

. Materials and methods
A workflow chart summarizing the conducted analyses in this study can be seen in Figure 1. . . Dataset and pre-processing RNA-seq data from 120 undissected fetal brain tissues collected from elective abortions were obtained from a public repository (O'Brien and Bray, 2018). The libraries were prepared from total RNA using the TruSEqStranded Total RNA Library Prep kit (Illumina), followed by ribosomal RNA depletion. A minimum of 100 million read pairs per sample were sequenced on Illumina HiSeq systems. Fetal sex was determined both by karyotyping and chromosomal sex genes expression. The dataset consisted of gene expression data from 70 males and 50 females with age ranging from 12 to 19 post-conception weeks (PCW), which were equally distributed for the known variables: age, RNA integrity number and read count (Supplementary Table S1).
The following analyses were conducted using the R language (v4.1.0). Counts for 29,875 genes were provided with betweensample normalization using the median of ratios method and variance-stabilizing transformation (Love et al., 2014), followed by a filter that removed low expressed genes. Additional details regarding the sample processing protocols and the RNA-seq preprocessing are published in the original article . As recommended for co-expression network construction (Langfelder and Horvath, 2008), genes with variance smaller than 0.1 were excluded. Hierarchical clustering based on the Euclidean distance of gene expression abundance demonstrated that one sample was an outlier. We removed this sample and further analyses were performed with the remaining 119 samples. Given its contribution to data variance, batch effect was corrected using LIMMA (v3.48.3) (Ritchie et al., 2015; Supplementary Figure S1). Gene annotation was converted to Gene Symbol with the package biomaRt (v2.48.3) (Smedley et al., 2009), and 18,904 genes remained.

. . Co-expression network modeling
The CEMiTool package (v1.16.0) (Russo et al., 2018) was used for the construction of co-expression networks and their clustering into co-expression modules (CMs). CEMiTool is an automated method that creates co-expression networks by intuitively selecting a soft-thresholding power based on the pairwise correlation between genes. This analysis was conducted on the full dataset, without sex stratification. Automatic filtering of genes with low variance using a model of negative binomial distribution to correct for mean-variance dependency and select variable genes and VST were turned off due to previous RNAseq data pre-processing steps. In order to evaluate the activity of each CM in males and females, samples were annotated according to sex and then submitted by CEMiTool to a gene set enrichment analysis (GSEA) (Subramanian et al., 2005). Genes from each individual CM found by CEMiTool had their expression normalized using z-scores, and these values were used to rank the expression in male and females to perform the GSEA. These z-scores are averaged and corrected by the number of genes to calculate a normalized enrichment score (NES). The result represents how the activity of each CM is within a phenotype, therefore some CMs might be more active in males or females. In order to consider a CM different between sexes, we have established a threshold based on NES distribution (absolute NES > 2).
Biological pathways were discovered applying overrepresentation analysis (ORA) on identified CMs using the R package clusterProfiler (v4.0.5) (Yu et al., 2012). Gene Ontology (v7.5.1), KEGG (v7.5.1) and Reactome (v7.4) gene set collection for enrichment were obtained from MSigDB , and the enrichment was considered significant if the Benjamini-Hochberg adjusted p < 0.0005 and the gene set included at least 10 genes. An interaction network was constructed by locating every gene in the PPI network obtained from the human dataset on StringDB (v11.5) (Szklarczyk et al., 2018). Genes with the highest number of interactions within a network are considered hubs. These genes were identified both in the co-expression network and in the PPI network by counting their number of interactions independently, and were labeled according to whether they were a hub in either network or in both.

. . Regulatory network inference
In order to build and analyze regulatory networks, three algorithms from the netZoo toolkit were used: PANDA, CONDOR and ALPACA. Gene regulatory networks were estimated using the Python (v3.7) implementation of PANDA (v0.8.1). PANDA gathers information from gene co-expression, PPI, and a prior TF motif-genes network to infer aggregate http://www.gsea-msigdb.org/gsea/msigdb/collections.jsp gene regulatory networks for each phenotype. Edge weights of PANDA networks are representative of the evidence of a regulatory interaction. PPI human dataset (v11.5) was inputted in the PANDA algorithm after genomic annotation conversion to Gene Symbol using clusterProfiler (Yu et al., 2012). The prior motif network was provided by PANDA developers and consisted of motif interactions of 730 TFs extracted from the Catalog of Inferred Sequence Binding (Weirauch et al., 2014). After removing the genes absent both in the motif dataset and in the gene expression dataset, the prior motif network contained regulatory information between 730 TFs and 15,150 genes. Only canonical TF-gene interactions were used in the analysis as recommended by Sonawane et al. (2017). . . Regulatory network comparison PANDA's male and female regulatory networks were initially submitted to the modularization algorithm CONDOR (v0.99.49) (Platig et al., 2016). CONDOR uses positive edge weights as evidence for the existence of an interaction between a TF and a gene, detecting bipartite community structures within each network using a maximum modularization approach (Barber, 2007).
Next, we used a graph-based approach called ALtered Partitions Across Community Architectures (ALPACA) (v0.99.49) (Padi and Quackenbush, 2018) that compares two networks and identifies de novo gene modules that best distinguish the networks. As described in the original article, ALPACA compares condition-specific networks to each other instead of to a random graph null model. In the comparison, one network is defined as the "baseline" network and the other is defined as the "perturbed." ALPACA defines a score denominated the "differential modularity score" that compares the density of modules in the "perturbed" network to the expected density in the matched "baseline" network. In our case, it allows us to contrast networks of male and female samples and partition the nodes into optimal differential modules (Ms). Since the 'baseline" community structure is the null model to which the "perturbed" network is compared, using the male network as baseline to compare to the female network could produce different results from the alternate comparison. For example, if the female network is used as the baseline, ALPACA compares the observed density of modules in the male network to the expected density in the female network, producing male-biased differential results. For this reason, we performed both analyses.
Differential modules were submitted to enrichment by ORA using the clusterProfiler package (Yu et al., 2012) was set to Benjamini-Hochberg adjusted p < 0.05. Represented Gene Ontology (GO) terms were selected based on their uniqueness and dispensability determined by REVIGO, a tool that calculates redundancy between a list of GO terms and identifies the most representative ones (Supek et al., 2011). Based on the enrichment p-value of each GO term and its level of informativeness, REVIGO selects GO terms that best describes the clusters of similar semantic terms, generating a reduced GO term list. For our analysis, we chose to generate small lists as a result (REVIGO reduction parameter of 0.05). Next, we excluded from the analysis GO terms from tissues or biological processes not related to our samples, as well as those with other GO terms in the same ontological tree with lower significance scores. The whole GO module enrichments and REVIGO results are attached at the Supplementary materials, and excluded REVIGO terms are marked gray.
. . Extracting di erential regulation drivers ALPACA also estimates how much a TF or a gene contributes to the differences in module density, calculating its "differential modularity score." Transcription factors with the highest score were extracted as described by Lopes-Ramos et al. (2021). Initially, scores were transformed to a log scale. The median and the interquartile range (IQR) for TFs and genes were calculated separately, given the structure of the bipartite network and number of interactions of each type of entity. TFs were considered significant if their differential modularity score deviated more than 1.5×IQR from the module median they participated in, while genes were selected when their score deviated more than 3×IQR.
. . Databases-specific enrichment analysis of ALPACA modules Genes from each differential module (M) were submitted to a modular single-set enrichment test (MSET) (Eisinger et al., 2013) using the whole network as background. To determine the over-representation of genes from each differential module in specific cell types, enrichment analysis was performed with a mid-gestational developing fetal brain single-cell dataset, fit to the age of our samples (Polioudakis et al., 2019). Associated neurodevelopmental disorders and psychiatric disorders module enrichment was conducted using curated gene sets for Alzheimer's disease (AD) (Jansen et al., 2019), ADHD (Demontis et al., 2018), major depressive disorder (MDD) (Howard et al., 2019), bipolar disorder (Stahl et al., 2019), ID (Ilyas et al., 2020), syndromic and non-syndromic ASD genes from SFARI , plus other independent ASD datasets (Sanders et al., 2015;Grove et al., 2019), schizophrenia (Lam et al., 2019), macrocephaly and microcephaly genes extracted from the testing panels by University of Chicago and Online Mendelian Inheritance in Man database , and cross-disorders gene lists (Deciphering Developmental Disorders Study, 2017;Lee et al., 2019).
MSET enrichment was also performed to assess the role of gonadal hormones in the regulation of ALPACA modules. In order to obtain sex steroid hormone target genes, gene set collections for the TFs estrogen and androgen receptors (ER and AR, respectively) were obtained from MSigDB (v7.5.1, c3.tft.v7.5.1.symbols.gmt) (Liberzon et al., 2015). All 5 motifs for the main estrogen receptors (ER_Q6, ER_Q6_01, ER_Q6_02, TGACCTY_ERR1_Q2, ERR1_Q2) were joined, resulting in 1,534 ER target genes. Androgen receptor targets consisted of 537 genes, obtained by joining 5 AR motifs (AR_01, AR_02, AR_03, AR_Q2, and AR_Q6). The MSET algorithm was employed with 10,000 permutations and the p-value was generated based on the over-representation of the input gene list in the dataset of interest when compared to the random sets generated by the permutations. This empirical p-value was adjusted using the Benjamini-Hochberg method.

. . Co-expression analysis reveals modules with sex-biased activity
Co-expression analysis with CEMiTool identified coexpression modules (CMs) of genes expressed from PCW 12 to 19 of the fetal brain development. These modules were also submitted to GSEA in order to determine sex activity and functional enrichment analysis based on Reactome, GO and KEGG.
Genes were grouped into 33 CMs based on the correlation of their normalized expression values, and 30 genes out of 18,904 were not correlated to any CM (Supplementary Figure S2). From all the CMs, a number of 13, 11 and 8 modules were significantly enriched for Reactome, GO and KEGG terms, respectively (adjusted p < 0.0005, and minimum of 10 genes) (Supplementary Table S2). Six out of 14 CMs that were considered different between sexes (absolute NES > 2) are significantly enriched for these databases (CM1, CM9, CM16, CM19, CM25, and CM31) (Supplementary Figure S3). Therefore, these following CMs will be discussed below.
Three selected modules were more active in females. CM1 is the largest module, composed of 5,306 genes. CM1 is more active in female samples (adjusted p < 0.001, NES = .

2.35) and enriched for
Reactome, GO and KEGG pathways associated with neurotransmitters, membrane receptors, ion channels, and signal transduction proteins, that act together in a complex molecular interaction network important to the synapses organization (Burke and Bender, 2019). This CM is also enriched for GO terms for cell maturation and neuron differentiation associated with synapses specification. CM16 (253 genes) is also more active in females (adjusted p < 0.001, NES = 2.88), and its genes are enriched for the citric acid cycle and respiratory electron transport pathways in both Reactome and GO databases. The last CM that is more active in females (adjusted p < 0.001, NES = 3.17), CM31 (78 genes) demonstrates significant enrichment for activation of HOX genes in hindbrain in the Reactome database, a term for which the module more active in males CM9 is also enriched. The highest difference in sex activity (adjusted p < 0.001, NES = 3.43) was observed in CM9 (654 genes), more active in males. Hub genes tend to play important roles in the network. In CM9, genes like BRCA1, CCNB2, CDCA8, and genes from the KIF family are central ( Figure 2A). Moreover, CM9 is significantly enriched for cell division processes, cell cycle repair and senescence and stress response according to all databases (see Figure 2B for Reactome results). Likewise, CM19 (229 genes) is more active in males (adjusted p < 0.001, NES = 2.83) and is consistently associated with ECM organization and ECM proteoglycans.
Altogether, CEMiTool results indicate 4 modules with sex differences in activity related to the following biological processes: (1) neuron differentiation and synapse specification; (2) mitochondria energy production; (3) cell proliferation and stress response; and (4) ECM organization. The first and second are more active in females, and the other 2 in males ( Figure 2C).

. . Regulatory networks community structure comparison
We modeled male and female PANDA networks integrating gene expression, transcription motif and PPI data. These networks indicate how strong are the evidences of regulatory interactions between TFs and genes for each sex in the form of bipartite networks. Then, we identified community structures of male and female PANDA networks with CONDOR. In order to identify transitions in their regulatory network structures, ALPACA was run twice, alternating the sex network used as baseline. Therefore, differential networks using the female network as baseline will be referred to as femalebaseline differential network, while results achieved with the male network as baseline will be referred to as male-baseline differential network. In our case, it is important to have in mind that the modularization of the female-baseline differential network result will give us, as explained on the Methods section, modules (Ms) that have a different density observed in the male network than the expected from the density in the female network. Therefore, it will give us male-biased modules and the associated results will be addressed as male-biased. Accordingly, using the male as baseline will give us differential network modules addressed as female-biased. Nevertheless, both results suggest differential patterns in gene regulation between male and female regulatory networks, but do not imply a direction in the difference like differential expression analysis or in activity like CEMiTool analysis. It also does not imply that a module is more or less regulated in any sex. ALPACA detected 7 differential modules when the female network was used as the baseline, which were determined male-biased modules (Supplementary Figure S4A). Each module was submitted to GO over-representation analysis using clusterProfiler, and 4 out of the 7 modules had significant enrichment (adjusted p < 0.05, gene set > 5) (Supplementary Table S3A). Next, we applied REVIGO to eliminate redundant GO terms, identifying a total of 118 GO terms that best described the 4 modules (Supplementary Table S3B). Male-biased M1 was enriched with GO terms associated with synapses (e.g. chemical synaptic transmission, inorganic ion transmembrane transport, and neurotransmitter transport); cell-cell adhesion, import into cell, cholesterol transport, regulation of lipid catabolic process, regulation of reactive oxygen species (ROS) metabolic process, extracellular matrix organization, cytokines, and inflammatory response; male-biased M2 with regulation of GTPase activity; male-biased M6 with different terms associated to cell cycle (e.g., mitotic cell cycle process, chromosome segregation, and cell cycle checkpoint signaling); and male-biased M7 with establishment of planar polarity, mitotic cell cycle, DNA repair, telomerase activity, mitochondrial translation, and gene expression. To visualize the results, we selected each male-biased module's top 10 REVIGO terms with the highest enrichment scores that best represented its category ( Figure 3).
Meanwhile, male-baseline differential network ALPACA analysis revealed 10 differential modules, considered femalebiased (Supplementary Figure S4B). Enrichment analysis revealed significant GO terms for 3 modules (adjusted p < 0.05, gene set > 5) (Supplementary Table S3C), and REVIGO reduced the number of terms to 84 (Supplementary Table S3D). Female-biased M1 was enriched with GO terms associated with cytokines production and immune system, reverse cholesterol transport, ROS metabolic process, metabolism (e.g., regulation of cellular ketone metabolic process, amine metabolic process, and regulation of nitric oxide metabolic process), extracellular modifications and tissue remodeling; female-biased M2 with GTPase activity and synaptic signaling; and female-biased M9 with process utilizing autophagic mechanism, and different terms associated with mitochondrial activity. The selected top 10 REVIGO terms can be seen in Figure 4.  In summary, ALPACA results point to regulatory differences between females and males in modules associated with: 1) synapse organization; 2) ECM modifications; 3) immune factors; 4) cell cycle and cell senescence; 5) metabolism and mitochondria function. It is important to note that the majority of these processes appear both in male-and female-biased modules, suggesting that they have differences in modularity comparing the observed density of nodes to the expected in both analyses. In addition, it is possible to observe that many of these biological processes were demonstrated to be enriched both in male-and femalebiased modules, indicating that there are sex-specific differences in modularity involved in these differentially regulated processes.

FIGURE
GO terms enrichment of modules from female-baseline di erential network resulted from ALPACA analysis of di erences in network topology from PANDA networks. GO terms were chosen according to their uniqueness, dispensability and significance based on REVIGO analysis. The color of the circle indicates the over-representation adjusted p-value, and the size represents the number of genes corresponding to the GO term.
. . Di erential modularity core transcription factors and genes The contribution of each TF and gene is calculated by ALPACA and defined as a score of "differential modularity." Top TFs were defined as those which scores were 1.5×IQR higher than the module median and top genes those which deviated 3×IQR.
The female-baseline differential network comprises 24 top TFs and 3 top genes. These TFs and genes are associated with the modules that have a different density in the male network than the expected based on the female network, and thereby are classified as male-biased (Table 1). Most of the top TFs are from the Ets family, such as ELK1, ERF and ERG, all located in the male-biased M7. These TFs have been implicated in biological functions like cellular proliferation, differentiation, tissue remodeling, and apoptosis (Sementchenko and Watson, 2000). The top TF was ETV1, a protein that is associated with neuron cell differentiation (Flames and Hobert, 2009;Abe et al., 2012), connection formation (Arber et al., 2000), and is regulated by the androgen receptor (Cai et al., 2007). Some top TFs from the male-biased M6 are from the CREB and ATF families, which are involved in pathways of cAMP-dependent phosphorylative activation of multiple kinases that are pertinent to neuron proliferation, survival, maturation, differentiation, and synaptic establishment (Alberini, 2009;Landeira et al., 2016). Lastly, top TFs like NFYA, NFYB, NFYC, and GABPA, in addition to the top gene SARS2 are all relevant to mitochondrial function and lipid metabolism (Yokogawa et al., 2000;Yang et al., 2014;Benatti et al., 2015).

FIGURE
Enrichment of modules from the male-baseline di erential network for GO terms selected by REVIGO. The intensity of the color of the circle indicates the over-representation adjusted p-value, and the size represents the number of genes corresponding to the GO term.
Interestingly, 13 TFs have a high score of differential modularity in both female-and male-baseline differential networks (Supplementary Table S5). Aforementioned top TFs from the Ets family are important to the female-biased M8 in a similar fashion to male-biased M7. Likewise, TFs from the CREB and ATF families are clustered in the male-biased M7, repeating their pattern in the female-biased M6. Integrating these results with CEMiTool co-expression modules, we observe that some of the shared top TFs are allocated in modules with Frontiers in Human Neuroscience frontiersin.org . /fnhum. .  Table S5). The top TF ERF belongs to the CM9, which is more active in males and related to cell cycle control, in consonance with the femalebiased M7 functional enrichment. In addition, top TFs ATF1, ATF2, and ETV1 are co-expressed in CM2, more active in males, while ATF7 and CREB5 are co-expressed in CM1 and more active in females.
. . ALPACA's di erential modules enrichment for databases of interest Using MSET, each ALPACA's differential module was tested for enrichment to genes expressed in mid-gestational cell types ( Figure 5 and Supplementary Table S6A). Since our dataset consists of non-dissected bulk RNA-seq, these results indicate that biological processes that are characteristic of certain cell types could be differentially regulated in ALPACA's differential modules.
Female-biased M2 and male-biased M3 are both enriched for genes expressed in differentiated or maturing excitatory neurons, which are generated during neurogenesis (Taverna et al., 2014). Female-biased M2 is also enriched for GO terms such as synaptic signaling and cell part morphogenesis, which are representative of neuronal activity in differentiated neurons and of dendritic and axonal growth in maturing neurons.
Another shared cell type enrichment pattern occurs in male-biased M7 and female-biased M9 for cycling progenitors (S phase) (FE = 1.8; adjusted p = 0.003, and FE = 1.58; adjusted p = 0.012, respectively) and intermediate progenitor (FE = 2.08; adjusted p < 0.001, and FE = 3.14; adjusted p < Frontiers in Human Neuroscience frontiersin.org . /fnhum. .  0.001, respectively). These modules are associated with energy production and cell cycle regulation. The cell cycle related male-biased M6 is enriched for genes expressed in outer radial glia (FE = 1.72; adjusted p = 0.008) and ventricular radial glia (FE = 2.37; adjusted p < 0.001), which are proliferative cells, as well as with cycling progenitors from G2/M (FE = 3.51; adjusted p < 0.001) and S phase (FE = 3.06; adjusted p < 0.001) and intermediate progenitors (FE = 2.22; adjusted p = 0.003). This enrichment pattern does not repeat in any female-biased module. Overall, cell type gene expression enrichment correlates with biological processes enriched within the same ALPACA module.
Lastly, we analyzed our ALPACA differential modules for enrichment of gene targets of sex steroid hormone receptors (Supplementary Table S6C). Male-biased modules M1 and M6 were enriched for estrogen receptor (ER) targets (FE = 1.27; adjusted p < 0.001, and FE = 1.46; adjusted p = 0.01, respectively). For female-biased modules, only M2 was enriched for ER targets (FE = 1.29; adjusted p < 0.001). Curiously, none of the modules were enriched for androgen receptors.

. Discussion
Our study used integrative approaches that account for biological complexity and takes advantage of its interacting entities to better understand sex differences in brain autosomal gene expression regulation. In this context, we applied network tools that take advantage of community information, since methods that identify changes in groups of nodes rather than individual edges or nodes have been demonstrated to find more robust differences between distinct phenotypes (Padi and Quackenbush, 2018). CEMiTool clusters genes based on the correlation of their expression, and could identify co-expression modules differentially activated in two conditions. Co-expressed genes often have related biological functions or participate in the same regulatory pathways (Langfelder and Horvath, 2008). Co-expression networks are built as a consequence of gene expression regulation, however they do not convey direct information about regulatory interactions. Therefore, we combined co-expression networks with regulatory networks inferred by the netZoo tools in order to extract this information. These complementary analyses also increase the reliability of our results that are achieved by the two independent tools, and thereby we expect similar biological processes to be pointed out as differentially co-expressed and regulated between sexes in both analyzes. Moreover, we observe that these differences in the regulation of biological processes are a consequence of the sex differences in modularity, and the addition of the CEMiTool results shows in which sex this biological process is more active.
In a recent review, Kostović et al. (2019) divided the fetal and perinatal developmental stages characterized by the occurrence and intensity of major histo-and neurogenetic events, originating transient patterns of microstructural features that will later differentiate into the multiple areas of the brain. The early fetal period (8-12.5 PCW) is characterized by three histo-and neurogenetic events: proliferation, migration and molecular specification, culminating in the establishment of the cortical plate. From 13 to 15 PCW, these events continue while cells begin to change aggregation patterns and migrate tangentially, causing the dispersion of cortical plate cells, reducing its density, and ultimately forming the cortical subplate. Through the midfetal period (15-23 PCW) four events are dominant: neuronal aggregation and cytoarchitecture changes, axonal outgrowth and ingrowth, dendritic differentiation, and molecular specification.  Moreover, neuronal differentiation, dendritic development, astroglia differentiation, gliogenesis and synaptic activity are also observed during the midfetal period.
In this context, this study contributes to the understanding of sex bias in brain development by capturing differences in pathways and biological processes in accordance Frontiers in Human Neuroscience frontiersin.org . /fnhum. . with the brain development events that have been described for the fetal period of our samples, which range from 12 to 19 PCW of age. Collectively, CEMiTools, and netZoo results showed sex differences mainly in processes associated with progenitor proliferation and survival as well as cell fate decisions pertinent to cell delamination and differentiation, cell migration and sinaptogenesis. Sex differences in cell proliferation is a consistent result along the conducted analyses. ALPACA's male-biased M6 is associated with cell cycle regulation while no female-biased module is, indicating that the regulatory patterns of this function in females are organized differently in males. This same module is enriched for multiple proliferative cell types as demonstrated by MSET analysis using mid-gestational brain gene sets, corroborating that these regulatory interactions are relevant to the neural progenitor activity and contribute to sex differences. Complementarily, CEMiTool demonstrates that cell cycle associated genes present the highest difference in activity between sexes, with higher mean expression in males. It has been demonstrated that prenatal brain growth in males is greater than in females from very early periods, and fetal head circumference measures are consistent with this concept (Gutiérrez-Adán et al., 2006;Broere-Brown et al., 2016). This higher proliferative profile has been well observed in cancer, resulting in higher incidence and mortality in males (Rubin et al., 2020). Abnormal neuronal proliferation has also been associated with NDs more common in males, such as ASD and ID (Ernst, 2016;Marchetto et al., 2016). Accordingly, male-biased M6 genes are enriched for two ASD gene sets, suggesting that this differential regulation in cell proliferation could be related to the male-bias in ASD.
We have also observed cohesive sex differences in the regulation of metabolism. Our discoveries show that genes related to the citric acid cycle, oxidative phosphorylation, and ATP synthesis are differentially regulated, while also more active in females. Metabolic processes are intrinsically involved in brain development because they regulate energy homeostasis, mediating nutrient availability, stress response and hormonal stimuli, and cell growth and division (Fritz and Fajas, 2010). The proliferative state of cells during neurogenesis is correlated with a decreased oxidative phosphorylation metabolism and an increase in the glycolytic activity, which is less efficient but produces ATP faster allowing for the augmented biosynthesis of macromolecules (Ježek et al., 2010). The glycolytic pathway includes the oxidation of fatty acids, which is a signal to induce cell stemness (Knobloch, 2017). It has been demonstrated that in normal conditions males favor the glycolytic pathway to generate ATP quickly, increasing their glucose and pyruvate uptake and producing more lactate, while females favor the pentose phosphate pathway (Rubin et al., 2020). Previous findings demonstrated that given an abundance of nutrients, normal proliferating cells can use glucose and glutamine as substrates to produce ATP and maximize their growth via PI3K/mTOR signaling pathways (Ray et al., 1995;Heiden et al., 2009). Research with murine hearts demonstrate that the PI3K/mTOR pathways are sexually dimorphic (Gürgen et al., 2013). TP53 gene regulation is more active in males according to our results. It is intricately associated with our top TFs in protein-protein interaction networks, and together with the top TFs from the MYC family it plays an important role in the glycolytic pathway (Whan Kim et al., 2004). MYC genes are known to stimulate glucose uptake, NADPH biosynthesis for redox homeostasis and glycolysis (Dang, 2010). In fact, the motif for MYC-associated zinc finger protein has been shown to be extremely differentially targeted between sexes across multiple tissues even when there are no differences in targeted gene expression (Lopes-Ramos et al., 2020). Therefore, we provide additional evidence to the theory that males naturally present a more proliferative metabolic profile, which seems to be partially defined by the activity of TFs from the MYC family.
Our results also point to a sex biased regulation of lipid catabolic processes, reverse cholesterol transport and amine and steroid metabolism. These lipid metabolic processes are essential for neural progenitor survival, proliferation and differentiation, in addition to ROS homeostasis, and mitochondrial biogenesis and function (Fame and Lehtinen, 2021). Cholesterol, for instance, is one of the most abundant lipids in neuron membranes and in the myelin sheaths, being essential for proper synaptogenesis, brain development and all sex steroid hormones biosynthesis (Orth and Bellosta, 2012). Consequently, the loss of the capacity to produce cholesterol in murine neural stem cells induces apoptosis, halting proliferation (Saito et al., 2009). Neural stem cells metabolism is primarily glycolytic, shifting to oxidative as they mature, and higher lipogenesis activity is essential to neuronal differentiation since lipids are necessary for the constitution of the membrane of neurons, their projections and vesicles (Knobloch, 2017). Lipogenesis and nucleotide synthesis are increased with the redirection of the glucose carbon flux to anabolic pathways caused by aerobic glycolysis, a process regulated by one of the top TFs we encountered, MLXIP (Fritz and Fajas, 2010). Indeed, females prioritize lipidic pathways, synthesizing and consuming more lipids (Ockner et al., 1979). These sex differences in cell metabolic profiles could be related to the higher proliferation of neural stem cells in males, while females shift earlier to an oxidative profile that differentiates progenitor cells into mature neurons. It is also notable that sex steroid hormones are intricatedly involved in metabolic differences through cAMP signaling pathways (Mayes and Watson, 2004), and their role in neuronal proliferation and differentiation has been extensively described (Karaismailoglu and Erdem, 2013;Miranda and Sousa, 2018).
In addition to the aforementioned metabolic functions of mitochondria in energy metabolism, apoptosis regulation and steroid hormone production, ALPACA also indicates regulatory differences in terms related to mitochondrial overall activity (mitochondrial gene expression, mRNA metabolic process, and mitochondrial translation) represented mainly in femalebiased ALPACA module M9. Indeed, females have higher mitochondrial activity than males in blood and muscle cells (Cardinale et al., 2018;Silaidos et al., 2018), and CEMiTool reveals that related genes are more active in female fetal brain samples as well. Studies demonstrated that mitochondria .
/fnhum. . from females produce more antioxidant proteins, therefore accumulating less ROS (Borrás et al., 2003). In consonance, glucocorticoid mediated ROS increase during gestational stress is higher in males, and sex-specific placental adaptation makes male fetuses keep growing even in adverse situations (Stark et al., 2011). ROS can have multiple sources in cells, including the activity of endothelial nitric oxide synthase, xanthine oxidase, membrane bound NADPH-oxidase and from the electron transport chain (Nayernia et al., 2014), which appeared in our results as differentially regulated. Mitochondrial functions have been associated with the coordination of cell proliferation, differentiation, migration and maturation in corticogenesis (Fame and Lehtinen, 2021), and ROS dynamics are essential to regulate cell fate (Inoue et al., 2016). Indeed, altered metabolic activities of ROS during neurogenesis are associated with neurodevelopmental disorders such as schizophrenia (Paulsen et al., 2012). There is evidence that aerobic glycolysis reduces ROS levels, and that their interplay is related to cell fate decision with the participation of MYC proteins (Rodic and Vincent, 2017). These intricate interactions are signaling factors for proliferation, differentiation, apoptosis, and autophagy when ROS is overabundant (Reczek and Chandel, 2017). Autophagic mechanisms appear as differentially regulated, and associated with female-biased modules. Male cells are normaly more autophagic than female cells during the prenatal period (Addis et al., 2014), and autophagy has been suggested to be regulated by sex dimorphic pathways such as PI3K/mTOR and MAPK (Shanware et al., 2013). At least in parts, these differences are attributed to differential regulation of the insulingrowth factor 1 (IGF-1), which has been demonstrated to be regulated in sexually dimorphic fashion and increases glucose uptake (Clifton et al., 2010). Our results suggest that this balance is differentially regulated between sexes, as seen in the enriched functions of female-biased M7 and male-biased M9. These modules are also enriched for several types of progenitor cell corroborating that their genes are relevant to proliferative functions. Altogether, these results could indicate that this complex metabolic conjecture is underlying higher male susceptibility to NDs in complicated gestations.
ALPACA's female-biased module M2 is associated with cell part morphogenesis and synaptic signaling, which are functions related to tissue remodeling. During the morphogenesis of the cerebral cortex, multiple cytoarchitecture transient structures are shaped and subsequently rearranged due to neural progenitor cells proliferation, migration and differentiation, culminating in tissue histogenesis (Kostović et al., 2019). As reviewed by Sokpor et al. (2022), this process requires delamination and depends on the dismantling of adherens junctions, modification of cell polarity, cytoskeleton remodeling, and reposition of organelles and signaling molecules after mitosis. Therefore, changes in ECM dynamics are indispensable to the regulation of neurogenesis, and improper control of these events affect brain development (Long and Huttner, 2022). In fact, delamination is suggested to be the most important mechanism to allow for cell migration and differentiation in the corticogenesis (Sokpor et al., 2022). ECM components are also relevant to neurogenesis due to their pro-proliferative signals through cell-cell integrins interaction (Kalebic and Huttner, 2020). Interestingly, CEMiTool module CM19 is more active in males and associated with ECM. Male-biased gene expression in developing human brain samples has been linked to ECM genes, which have been implicated in the etiology of autism (Ziats and Rennert, 2013). Both CEMiTool and ALPACA results corroborates sex differences in functions pertinent to these ECM related processes, such as ECM organization, regulation of cell-cell adhesion, establishment of planar polarity and tissue remodeling. Sex differences encountered in synaptic signaling and immune response might reflect the relevance of neurotransmitters and cytokines as signaling molecules that control processes during brain development (Deverman and Patterson, 2009;Ojeda and Ávila, 2019). These processes are enriched in both female-and male-biased ALPACA modules and are more active in females in the CEMiTool module CM1. Neurotransmitters are essential for cell migration due to its chemotropic effects on migrating neurons (Behar et al., 1996), and there is abundant evidence that neurotransmitters and their receptors modulate gene expression through calcium and CREB signaling pathways impacting on neurogenesis, neuronal differentiation and programmed cell death (Jagasia et al., 2009;Ascenzi and Bony, 2017). Studies with animal models have demonstrated neurochemical differences between males and females, mostly due to sex hormones action (Uhl et al., 2022), and our results corroborate the hypothesis that during corticogenesis estrogen is a relevant signaling molecule (Denley et al., 2018). In a similar fashion, cytokines are critical signaling factors for the renewal of neural stem cells, regulation of differentiation between neurons and glial cells, determination of neuronal identity, and chemotropism, modifying cell migration and axon guidance (Deverman and Patterson, 2009). One of the differentially regulated cytokines between males and females we identified, interleukin-6 is temporally and spatially controlled in order to guarantee proper brain development, participating in multiple pathways of neuronal survival and differentiation (Kummer et al., 2021). For instance, this cytokine can promote ambiguous results: activating STAT-3 pathway, it induces gliogenesis, and through the MAPK/CREB pathway it can stimulate neurogenesis (Erta et al., 2012). Interestingly, women carrying female fetuses produce more IL-6, TNF-α and IL-1β (Mitchell et al., 2017).
A previous study with the same dataset identified differences in sex expression , and demonstrated that genes more expressed in males are enriched in neural progenitor cells marker genes, while genes more expressed in females are enriched in Cajal-Retzius cells and glia. We have identified enrichment for neural progenitors and glial cells in sex-biased Frontiers in Human Neuroscience frontiersin.org . /fnhum. .
modules, but Cajal-Retzius cells markers were not included in the single-cell database. However, we would like to highlight that unlike in differential expression analysis, sex-biased modules only indicate that regulatory patterns are different and do not convey a direction of the differences. For this reason, we can not assume that any sex is enriched for any cell type, but we can infer that sex differences in regulatory patterns are relevant to cell type processes. Moreover, the differences between the two studies could also be attributed to methodological differences and the usage of distinct databases. The same limitation applies to the enrichments found for genes in disorders databases. A transcriptomic study identified that genes that are more expressed in fetal male brains are enriched for disorders like schizophrenia, Alzheimer's disease and autism (Shi et al., 2016). In our analyses, we found sex differences in regulatory patterns of modules enriched for intellectual disability and mainly ASD. Once again, our methods do not imply that the enrichment is specific to any sex. Besides, our dataset contains a larger number of samples and a narrower and more specific time window of the fetal period, which could explain why some enrichments were not reproduced.
Altogether, our study provides evidence that autosomal genes related to biological functions pertinent to neurodevelopment are differentially regulated between sexes. The mechanisms that are possibly responsible for these differences are the gonadal hormones, the sex chromosomes, and the sex-specific epigenetic events (McCarthy, 2016;Maschietto et al., 2017;Tahira et al., 2018). Our findings reinforce that CREB and MYC participate directly in the pathways that differentiate males from females during development (Auger, 2003), and that gene expression regulation of genes involved in metabolism, proliferation, and delamination contribute to neurodevelopmental differences between males and females. Those processes are associated in a complex and intricate mechanism relevant to neurodevelopmental disorders sex bias, suggesting that multiple factors with evidences of sex differentiation interact and regulate development. Important limitations from this study are the usage of bulk tissue samples and a temporal interval, which narrow the spatial and temporal understanding of the development dynamics, and the fact that our results are inferences based on bioinformatics tools and molecular biology assumptions, requiring further exploration of these regulatory interactions with practical experiments. Nevertheless, we believe our study lays a path for the connection of sex differences in multiple pathways with interactions that had not been established yet. These findings could help to understand how these differences in brain development impact on sex-biased neuropsychiatric disorders and provide targets for sex-specific terapeutic interventions.

Data availability statement
Publicly available datasets were analyzed in this study. This data can be found here: https://doi.org/10.6084/m9.figshare. 6881825.v1. No new custom software or code has been used in this paper. Codes (R and Python programming language) are all based on pre-existing, peer-reviewed published softwares for RNA-Seq and network analyses, as described above. These codes are available at https://github.com/victortoledo/frontiers_ sexdiff_fetal_brain.git, along with some of the multiple data that were integrated, in order to ensure clarity and reproducibility.

Author contributions
VT, AT, and HB designed the study and interpreted the results. VT, AF, AB, and AT conducted the analyses. VT and HB drafted the manuscript. All authors contributed with the further writing of the manuscript and approved the final manuscript.

Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnhum. 2022.955607/full#supplementary-material SUPPLEMENTARY FIGURE S Singular Value Decomposition Analysis (SVD) of top principal components correlated to the samples covariates before and after batch correction using LIMMA, demonstrating the relevance of the batch e ect for data variance.

SUPPLEMENTARY FIGURE S
CEMiTool co-expression modules and their sizes. Not correlated genes are represented by the bar at the bottom.

SUPPLEMENTARY FIGURE S
Sex activity of genes from each of the modules of co-expression detected by CEMiTool. Modules with a significant enrichment (adjusted p < . ) of gene expression of either sexes calculated by GSEA and an absolute NES > are highlighted and are considered relevant. In red, modules that are more active in females; in blue, modules more active in males.

SUPPLEMENTARY FIGURE S
ALPACA's female-and male-baseline di erential module sizes and its composition by TFs (cyan) and genes (purple).

SUPPLEMENTARY TABLE S
Covariates and principal components associated with samples as described by O'Brien and Bray ( ).

SUPPLEMENTARY TABLE S
(A) CEMiTool module enrichment by over-representation analysis using Reactome pathways as background. (B) CEMiTool module enrichment by over-representation analysis using GO terms as background. (C) CEMiTool module enrichment by over-representation analysis using KEGG pathways as background.

SUPPLEMENTARY TABLE S
(A) ALPACA's male-biased modules enrichment by over-representation analysis of GO terms. (B) GO terms selected with REVIGO using genes from ALPACA's male-biased modules. (C) ALPACA's female-biased modules enrichment by over-representation analysis of GO terms. (D) GO terms selected with REVIGO using genes from ALPACA's female-biased modules.

SUPPLEMENTARY TABLE S
Female-biased TFs and genes with the highest contribution to di erential modularity.

SUPPLEMENTARY TABLE S
Shared top TFs between male-and female-biased regulation modules and their membership in CEMiTool co-expression modules.

SUPPLEMENTARY TABLE S
(A) MSET enrichment results for ALPACA's di erential modules using mid-gestational cell type gene lists as background. (B) MSET enrichment results for each ALPACA's di erential modules using neuropsychiatric disorders associated genes as background. (C) MSET enrichment results for ALPACA's di erential modules using lists of gene targets of sex steroid hormones receptors as background.