Sex differences in tissue gene regulation — contributions of three classes of sex-biasing factors

Sex differences in tissue gene regulation—contributions of three classes of sex-biasing factors Montgomery Blencowe1,2, Yuichiro Itoh1.,3,7, Xuqi Chen1,3, Yutian Zhao1, Yanjie Han1, Benjamin Shou, Rebecca McClusky1,3, Karen Reue4, Arthur P. Arnold 1,2,3*, Xia Yang1,2,4,6* 1Department of Integrative Biology & Physiology, University of California, Los Angeles, California USA 2Interdepartmental Program of Molecular, Cellular and Integrative Physiology, University of California, Los Angeles, California USA 3Laboratory of Neuroendocrinology of the Brain Research Institute, University of California, Los Angeles, California USA 4Department of Human Genetics and Molecular Biology Institute, David Geffen School of Medicine, University of California, Los Angeles, California USA 5Molecular Biology Institute, University of California, Los Angeles, California USA 6Institute for Quantitative and Computational Biosciences, University of California, Los Angeles, California USA 7Department of Neurology, David Geffen School of Medicine, University of California, Los Angeles, California USA Corresponding Authors Arthur Arnold, Ph.D. arnold@ucla.edu Department of Integrative Biology & Physiology, UCLA 610 Charles E Young Drive South Los Angeles CA 90095-7239 Xia Yang, Ph.D. xyang123@ucla.edu Department of Integrative Biology & Physiology, UCLA 610 Charles E Young Drive South Los Angeles CA 90095-7239


INTRODUCTION
The incidence, progression, clinical manifestation, and genetic risks of many diseases, such as metabolic disorders including obesity, non-alcoholic fatty liver disease, and diabetes, differ between females and males (Arnold 2010;Clayton and Collins 2014;Rask-Andersen et al. 2019), which indicates that one sex may have endogenous protective or risk factors that could become targets for therapeutic interventions. Current sexual differentiation theory suggests that three major classes of factors cause sex differences (Arnold 2009;Arnold 2012;Arnold et al. 2013;Schaafsma and Pfaff 2014;Arnold 2017). First, some sex differences are caused by different circulating levels of ovarian and testicular hormones, known as "activational effects". These differences are reversible because they are eliminated by gonadectomy of adults. Second, certain sex differences persist after gonadectomy in adulthood and represent the effects of permanent or differentiating effects of gonadal hormones, known as "organizational effects," that form during development. A third class of sex differences are caused by the inequality of action of genes on the X and Y chromosomes in male (XY) and female (XX) cells, and are called "sex chromosome effects".
To date, few studies have systematically evaluated the relative size and importance of each of these three classes of factors acting on phenotypic or global gene regulation systems (Arnold 2019). For instance, the activational effects of hormones have been established as a significant contributor to sexual dimorphism in metabolic diseases, with additional evidence pointing to sex chromosome effects on obesity and lipid metabolism (Chen et al. 2012;Link et al. 2015;Link et al. 2020). Previous studies have also emphasized the importance of organizational or activational hormone effects that contribute to sex differences in gene expression in the liver (Mode and Gustafsson 2006;van Nas et al. 2009;Waxman and Holloway 2009;Sugathan and Waxman 2013;Zheng et al. 2018). However, the tissue-specific contributions and the interactions of activational, organizational, and sex chromosome effects on gene regulation are poorly investigated.
Here we conduct a systematic investigation to tease apart the three sex-biasing factors in gene regulation (Figure 1). We used the Four Core Genotypes (FCG) mouse model, in which the type of gonad (ovary or testis) is independent of sex chromosome complement (XX or XY) (De Vries et al. 2002;Burgoyne and Arnold 2016). The model is a 2x2 comparison of the effects of sex chromosome complement (XX vs. XY) and type of gonad (ovaries Sry which is usually located on the Y chromosome was deleted (a spontaneous deletion) and inserted as a transgene onto to chromosome 3. B. Production of gonadal male with 4 possible offspring types. A gonadal male mouse with XY -Sry + produces 4 gametes which can result in the four core genotypes (FCG). C. FCG. Four types of mouse offspring were produced (two gonadal males and two gonadal females): XY -Sry + (XYM), XXSry + (XXM), XX (XXF), XY -(XYF). D. Three varieties of hormone replacement for each group. Each of the four core genotypes underwent a gonadectomy (GDX) at day 75 and were implanted with a capsule that contained either estradiol, testosterone or blank. E. Tissue removal. Two metabolically relevant tissues were selected for removal, the liver and the inguinal adipose tissue, then RNA isolation was completed. F. Gene expression data generation and quality control. Using an Illumina microarray, we measured the transcriptome and then carried out a principal component analysis (PCA) to identify outliers and global patterns. G. Bioinformatics analyses. Differentially expressed genes (DEGs) influenced by individual sexbiasing factors were identified using 3-way ANOVA (chromosomal, gonadal and hormonal effects), 2-way ANOVA (gonadal and chromosomal effects under each hormone condition), and a 1-way ANOVA (estradiol and testosterone treatment effects in individual genotypes). Gene coexpression networks were constructed using MEGENA and differential coexpression modules (DMs) affected by individual sex-basing factors were identified using 3-way, 2-way, and 1-way ANOVAs. DEGs and DMs were analyzed for enrichment of functional categories or biological pathways. The relevance of the DEGs to human disease relevance was assessed via integration with human genome-wide association studies (GWAS) for >70 diseases using the Marker Set Enrichment Analysis (MSEA). Transcription factor analysis and gene regulatory network analysis were additionally conducted on the DEGs derived from the one-way ANOVA.
Using the FCG model, we assessed the role of the three sex-biasing factors on gene expression, molecular pathways, and gene network organization in the liver and adipose tissue, which are central tissues for metabolic and endocrine homeostasis, with adipose tissue additionally contributing to immune functions. We found that all three factors influence gene regulation involved in immune response, cell development, and metabolic processes, but to different degrees and with tissue and molecular specificity. We further integrated sex-biased genes and networks influenced by each sex-biasing factor with diverse diseases in human populations to predict the phenotypic consequences of each sex-biasing factor. This analysis revealed extensive associations of the sex-biased genes with coronary artery disease (CAD), type 2 diabetes (T2D), autoimmune disorders, and a variety of lipid traits. Our study provides a roadmap of the relationship between hormonal and genetic origins of biological sex and prevalent human diseases.

Overall study design
We used FCG mice on a C57BL/6J B6 background. In FCG mice, the Y chromosome (from strain 129) has sustained a spontaneous deletion of Sry, and a Sry transgene is inserted onto chromosome 3 (Burgoyne and Arnold 2016). Here, "male" (M) refers to a mouse with testes, and "female" (F) refers to a mouse with ovaries. FCG mice include XX males (XXM) and females (XXF), and XY males (XYM) and females (XYF; Figure 1). A total of 60 FCG mice were gonadectomized (GDX) at 75 days of age and implanted immediately with medical grade Silastic capsules containing Silastic adhesive only (blank control; B) or testosterone (T) or estradiol (E). This study design produced 12 groups (Figure 1), with 4 groups of FCG mice (XXM, XXF, XYM, XYF) and each group subdivided into B or T or E based on hormonal treatment: XXM_B, XXM_T, XXM_E, XYM_B, XYM_T, XYM_E, XXF_B, XXF_T, XXF_E, XYF_B, XYF_T, XYF_E (4-6/per genotype/per treatment). Liver and inguinal adipose tissues were collected 3 weeks later for transcriptome analysis. The design allowed detection of differences caused by three factors contributing to sex differences in traits. (1) "Sex chromosome effects" were evaluated by comparing XX and XY groups. (2) "Gonadal sex effects" were determined by comparing mice born with ovaries vs. testes. Since mice were analyzed as adults after removal of gonads, the gonadal sex effects represent organizational (long-lasting) effects of gonadal hormones, such as those occurring prenatally, postnatally, or during puberty. This group also includes effects of the Sry gene, which is present in all mice with testes and absent in those with ovaries. Any direct effects of Sry on non-gonadal target tissues would be grouped with effects of gonadal sex. (3) "Hormone treatment effects" refers to the effects of circulating gonadal hormones (activational effects) and were evaluated by comparing E vs. B groups for estradiol effects, and T vs. B groups for testosterone effects. We then asked which individual genes in liver and adipose tissues were affected by adult hormone level, gonadal sex, and sex chromosome complement, as well as interactions between these factors, using a 3-Way ANOVA (3WA). Tens to thousands of differentially expressed genes (DEGs) were identified in liver (Figure 2, Table 1) and adipose tissue (Figure 2, Table 2). In both tissues, hormonal treatments affected the largest numbers of genes, followed by fewer genes that were responsive to gonadal/organizational effects or sex chromosome complement. Testosterone treatment in the liver induced the largest number of DEGs (Figure 2A), whereas in inguinal adipose tissue estradiol treatment affected the greatest number of DEGs ( Figure 2D). These results support tissue-specific sensitivity to different hormones.

Figure 2. Bar graphs (A-F) and heatmaps (G-H)
representing the number of DEGs for each sex-biasing factor and differential co-expression modules from a 3-way, 2-way, and 1-way ANOVA, respectively. Each bar graph represents the number of DEGs based on each specific statistical analysis at FDR<0.05. A, D represent results from 3-way ANOVAs run separately in testosterone vs. blank groups, and estradiol vs. blank groups to examine hormone, gonad, and sex chromosome effects as well as interaction terms. B and E represent results from 2-way ANOVAs with factors of gonadal sex and sex chromosomes as well as interaction terms, run separately on data from testosterone (T), estradiol (E), and blank (B) treatment groups. C and F represent results from 1-way ANOVAs testing effects of hormonal treatments (vs. Blank) in each of the four genotypes for liver and inguinal adipose tissue. In A and D, bars are colored according to the hormonal treatment variable -data from estradiol vs. blank (pink), or data from testosterone vs. blank (blue). In B and E, colors represent the hormonal treatment condition (testosterone groups blue, estradiol groups pink, and blank groups white) for each of the separate 2-way ANOVAs. In C and F, colors show effects of testosterone vs blank (blue) or estradiol vs blank (pink) in each of the four genotypes. Horm = Hormone, Chr = Sex Chromosome, M = Testes/Sry present, F = Ovaries present, no Sry. G represents the heatmap for liver. H represents the heatmap for adipose tissue. Each heatmap shows results from 1-way ANOVA, 2-way ANOVA and 3-way ANOVA for each test factor hormone (H), chromosome (C), and gonad (G) when treated with testosterone (T), estradiol (E) and blank (B). Interaction terms among H, C, and G were also tested. For instance, C:G:H indicates the interaction term among the 3 factors in 3-way ANOVA. Each module was annotated with canonical pathways from GO and KEGG. When a module is not annotated with a pathway name, it did not show significant enrichment for genes in any given pathway tested. Colors correspond to the statistical significance of the effects of sex factors on modules in the form of -log10(FDR).  Next, we explored the sex chromosome and gonadal effects within each hormonal treatment group using a 2-way ANOVA (2WA). In the liver, the organizational effects of gonad type were strongest in gonadectomized mice without hormone replacement (blank group) ( Figure 2B). By contrast, in adipose tissue the gonadal sex effect was most prominent in the estradiol treated groups (Figure 2E), suggesting that estradiol levels augment the enduring differential effects of gonads on the adipose tissue transcriptome. Sex chromosome effects were limited regardless of hormonal treatment status, but captured specific genes important in inflammatory signaling and metabolic pathways.
Lastly, we examined the effects of testosterone and estradiol within individual genotypes using a 1-way ANOVA (1WA) followed by post-hoc analysis. More liver genes were affected by testosterone than by estradiol regardless of genotype, although XYM liver appeared to be less responsive to testosterone than liver from other genotypes ( Figure 2C). By contrast, in adipose tissue, estradiol affected more DEGs in XX genotypes (XXM and XXF) than in XY genotypes (XYM and XYF), whereas testosterone had minimal impact on adipose tissue gene expression in all four genotypes ( Figure 2F). These results further support tissue-specific effects of estradiol in adipose tissue and testosterone in liver, and indicate that activational effects of hormones also depend on sex chromosome complement and hormonal history (gonadal sex) of the animal. These interactions among the sexbiasing factors were further supported by the numerous DEGs with significant effects from the interaction terms in the ANOVA analyses (Supplementary Table 1). For instance, 31 DEGs were affected by interactions between estradiol and gonad type in adipose tissue, and 2-10 DEGs showed interaction effects between pairs of sex-basing factors in liver (FDR<0.05). In particular, Mpp5 (Membrane Palmitoylated Protein 5, important in cell junction organization and PI3K-AKT signaling) was found to be affected by hormone, gonadal type, and sex chromosome as well as interactions between any two factors and the interaction among all three factors in adipose.

Genes and pathways affected by hormonal treatment
In the liver, the 3WA analysis showed that testosterone treatment induced the greatest number of DEGs with 1378 compared to 333 DEGs from estradiol treatment ( In contrast to liver, we found that the effect of estradiol treatment was more profound (2029 DEGs Overall, both estradiol and testosterone affected genes involved in metabolism, development, and immune function. However, estradiol primarily affected these processes in the adipose tissue, whereas testosterone exhibited influence in the liver. Additionally, the activational effects of hormones depended on the chromosomal and gonadal sex contexts as reflected in the different DEGs between genotypes.

Genes and pathways affected by gonadal sex
In the liver, 3WA analyses revealed 93 DEGs influenced by gonadal sex when testosterone and blank treatment groups were considered, and 209 DEGs in the analysis of estradiol and blank groups ( Table 1). These genes were enriched for immune/defense response and lipid metabolism pathways. To further tease apart the effect of gonadal sex in individual hormonal treatment groups, we conducted a 2WA and found that gonadal sex has the strongest influence on inflammatory and metabolism genes in the absence of hormones (blank group; 115 DEGs), whereas the gonadal sex effects were reduced by estradiol treatment (53 DEGs) and minimized in the testosterone group (9 DEGs; Table 1; Figure 2B).
For the inguinal adipose tissue, gonadal sex showed strong effects on gene expression, with more than twice as many DEGs as in liver tissue in the 3WA analysis (Figure 2A vs. 2D). Genes involved in developmental processes were enriched in these DEGs. Further dissection of the gonadal sex effect in individual hormonal treatment groups in a 2WA analysis showed that the effects of gonadal sex were strongest in the estradiol treatment group (400 DEGs), followed by testosterone group (161 DEGs), and lastly by the blank group (70 DEGs; Figure   2E). In addition to affecting development related genes, gonadal sex also affected genes involved in arginine and proline metabolism in the estradiol treatment group, and cancer-related genes in the testosterone treatment group.
These results support the importance of gonadal sex in regulating development processes in both tissues.
However, in the liver, hormonal treatments minimized the effects of gonadal regulation of gene expression, whereas in the adipose tissue, hormones amplified the gonadal influence on gene expression to regulate additional metabolism genes (by estradiol) or cancer-related genes (by testosterone). In both tissues, the gonadal sex effect was more prominent in the estradiol-treated group than in the testosterone-treated group ( Figure 2B vs. 2E).

Genes and pathways affected by sex chromosome complement
In both the 3WA and 2WA analyses, ten or fewer genes were found to be significantly affected by sex chromosome complement at FDR < 0.05 in the liver (Table 1; Figure 2C) and 10-22 DEGs were influenced by sex chromosomes in the adipose tissue (Table 2; Figure 2F). Not surprisingly, these genes were mainly sex chromosome genes known to exhibit sex differences in gene expression levels, including Xist, Ddx3y, Kdm6a, Hccs, Cited1, Tlr7 and Eif2s3x/y (Chen et al. 2012;Berletch et al. 2015;Golden et al. 2019;Itoh et al. 2019). However, autosomal genes were also influenced by sex chromosome type in both liver (e.g., Ntrk2 and H2-dmb1) and adipose tissue (e.g., Mpp5, Rbm35a, and Dnaic1). Overall, the genes influenced by sex chromosome complement are involved in diverse processes including inflammation/immune response (Tlr7, H2-dmb1, Cited1), GPCR signaling (Rbm35a), metabolism (Hccs), and cell junction organization (Mpp5).

Coexpression modules affected by each sex-biasing factor
The above DEG analyses focused on genes that were individually influenced by sex-biasing factors. Sets of genes that are highly coregulated or co-expressed can offer complementary information on coordinated gene regulation by sex-biasing factors that might be missed by the DEG-based analyses. To this end, we constructed gene coexpression networks for each tissue using MEGENA and identified 326 liver and 131 adipose coexpression modules. The first PCs of the coexpression modules, each composed of coregulated genes, were assessed for influence by sex chromosome, gonadal sex, and hormonal treatment factors using 3-, 2-, and 1-WAs.
In the heatmaps summarizing the significance of impact of each sex-biasing factor on individual modules ( Figure 2G, 2H), we confirmed the large effect of hormonal treatment in regulating modules enriched for diverse biological pathways. In the liver, testosterone affected modules involved in metabolism (RNA, lipid, protein), development, protein assembly, chemical response, immune system (inflammation, adaptive immune response), apoptosis, and transcription/translation. In adipose tissue, estradiol influenced modules related to focal adhesion, development, metabolism (protein, lipid, oxidative phosphorylation), immune system (complement and coagulation), and translation.
Gonadal sex also showed considerable influence on liver modules related to protein metabolism/assembly, development, stress/immune response, apoptosis, and transcription/translation regulation, whereas in adipose tissue gonadal sex mainly affected developmental and focal adhesion processes, and to a lesser degree, lipid metabolism, biological oxidation, and intracellular signaling modules ( Figure 2G).
The coexpression network analysis also confirmed the limited effect of sex chromosomal variation on altering coexpression modules (Figure 2G, 2H). However, in adipose tissue, sex chromosomes showed weak effects on modules related to lipid metabolism and intracellular signaling when estradiol and blank groups were considered, but not when the testosterone group was included ( Figure 2H).
Overall, the gene coexpression network analysis offered clearer patterns of tissue-specificity and functional specificity of each sex-biasing factor compared to the DEG-based analysis.

Effect of hormonal treatment on gene expression direction across genotypes
Due to the dominant effect of hormonal treatment as compared to gonadal sex or sex chromosome differences based on the above analyses, we further investigated the differences between testosterone and estradiol treatments in terms of the gene sets they target and the direction of gene expression change within and between tissues.

Overlapping DEGs between testosterone and estradiol treatment
Between testosterone and estradiol DEGs informed by 3WA (Supplementary Figure 2), 226 overlapped in the liver and 383 overlapped for adipose tissue. However, estradiol DEGs in individual genotypes did not overlap with those caused by testosterone in 1WA (p<0.05) in any genotype group (Figure 3). In particular, for the XYF mouse we found no overlapping DEGs in either the liver or adipose DEGs between testosterone and estradiol  The bar graphs focused on the DEGs that passed an FDR < 0.05 and were overlapping between testosterone and estradiol treatment for each genotype and tissue. To understand the effects of each hormone, we plotted the log2 fold change of the hormonal effects. The Venn diagrams showcase comparison of DEGs of T effect vs E effect, as well as the top 5 up and down regulated genes for T or E in liver or adipose tissue for each genotype. There was no statistically significant overlap between any comparisons in the Venn diagrams.
Another observation is that DEGs in XXF adipose and XXM liver have noticeably large (>250 fold) increases in gene expression affected by the two hormones (Figure 3A; 3F). These include Lcn13 (Lipocalin-13, involved in glucose metabolism) and Slco1a1 (Solute Carrier Organic Anion Transporter 1A1, important in Na + independent transport of anions such as prostaglandin E2 and taurocholate) in XXM liver, and Mpp5 (Membrane Palmitoylated Protein 5, important in cell junction organization and PI3K-AKT signaling) and Rfwd3 (Ring finger and WD repeat domain 3, important in ligase activity and p53 binding) in XXF adipose.

Top DEGs specific for estradiol and testosterone for liver and adipose tissue
Besides the overlapping DEGs between hormones, we investigated the top DEGs that were specific for In adipose tissue, the top DEGs for testosterone and estradiol generally varied across genotypes. Limited consistency in testosterone DEGs was seen for Mt2 (Metallothionein-2, important in the detoxification of heavy metals) between XXM and XYM, and for estradiol treatment, consistency was seen between XXM and XYM for Hpca (Hippocalcin, involved in calcium binding) and between XXF and XYM for Thbs1 (Thrombospondin 1, important in platelet aggregation).

Transcription factor (TF) network analysis
To understand the regulatory cascades that explain the large numbers of sex-biased genes affected by hormone treatments, we performed TF analysis using as input DEGs that passed an FDR<0.05 from 1WA ( Table  1; Table 2). For the testosterone liver DEGs, we identified 67, 66, 60 and 62 TFs for XYM, XXM, XXF and XYF respectively (Figure 4A-D; Supplementary Table 2). As expected, we captured gonadal hormone receptors including Androgen Receptor (AR) as a highly ranked TF in all genotypes and estrogen receptors (ESR1, ESR2, ESRRA) to be TFs with lower rank. We also found NR3C1 (Nuclear receptor subfamily 3; the glucocorticoid receptor important for inflammatory responses and cellular proliferation) to be among the top 5 TFs for all four genotypes and the top-ranked TF for XXF and XYF, which is consistent with a female bias for this TF (Oliva et al. 2020). A number of circadian rhythm TFs were found throughout all genotypes in the liver including CRY1, CRY2, PER1, and PER2, which is consistent with sex differences in body clocks (Anderson and FitzGerald 2020).
Additional consistent TFs for testosterone effect in liver across multiple genotypes, where sex bias has been documented to some extent previously, include FOXA1/2, XBP1, HNF4A, SPI1, and CTCF.  ) and humans (Oliva et al. 2020).

Gene regulatory network analysis
An alternative and complementary approach to the TF analysis above is to utilize a gene regulatory network approach to decipher the key drivers (KDs) that may drive sex-biased gene alterations in each genotype based on the DEGs found in 1WA ( Table 1; Table 2). These KDs may or may not be TFs.
In the liver (Figure 4I), we saw a large overlap between the KDs predicted across all four genotypes for testosterone treatment, but no KDs captured for estradiol treatment. We identified Cyp7b1, which is important in converting cholesterol to bile acids and metabolism of steroid hormones, as among the top 5 KDs for all genotypes.
Mgst3 (involved in inflammation), C6 and C8b (complement genes), and Ces3b (xenobiotics detoxification) were top 5 KDs for 3 of the 4 genotypes ( Figure 4I). We also identified KDs specific to particular genotypes (Supplementary Table 4) such as Es31 (xenobiotics detoxification) for female gonads, Slc22a27 (anion transport) for XXF, Serpina6 (inflammation) for XYF, and Hsd3b5 (steroid metabolism) for male gonads. Among these KDs, Slc22a27 was previously found to be expressed predominantly in females and Hsd3b5 and Cyp7b1 were male specific (Adams et al. 2015), thus agreeing with our results.
For estradiol, 31 KDs were found for adipose tissue DEGs from the XXF, XXM and XYM genotypes ( Figure 4J; Supplementary Table 5). The KDs included Mrc1 (response to infection), which is the only overlapping top KD between genotypes XXF and XXM. KDs that were more highly ranked for XXM but still statistically significant in XXF included genes involved in extracellular matrix organization (Prxx2, Mfap2, Col1a2, and Gas7), and those specific to XXF are relevant to lipid synthesis/metabolism (Tbxas1, Pla1a) and immune function (Emr1 and Ccdc109b). Irf7 is the only KD for XYM, which has been recently suggested to be a TF in adipocytes with roles in adipose tissue immunity as well as obesity (Kuroda et al. 2020).

Disease association of the genes affected by sex-biasing factors
Finally, to test the disease relevance of the genes affected by sex-biasing factors, we used a marker set enrichment analysis (MSEA) to detect whether the DEGs highlighted in 1WA overlap with genes previously identified to have SNPs associated with diseases/pathogenic traits by GWAS (Figure 5). Of the 73 disease/traits screened for which full GWAS summary statistics is available (Supplementary Table 6), we focused on two broad categories, "cardiometabolic" (Figure 5A;5B) and "autoimmune" (Figure 5C;5D) Figure 5A-B). The autoimmune category included Irritable Bowel Disease (IBD), Ulcerative Colitis (UC), Crohn's Disease (CD), and Type 1 Diabetes (T1D) (Figure 5C; 5D). For hormone DEGs, we focused on those that are directly relevant to the general human population to understand how testosterone or estradiol can affect disease outcomes on XYM (physiological males) or XXF (physiological females).

Figure 5. Bar graphs showing enrichment of the hormone DEGs and gonadal DEGs for known cardiometabolic and autoimmune diseases based on MSEA analysis. A. Bar graph for association of adipose
DEGs from the post hoc one-way ANOVA for ten cardiometabolic diseases/traits. B. Bar graph for association of liver DEGs from the post hoc one-way ANOVA with cardiometabolic diseases/traits. C. Bar graph for association of adipose DEGs from the post hoc one-way ANOVA for five autoimmune diseases. D. Bar graph for association of liver DEGs from the post hoc one-way ANOVA for five autoimmune diseases. E. Bar graph for association of adipose DEGs from the two-way ANOVA for ten cardiometabolic diseases/traits. F. Bar graph for association of liver DEGs from the two-way ANOVA f ten cardiometabolic diseases/traits. A-D. DEGs at an FDR <0.05 derived from the posthoc one-way ANOVA were tested against genetic association signals with cardiometabolic and autoimmune diseases and traits. E and F DEGs at an FDR < 0.05 were tested against genetic association signals with cardiometabolic diseases. Dotted line signifies FDR <0.05 and * denotes enrichment minimally below the FDR< 0.05 cutoff.

Disease association for hormone DEGs
When adipose tissue and cardiometabolic diseases were considered ( Figure 5A; Supplementary estradiol DEGs in XYM are enriched for CD and T1D, and those in XXF for CD and IBD. UC appears to be a unique association with testosterone DEGs whereas CD is associated with adipose DEGs of both hormones.
When liver tissue and cardiometabolic traits were analyzed ( Figure 5B; Supplementary Table 9), we see low level enrichment but specificity of testosterone DEGs for both T2D and LDL, whereas no estradiol DEGs were enriched for any of the ten traits. For the association of liver DEGs with autoimmune disease (Figure 5D; Supplementary Table 10), we found that testosterone DEGs in XYM to be related to UC and T1D, whereas DEGs in XXF to be related to IBD, UC and T1D. Interestingly, estradiol DEGs from XYM had no association with autoimmune diseases but DEGs in XXF were associated with all autoimmune diseases. Taken together, liver DEGs in XXF (gonadal females) affected by both testosterone and estradiol were enriched for autoimmune diseases, implying the modifiable effect of either hormone in these diseases in females.
Overall, the disease association analyses identified adipose DEG sets altered by both hormones for both cardiometabolic and autoimmune processes. For liver DEGs, the most significant associations were with autoimmune diseases and subtle T2D and LDL associations were identified for liver testosterone DEGs. It is unclear, however, whether the associations would be pathogenic or protective.

Disease association for gonadal sex DEGs
We also used MSEA to detect whether gonadal DEGs highlighted in 2WA (FDR<0.05) overlap with genes previously identified to have SNPs associated with diseases/pathogenic traits by GWAS (Figure 5E-F). For adipose tissue, the gonadal DEGs were particularly enriched in cardiometabolic disease/trait GWAS when on an estradiol background, capturing significance in CAD, T2D, Childhood BMI and HDL traits. Gonadal DEGs on a testosterone background only showed association with HDL and no disease association was captured when no hormone (blank) is the background (Supplementary Table 11). In regard to the liver, gonadal DEGs on an estradiol background were enriched for TC, TG, and LDL GWAS signals (Supplementary Table 12).

Disease association for sex chromosome DEGs
Due to the low number of DEGs captured for the sex chromosome effect, no enrichment results are possible through MSEA, therefore we queried whether sex chromosome DEGs have been previously implicated in human diseases by overlapping the DEGs (FDR<0.05 from 2WA) with candidate genes from the GWAS catalog for 2203 traits. Both adipose tissue and liver DEGs demonstrating sex chromosome effects overlapped with GWAS candidates for both cardiometabolic and autoimmune disease (Supplementary Table 13) as well as various other diseases (Supplementary Table 14).

Disease association for DEGs showing interactions among sex-biasing factors
Similarly, overlapping DEGs that showed interactions among the sex-biasing factors (Supplementary Table 1) with GWAS candidate genes, demonstrated that both liver and adipose DEGs affected by interaction between gonad and hormone, and between sex chromosome and gonad overlapped with numerous autoimmune and cardiometabolic traits and disease (Supplementary Table 15).

DISCUSSION
The variation in physiology and pathophysiology between sexes is established via the modulatory effects of three main classes of sex-biased agents: sex chromosome complement (XX vs. XY), long-lasting organizational effects of hormones secreted by the gonads at critical periods such as prenatally or at puberty, and the acute, transient activational effects of circulating hormones that wax and wane at various life stages (Arnold 2009). The manifestations of these sex-dependent modulators impact disease incidence and severity, including metabolismrelated diseases and autoimmune diseases (Voskuhl and Gold 2012;Link and Reue 2017). In this study, we separated the effects of these sex-biasing components using the FCG model, thus enabling the analysis of each contributing factor as well as their interactions in altering gene expression in inguinal adipose and liver tissues, which are relevant in systems metabolism and immunity (Figure 1).
Our data revealed distinct patterns between tissues in the relative contribution of each sex-biasing factor to gene regulation ( Table 1; Table 2). In particular, the liver transcriptome is mainly affected by acute effects of testosterone, followed by acute effects of estradiol, organization effect of gonadal sex, and sex chromosome complement, whereas inguinal adipose gene expression is primarily regulated by acute effects of estradiol, followed by gonadal sex, acute effects of testosterone, and sex chromosome complement. The genes and pathways regulated by the sex-biasing factors are largely different between factors, although metabolic, developmental, and immune functions can be regulated by both activational effects of sex hormones and gonadal sex (organizational effects).
Sex chromosome effects were primarily associated with genes that reside on X and Y chromosomes, along with a handful of autosomal genes involved in inflammation and metabolic processes that are downstream of the sexbiasing effects of X and Y genes. Lastly, the liver and adipose tissue genes affected by the sex-biasing factors were found to be downstream targets of numerous TFs and network regulators, not just the sex hormone receptors, and show association with human cardiometabolic and autoimmune diseases (Figure 6).

Figure 6. Study Summary.
Utilizing the FCG model, we separated the three major classes of sex-biased agents and uncovered their relative contribution to transcriptional alterations in the liver and adipose tissue, the resulting biological processes enriched and finally the diseases associated.
Previously, sex differences in the liver transcriptome have been largely attributed to sex differences in the circadian rhythm and levels of Growth Hormone, which are established because of perinatal organizational masculinization of hypothalamo-pituitary mechanisms controlling Growth Hormone (Mode and Gustafsson 2006;Waxman and Holloway 2009;Sugathan and Waxman 2013). Genes regulated in this manner would be expected to appear in the gonadal effect DEGs. Our results suggest, however, that the acute activational effects of gonadal hormones might be a more important influence, because of the larger number of testosterone or estradiol DEGs compared to gonad DEGs. Our results support previous evidence that removal of gonadal hormones in adulthood eliminates most sex differences in mouse liver gene expression (van Nas et al. 2009;Norheim et al. 2019), and that liver-specific knockout of estrogen receptor alpha or androgen receptor altered genes that underlie sex differences in the liver transcriptome (Zheng et al. 2018). It is possible that the effects of gonadal steroids during adulthood are required for some of the organizational effects of testosterone mediated via Growth Hormone action. In contrast to liver, gonadectomy does not eliminate sex differences in the adipose transcriptome (Norheim et al. 2019), which agrees with our finding that the organizational effects of gonads play a strong role, in addition to estradiol, in adipose gene regulation.
Although we observed a large effect of hormonal treatment relative to gonadal sex and sex chromosome complement for both tissues (Figure 2A-F), the striking tissue-specificity for each of the sex-biasing factors observed here highlights that individual tissues have unique sex-biased regulatory mechanisms. Circulating levels of testosterone play a major role in determining sex differences in liver gene expression patterns (Figure 2A), confirming earlier studies (van Nas et al. 2009;Norheim et al. 2019), whereas levels of estradiol as well as gonadal sex are the main factors to drive sex differences in adipose transcriptome ( Figure 2D). For adipose tissue, the stronger estradiol effect than testosterone is consistent with evidence that ERs are expressed at higher levels than ARs in adipose tissue, as well as having higher overall expression in females (Shen and Shi 2015a). In the liver, ERs are equally distributed between males and females but are generally expressed at a low level (Eisenfeld and Aten 1987). However, ARs are present at much higher concentrations (Eisenfeld and Aten 1987) and higher expression levels in males compared to females (Shen and Shi 2015a); this agrees with our TF analysis ranking ARs more highly than ERs in the liver for all genotypes.
We found that the gonadal sex factor primarily affects developmental pathways, cell adhesion, and metabolic pathways in adipose (Table 2; Figure 2H), which corroborates past evidence indicating that early gonadal sex status and associated hormonal release play critical roles in the development of sex differences and disease outcomes (Leung et al. 2004;Varlamov et al. 2012;Shen and Shi 2015b). Interestingly, the gonadal sex factor on gene expression is more prominent in mice treated with estradiol rather than testosterone, suggesting that estradiol may amplify, and testosterone may minimize, the effect of gonadal sex on gene regulation in adipose ( Figure 2B; 2E).
Compared to the organizational gonadal sex effects and activational hormone effects, the sex chromosome effects were minimal, and no coherent pathways were found for the sex chromosome-driving DEGs ( Table 1; Table   2) or co-expression modules (Figure 2G-H). The DEGs include those known to escape X inactivation (Kdm6a, Eif2s3x, Ddx3x) (Chen et al. 2012;Berletch et al. 2015) and their Y paralogues (Eif2s3y, Ddx3y). The X escapees are expressed higher in XX than XY cells, causing sex differences in several mouse models of metabolic, immune, and neurological diseases (Kaneko and Li 2018;Itoh et al. 2019;Davis et al. 2020;Link et al. 2020). Comparison between XXF and XYF revealed that immune and arrhythmogenic right ventricular cardiomyopathy (ARVC) related pathways are enriched among the DEGs affected by sex chromosome when gonadal hormonal levels are absent (blank group) in these mice that were never masculinized by organizational effects of testicular androgens.
This may be relevant to the postmenopausal state in women, which is characterized by an increase in risk of heart and autoimmune disease, potentially driven by sex chromosomal effects.
As our comparative analysis of the three classes of sex-biasing factors clearly determined that the activational effects of gonadal hormones are the dominant factors, we further investigated the main players that may regulate the sex-biased genes, using a gene regulatory network analysis and a TF analysis, revealing both expected and novel findings. In concordance with the importance of hormonal effects and consistent with recent studies searching for tissue-specific sex bias Oliva et al. 2020), TFs for hormone receptors (AR, and ESR1/2) were captured in the majority of genotypes (Figure 4A-4H). Beyond the major hormonal receptors, within the liver numerous circadian related TFs were captured (PER1, PER2, CRY1 and CRY2). Although it is known that males and females have differing biological clocks (Anderson and FitzGerald 2020), the contribution of hormones particularly in this rhythm is far from fully elucidated and our findings support that hormones need to be taken into account in liver circadian rhythm studies. In adipose tissue for estradiol treatment, however, we found that surprisingly the XXF genotype has no significant signal for ERs, which may imply that estradiol's major contribution in adipose gene regulation is more importantly through TFs such as H2AZ, which have been shown to be essential for estrogen signaling and downstream gene expression (Gévry et al. 2009). In addition to TFs, we utilized a GRN analysis, revealing non-TF regulators. For the liver GRN (Figure 4I), key driver genes for testosterone DEGs are involved in immune processes (Mgst3, C6, C8b), steroid metabolism (Cyp7b1 and Hsd3b5), and xenobiotic detoxification (Ces3b and Es31). In adipose tissue (Figure 4J), the network is dominated for estradiol signals likely due to the larger number of DEGs as compared to testosterone treatment. Here, we see far fewer shared key drivers across genotypes relative to the results in the liver with testosterone treatment, indicating that estradiol has more finely tuned interactions with the gonadal sex and sex chromosome genotypes than the broad effect of testosterone. The role of the numerous liver and adipose regulators uncovered in our analysis warrant further investigation.
Lastly, to provide context to the health relevance of the liver and adipose sex-biasing DEG sets, we looked for GWAS association of these genes with human diseases/traits. We found that hormone-affected genes in adipose tissue were enriched for genetic variants associated with numerous cardiometabolic diseases/traits, but the enrichment was weaker for the liver DEGs (Figure 5A; 5B). CAD and its related traits (HDL, LDL, TG, TC, BMI and fasting glucose) showed significant and widespread enrichment in adipose DEGs for both testosterone and estradiol. Another important area of sex difference is found within autoimmunity, which occurs more in females (Mauvais-Jarvis et al. 2020). While both adipose and liver DEGs from multiple hormone-genotype combinations were enriched for autoimmune diseases, the liver DEGs, particularly those from the XXF genotype, had more prominent autoimmune association. Beyond the hormonal DEG enrichment in human disease/trait, we also found that DEGs caused by gonad type from both adipose and liver are involved in cardiometabolic disease (Figure 5E; 5F). Finally, despite minimal DEGs captured for the sex chromosome effect as well as the interactions between each sex biasing factors, we found overlap of these DEGs with various disease traits. The DEGs underlying disease associations may explain the differential susceptibility of males and females to these major diseases, and warrant further investigation to distinguish risk versus protection through the genes identified in this study.
The analyses presented in this study show an extensive dissection of the relative contribution of three classes of sex biasing factors on liver and adipose gene expression, their associated biological processes and regulators, and their potential contribution to disease. Despite retrieving numerous new insights, we acknowledge the following limitations. First, gonadectomy and subsequent treatment of hormones may have caused activational effects that do not match the effects of endogenous physiological changes in the same hormones, leading to more predominant activational effects being observed. Second, the comparison of mice with testes vs. ovaries does not map perfectly onto mice that had organizational effects of testicular vs. ovarian secretions because of the potential effects of the Sry transgene, which was present in tissues only of mice with testes. Lastly, only liver and inguinal adipose tissues were investigated, and other tissues warrant examination in future studies.
Overall, our data revealed tissue-specific differential gene expression resulting from the three sex-biasing factors, thereby teasing apart their relative contributions to the differential expression of key genes in a variety of clinically significant pathways including metabolism, immune activity, and development. Importantly, in addition to establishing the critical influence of hormones and their effect on the transcriptome in a tissue specific manner, we also uncovered and highlighted the underappreciated role of the sex chromosomal effect and organizational gonadal effect as well as interactions among sex-biasing factors in global gene regulation. Our findings offer a comprehensive understanding of the origins of sex differences, and each of their potential associations with health and disease.

Overall Study Design
In this study we used FCG mice on a C57BL/6J B6 background (B6.Cg-TgSry2Ei Srydl1Rlb/ArnoJ, Jackson were euthanized 3 weeks later; liver and inguinal adipose tissues were dissected, snap frozen in liquid nitrogen and stored at -80°C for RNA extraction and Illumina microarray analysis.

Genotyping
DNA was extracted from tails or ears using Chelex resin (Bio-Rad, Hercules, CA). The genotype of mice was determined by PCR using the primers described in (Itoh et al. 2015) and (Burgoyne and Arnold 2016).

RNA isolation, microarray hybridization, and quality control
RNA from liver and inguinal adipose tissue was isolated using Trizol (Invitrogen, Carlsbad, CA). Individual samples were hybridized to Illumina MouseRef-8 Expression BeadChips (Illumina, San Diego, CA) by Southern California Genotyping Consortium (SCGC) at UCLA. Principal Component Analysis (PCA) was used to identify three outliers among the adipose samples, which were removed from subsequent analyses.
The 3WA tested the general effects of 3 factors of sex chromosomes, gonad, and hormonal treatments, as well as their interactions. Two sets of 3WA were conducted, one comparing T vs. B, and the other E vs. B, to discriminate separate E or T effects. The 2WA tested the effects of sex chromosomes and gonads as well as their interactions within each hormonal treatment group (T, E, or B) separately. For 1WA, we tested the effects of T (comparing T vs. B) and E (comparing E vs. B) within each genotype. Multiple testing was corrected using the Benjamini-Hochberg (BH) method, and significance level was set to FDR <0.05 to define significant DEGs. Suggestive DEGs at p<0.01 and p<0.05 were also retrieved for analyses that are less sensitive to DEG cutoffs, such as pathway analysis.

Co-expression network construction and identification of co-expression modules affected by individual sexbiasing factors
As co-expression networks can reveal unique biology that cannot be retrieved by DEG analysis (Huan et al. 2015;Cordero et al. 2019;Zhao et al. 2019), we used the Multiscale Embedded Gene Co-expression Network Analysis (MEGENA) (Song and Zhang 2015), a method similar to WGCNA (Langfelder and Horvath 2008), to recognize modules of co-expressed genes affected by the three different sex-biasing factors (details in Supplementary Methods). The unique strength of MEGENA compared to WGCNA is that it allows a gene to be in multiple modules and produces smaller and more functionally coherent modules. The influence of each sex-biasing factor on the resulting modules was assessed using the first principle component of each module to represent the expression of that module, followed by 3WA, 2WA, 1WA tests and FDR calulation as described under the DEG analysis section to identify differential modules (DMs) at FDR <0.05 that are influenced by the various sex-biasing factors.

Annotation of the pathways over-represented in the DEGs and DMs
For each of the DEG sets and DMs that were significantly affected by any of the three sex-biasing factors, we conducted pathway enrichment analysis against Gene Ontology (GO) Biological Processes and KEGG pathways using Fisher's exact test, followed by BH FDR estimation. Pathways that had >5 overlapping genes with the DEGs or DMs and FDR < 0.05 were used to annotate the functions of the DEGs or modules.

Gene regulatory network analysis
To predict potential regulators of the sex-biased DEGs, we used the Weighted Key Driver Analysis (wKDA) function of the Mergeomics pipeline (Shu et al. 2016) and liver and adipose Bayesian networks. In brief, the Bayesian networks were built from multiple large human and mouse transcriptome and genome datasets (Yang et al. 2006;Wang et al. 2007;Emilsson et al. 2008;Schadt et al. 2008;Tu et al. 2012). To identify the hub genes or key driver genes within these tissue-specific networks, we searched for network genes whose neighbors were enriched for DEGs from the corresponding tissues. Bayesian network key driver genes were considered significant if they passed an FDR<0.05. The Bayesian networks of top key drivers were visualized using Cytoscape (Shannon et al. 2003).

Transcription Factor (TF) Analysis
To predict TFs that may regulate the sex-biased DEGs sets, we used the Binding Analysis for Regulation of Transcription (BART) computational method . We followed the tool's recommendation of a minimum of 100 DEGs as input and an Irwin-Hall p-value cut off (p < 0.01) for identify TFs.

Maker Set Enrichment Analysis (MSEA) to connect sex biasing DEGs with human diseases or traits
To assess the potential role of the DEGs affected by each of the sex-biasing factors in human diseases, we collected the summary statistics of human GWAS for 73 diseases or traits that are publicly available via GWAS catalog (MacArthur et al. 2017). We used the MSEA function embedded in Mergeomics (Shu et al. 2016) to assess the association of the DEGs with human diseases using a Chi-square like statistic (details in Supplementary Methods).