Differential causal networks highlight sex-based differences in human tissues.

Figure 2 illustrates the number of differentially expressed genes and the corresponding number of filtered genes, considering only the top 20 tissues. Pancreas tissue exhibited the highest number of sex-specific differentially expressed genes, followed by subcutaneous adipose tissue. In contrast, brain-associated tissues displayed substantially fewer differentially expressed genes, with blood tissue showing the lowest count overall, as detailed in the Table 1 . Comparison of differentially expressed and filtered genes between sexes in the top 20 tissues, excluding Y chromosome and low-expression genes. Data summary showing the number of selected samples and the corresponding identified differentially expressed genes for each tissue type analyzed, sorted in descending order Focusing on DCN results, DCNs facilitate the exploration of causal edges that differ between two conditions, in this case, sex. The numbers of edges in the three types of computed DCNs are summarized in Fig. 3 . Dot plot illustrating the number of edges across various tissues, categorized by DCN Type: Symmetric (first line), F–M (second line), and M–F (third line), with dot size indicating edge count and tissues sorted by total descending order based on the total number of DCN edges. Since the symmetric difference represents the union of results achieved by swapping the minuend network and the subtrahend network, the number of edges in all symmetric DCNs equals the sum of the other two DCNs for each tissue, respectively. Furthermore, given the non-zero differences, it can be stated that the causal relationships between the differentially expressed genes differ between men and women, i.e. there are edges that are present in the male sex, but not in the female sex and vice versa. Additionally, to complete the analysis, functional and gene enrichment analyses were performed. Functional enrichment allows for the contextualization of the evidence found, helping to identify significant pathways in which nodes are involved. On the other hand, gene enrichment analysis may reveal the functional significance of GE changes. The use of both STRING and classical gene enrichment analyses within Cytoscape is beneficial, as they are complementary. Although STRING focuses on protein-protein interactions, hence functional networks and possible therapeutic targets, classical gene enrichment identifies overrepresented pathways and biological functions associated with particular genes, thus contextualizing findings. Together, they enhance the understanding of complex biological mechanisms and support the validation of results, making them valuable tools for comprehensive genomic analyses. The enrichment analysis shows an overlap between known biological pathways and those highlighted through the structure of DCNs, underscoring the biological relevance of the inferred causal relationships. The Fig. 4 illustrates the frequency of pathways enriched in more than 5 tissues, taking into account all analyzed tissues in all DCNS. Bubble plot showing pathways occurring more than five times, with bubble size and a purple-to-yellow gradient reflecting frequency. Haptoglobin-hemoglobin complex, Extracellular region pathways followed by Immunoglobulin, B cell receptor (BCR) signaling and Blood microparticle pathways were the most frequent. Extracellular pathways involve various mechanisms outside the cells, implicating extracellular vesicles, exocytosis and cell signaling, which are crucial for cell communication [ 26 ], regulation, and homeostasis. Immunoglobulin pathways include a class of antibodies for infection control being part of the human immune system. Immunoglobulins are vital to the immune system, functioning as both stimulators and inhibitors of inflammation, and their positive effects are noted across multiple diseases [ 27 ]. BCR signaling pathways are important for the activation of B cells with subsequent immune response (adaptive immunity). In addition, BCR pathway is central to the onset of B-cell malignancies and consequently, could be a therapeutic target [ 28 ]. Blood microparticles, small vesicles released from cells into the bloodstream are involved in pathways related to inflammation, coagulation, cellular communication and apoptosis. High level of blood microparticles are usually related to systemic inflammation and thrombotic disorders [ 29 ]. In the Fig. 5 , the frequency of the most represented gene enrichment pathways, identified in more than 5 tissues, is reported, considering all the analyzed tissues in all the DCNs. Bubble plot showing GE pathways occurring more than five times, with bubble size and a purple-to-yellow gradient reflecting frequency. The most frequent pathways refer to the inactivation of the sex chromosome, the binding of RNA with poly(C) sequences, and the HOXA10 gene, which has been shown to be involved in developmental processes and is regulated by sex steroids [ 30 ]. The motif NRTCGTAANN is a specific sequence pattern recognized in these latter pathways. The inactivation of one X chromosome in females (XX) serves to balance the dosage of sex chromosomes [ 31 ]. Therefore, the presence of these pathways could serve as evidence of the validity of the proposed methodology, demonstrating its ability to identify different relationships between men and women. The poly(C) binding proteins are involved in mRNA stabilization, translational activation, and translational silencing, for which the pathways that comprise them are quite complex for the role they play [ 32 ]. HOXA10 gene with different motifs is associated to pathways for gene regulation, cell migration and cancer progression. In particular, there is an involvement of these gene pathways with both women endometriosis [ 33 ] and prostate cancer [ 34 ]. Further details on the pathways identified in the DCNs of individual tissues are provided in the following subsections and in Supplementary material. The top four significant pathways (with the lowest FDR values for pathways identified by String enrichment and the lowest p-values for pathways identified by gene enrichment) highlighted by the applied enrichments are shown. The representations of all tissues DCNs are reduced for better visualization. Only the top 2 genes (nodes) with the highest degree and their 1-hop neighbors for both incoming and outgoing edges are considered. An exception is made for tissues related to the brain, kidney cortex and blood tissues given the small size of the networks. Pancreatic tissue exhibited the highest number of differentially expressed genes between sexes. As a key component of the endocrine system, the pancreas is implicated in diseases such as diabetes mellitus and cancer, where sex-specific GE may influence disease susceptibility and immune responses [ 35 , 36 ]. Investigating the underlying causal relationships and associated pathways is therefore of particular interest. Pathway enrichment analysis, summarized in Table 2 , revealed that the most significantly affected pathways in the DCNs are immune-related, consistent with the central role of inflammation and immune dysfunction in both diabetes and pancreatic cancer [ 37 , 38 ]. Pancreas significant pathways In addition, gene enrichment analysis indicates which genes belong to interrelated functions, as illustrated in the Table 3 . In fact, as already shown in the pathway enrichment, genes related to the immune response turn out to be among the most significant (lower p-value), and, as mentioned, have a non-secondary role in pancreatic-related diseases. Given the size of the DCNs obtained for this tissue, having the largest number of differentially expressed genes, reduced DCNs have been depicted in the Figs 6 – 8 . Pancreas gene enrichment Male-Female DCN for pancreas tissue. Female-Male DCN for pancreas tissue. Symmetric DCN for pancreas tissue. Adipose tissue, composed predominantly of adipocytes, was analyzed in two forms: subcutaneous and visceral (omentum). Pathway enrichment analysis of subcutaneous adipose tissue identified three principal pathways linked by their roles in inflammation and lipid metabolism ( Table 4 ). Serum amyloid A proteins, key mediators of the acute-phase response, also regulate cell–cell communication and modulate inflammatory, immunologic, neoplastic, and protective processes [ 39 ]. Chronic activation of the acute-phase response by adipose tissue has been implicated in the elevated inflammatory state of diabetes and its cardiovascular sequelae [ 40 ]. Adipose subcutaneous significant pathways Additionally, pathways related to high-density lipoprotein (HDL) metabolism were enriched. Given the critical role of HDL in cholesterol homeostasis and cardiovascular protection [ 41 ], and the importance of adipose tissue in HDL dynamics, these findings offer insights into potential therapeutic targets aimed at enhancing cardiometabolic health. Gene enrichment for adipose subcutaneous tissue ( Table 5 ) shows the presence of pathways related to immune response and inflammatory processes. Moreover, the pathway related to the binding and the uptake of ligands plays a crucial role in lipid metabolism by mediating the uptake and binding of various ligands, including lipoproteins and fatty acids, in adipose tissue. These receptors are critical for maintaining metabolic homeostasis and facilitating lipid transport [ 42 ]. Given the size of the DCNs obtained, reduced DCNs have been depicted in the Figs 9 – 11 . Adipose subcutaneous gene enrichment Male-Female DCN for adipose subcutaneous tissue. Female-Male DCN for adipose subcutaneous tissue. Symmetric DCN for adipose subcutaneous tissue. Table 6 shows significant pathways identified by String analysis applied to adipose visceral omentum tissue. In addition to pathways related to genes involved in sex differentiation, homeobox proteins pathway has also been identified. It involves transcription factors for cell differentiation, development, and regulation also in the early stages of embryonic development [ 43 ]. Moreover, genes involved in the lipid response were detected in the DCNs of the adipose visceral omentum tissue. Reduced DCNs for adipose visceral omentum tissue are illustrated in the Figs 12 – 14 . Adipose visceral omentum significant pathways Male-Female DCN for adipose visceral omentum tissue. Female - Male DCN for adipose visceral omentum tissue. Symmetric DCN for adipose visceral omentum tissue. The lungs, which are central to the respiratory system, demonstrated a significantly lower number of differentially expressed genes compared to other analyzed tissues. This finding may be partially due to the prolonged ischemic time commonly associated with GTEx lung samples, which can result in increased RNA degradation and, subsequently, reduced RNA quality [ 44 ]. Reduced directed CNs (DCNs) for lung tissue are shown in Figs 15 – 17 . Male-Female DCN for lung tissue. Female-Male DCN for lung tissue. Symmetric DCN for lung tissue. STRING functional enrichment analysis for lung tissue did not identify significant pathways. However, gene enrichment analysis ( Table 7 ) revealed genes implicated in maintaining X-chromosome-related transcriptional suppression (such as HNRNPU, a gene encoding for the nuclear ribonucleoprotein U). The enrichment results suggest the presence of key gene regulatory processes and RNA interactions in lung tissue. Notably, the transcription factor HOXA10, frequently upregulated in lung adenocarcinoma and implicated in cancer progression [ 45 , 46 ], emerged alongside processes related to sex chromosome inactivation and RNA binding, pointing to potential mechanisms underlying lung-specific molecular regulation. Lung tissue gene enrichment Lung tissue samples from the GTEx project are known to exhibit reduced RNA integrity due to extended ischemic intervals, potentially leading to RNA degradation and introducing systematic biases in transcriptomic measurements. In this study, additional correction procedures beyond the standard GTEx quality control pipeline were not implemented, as the primary objective was to facilitate exploratory and comparative analyses across tissues under consistent preprocessing conditions. Nevertheless, the notably lower number of differentially expressed genes and the increased sparsity of the inferred CNs in lung tissue align with the anticipated effects of compromised RNA quality. Despite the relatively low number of differentially expressed genes, causal discovery was performed using a robust ensemble strategy, and differential edges were derived from consistent differences between male and female network topologies. This supports the biological relevance of the identified sex-specific regulatory patterns in this tissue. GE from 11 brain regions (the amygdala, anterior cingulate cortex, caudate, cerebellar hemisphere, cerebellum, cortex, frontal cortex (BA9), hippocampus, hypothalamus, putamen, and cervical spinal cord -C1-) were analyzed. Pathway enrichment analyses were conducted for each region. Results for the cortex, cerebellum, cerebellar hemisphere, caudate, hypothalamus, frontal cortex (BA9), cervical spinal cord (C1), and amygdala are summarized in Table 8 . Reported p-values represent the mean values across these tissues. Brain tissues gene enrichment The brain hippocampus tissue showed only genes related to the pathways involving inactive sex chromosomes, whereas the brain putamen basal ganglia exhibited different pathways, as depicted in Table 9 . All the DCNs for the brain tissues are illustrated in Figures 18 – 50 . Male-Female DCN for brain amygdala tissue. Female-Male DCN for brain amygdala tissue. Symmetric DCN for brain amygdala tissue. Male-Female DCN for brain anterior cingulate tissue. Female-Male DCN for brain anterior cingulate tissue. Symmetric DCN for brain anterior cingulate tissue. Male-Female DCN for brain caudate basal ganglia tissue. Female-Male DCN for brain caudate basal ganglia tissue. Symmetric DCN for brain caudate basal ganglia tissue. Male-Female DCN for brain cerebellar hemisphere tissue. Female-Male DCN for brain cerebellar hemisphere tissue. Symmetric DCN for brain cerebellar hemisphere tissue. Male-Female DCN for brain cerebellum tissue. Female-Male DCN for brain cerebellum tissue. Symmetric DCN for brain cerebellum tissue. Male-Female DCN for brain cortex tissue. Female-Male DCN for brain cortex tissue. Symmetric DCN for brain cortex tissue. Male-Female DCN for brain frontal cortex tissue. Female-Male DCN for brain frontal cortex tissue. Symmetric DCN for brain frontal cortex tissue. Male-Female DCN for brain hippocampus tissue. Female-Male DCN for brain hippocampus tissue. Symmetric DCN for brain hippocampus tissue. Male-Female DCN for brain hypothalamus tissue. Female-Male DCN for brain hypothalamus tissue. Symmetric DCN for brain hypothalamus tissue. Male-Female DCN for brain putamen basal ganglia tissue. Female-Male DCN for brain putamen basal ganglia tissue. Symmetric DCN for brain putamen basal ganglia tissue. Male-Female DCN for brain spinal cord cervical tissue. Female-Male DCN for brain spinal cord cervical tissue. Symmetric DCN for brain spinal cord cervical tissue. Brain putamen basal ganglia tissue gene enrichment Reduced DCNs for the heart left ventricle and atrial appendage tissues are illustrated in Figs 51 – 56 . In the left ventricle, gene enrichment analysis identified pathways associated with protein synthesis (poly(C) RNA binding), cell proliferation (HOXA10 regulatory pathway), and X chromosome inactivation ( Table 10 ), similarly observed in some brain tissues. For the atrial appendage, pathways related to adipocyte biology were enriched ( Table 11 ), consistent with evidence that abnormal adipose tissue accumulation can impact cardiac function and promote thrombus formation [ 47 ]. Male-Female DCN for heart atrial appendage tissue. Female-Male DCN for heart atrial appendage tissue. Symmetric DCN for heart atrial appendage tissue. Male-Female DCN for heart left ventricle tissue. Female-Male DCN for heart left ventricle tissue. Symmetric DCN for heart left ventricle tissue. Heart left ventricle tissue gene enrichment Heart atrial appendage tissue significant pathways

Gene counts and gene TPMs (transcripts per millon) for 59 033 genes across 40 human tissues were extracted using data from the GTEx project (v10 data release) [ 11 ]. For this analysis sex-specific tissues and tissues with few samples were excluded. GE data from 981 subjects (654 males and 327 females) were included in the analysis. Gene counts measure a gene’s expression level in a given sample and are often used in DEA [ 12 ]. Gene TPMs, on the other hand, represent a normalized measure of GE and are helpful in comparing the relative abundance of different genes [ 13 ]. Based on this distinction, gene counts were used for sex gene differentiation, and gene TPMs, filtered by differentiated genes, were subsequently utilized to construct DCNs. Genes expressed on the Y chromosome were also identified and excluded from the analysis. These genes were retrieved from the NCBI Gene Database [ 14 ] using the query Homo sapiens Y chromosome . Both gene symbols and their aliases were considered during this extraction. The exclusion of Y chromosome-related genes in this analysis is due to the need to prevent the differential analysis from being distorted by sex-specific effects. Y-chromosome genes are intrinsically expressed only in males, and their inclusion would lead to obvious and biologically trivial differences in expression and causal structure. The extracted data’s pre-processing phase was divided into two main steps: the identification of differentially expressed genes using gene counts and the filtering and cleaning of TPM data. Subjects with unique values across all fields and complete data in the age and gender fields were included in the study. Duplicate entries and records with incomplete data in the specified fields were excluded. The numbers of males and females were balanced by randomly sampling the larger set to match the size of the smaller set, ensuring that differences in sample sizes between sexes would not affect subsequent analyses. The PyDESeq2 library [ 15 , 16 ] was used for DEA between sexes. Genes showing differential expression were defined as those with an adjusted p-value below 0.05 and an absolute fold change greater than 1. Bootstrap resampling involved randomly selecting 90% of the samples with reinsertion during each iteration, ensuring that the proportion of males and females remained constant. Across 100 iterations, the frequency with which each gene was identified as significant was tracked. Only genes appearing as significant in at least 75% of iterations were considered robust and included in the final results. As mentioned earlier, a balanced number of samples from males and females was selected prior to the gene differentiation phase to mitigate the potential influence of sample size disparities between sexes. In addition, a possible distribution shift effect could, in theory, influence the causal discovery phase, although in practice, thanks to the rigorous selection and pre-processing procedures of GTEx data, the risk of introducing bias is very low. The subsequent preprocessing of TPM data was conducted using the differentially expressed genes identified in the previous stage as a filter. The data were subsequently categorized into two distinct groups (males and females) and further refined by excluding genes that had zero counts in more samples than the total number of samples in each dataset minus ten. As a result, CNs for men and women will share the same number and type of genes, but will show different causal edges between them, allowing the DCNs computation. A DCN is produced by subtracting one CN from another or performing a symmetric difference between two CNs with the same nodes but different edges. The PC (Peter–Clark) algorithm [ 17 ] was applied to perform the causal discovery phase on filtered TPMs data, obtaining four CNs for each tissue: two for males (one for each age range) and two for females (one for each age range). As a next step, as proposed in [ 6 ], three types of differences between male and female CNs for each tissue were computed: Female CN is subtracted from the male CN: the (Asymmetric) DCN contains all the edges of male CN that are not in female CN; Male CN is subtracted from the female CN: the (Asymmetric) DCN contains all the edges of female CN that are not in male CN; Symmetric difference between male and female CN is constructed: the (Symmetric) DCN contains all the edges of the female CN that are not in the male CN and all the edges of the male CN that are not in the female CN. In other words, it is the union of all the edges in the female CN and the male CN described above. In summary, the investigation of sex-specific mechanisms was conducted using three distinct DCN types constructed based on the presence or absence of causal edges inferred independently in male and female subgroups using the stable PC algorithm. Specifically, edges classified as M–F represent causal interactions detected only in the male subgroup, indicating male-dominant mechanisms; F–M edges, conversely, denote interactions unique to females; while symmetric edges correspond to the union of both sex-specific interactions. In addition, considering the definitions proposed in [ 18 ], the DCNs performed in this study consider the difference in support: where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathit{E^{\textrm{diff}}}$\end{document} is the edge difference set, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $A^{(1)}$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $A^{(2)}$\end{document} denote the adjacency matrix respectively of the graphs \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $G^{(1)}$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $G^{(2)}$\end{document} of which the difference is performing [ 18 ]. This definition considers structural differences between networks, in line with what has been proposed in [ 6 ]. PC algorithm (from the names of its creators, P.S. and C.G.) was used to derive the structure of CNs [ 17 ]. Causal directionality was inferred using the PC algorithm, which determines causal structure by testing for conditional independence rather than relying on correlations alone [ 17 ]. Assuming causal sufficiency, faithfulness, and the Causal Markov condition, the algorithm identifies v-structures and applies orientation rules to resolve the direction of causal edges. This approach enables the discrimination of genuine causal relationships from spurious associations. The CNs were obtained using the gCastle Python package [ 19 ], which implements different algorithms for PC. In this work, the stable variant of the PC algorithm (variant = “stable”) was adopted. It is an efficient method for causal discovery that eliminates order-dependencies in its implementation [ 20 ]. The output of the PC algorithm can be order-dependent, i.e. the order in which the variables are presented can affect the results. The authors of [ 20 ] proposed the stable version of it with several modifications to the algorithm to remove this dependency in the various stages of the algorithm and ensure reliable results. The default parameter settings were kept: ci test = “fisherz” (independent test), alpha = 0.05 (significance level), and no background information was added (priori knowledge: None). Additional evidence supporting this approach is available in the supplementary material. For each tissue, each of the two groups (men and women) is further divided into two age sets: subjects aged between 20 and 49 years and subjects aged between 50 and 79. Therefore, a CN is generated. To enhance the robustness of the resulting CNs, the PC algorithm is applied repeatedly on each group (100 iterations). Only causal edges present in at least 70% of the iterations were considered eligible as causal relationships between genes. The result of this step is four CNs that are subsequently aggregated by combining the relationships found for the two age ranges while maintaining the distinction between sexes. As the last step of this phase, the DCNs were calculated using the NetworkX Python library [ 21 ] to identify the differences in causal relationships between two conditions. Pathway enrichment was performed by utilizing the KEGG pathway database [ 22 ] and the Biological Process category of Gene Ontology [ 23 ], providing a way to analyze and interpret the functions of genes based on their annotations [ 24 ]. Functional enrichment analysis was conducted using the STRING enrichment app of Cytoscape software, while gene enrichment analysis was performed on gProfiler [ 25 ].

As widely discussed, this study underscores the powerful application of DCNs in analyzing human tissues to uncover sex-specific causal relationships among genes. By leveraging the comprehensive GTEx dataset and selecting 40 tissues, the work highlighted significant GE differences between males and females. These results are reinforced by the use of bootstrap resampling in gene differential analysis and multiple iterations in CN generation. The insights gained emphasize the critical role of sex as a determinant in understanding the underlying mechanisms of both healthy and pathological states, paving the way for the development of targeted therapies and personalized medicine. While sex emerged as the primary factor influencing causal relationships, age-related differences were also considered and accounted for, adding an additional layer of validation to our results. The biological insights provided by enrichment analysis following the DCN analysis offer a solid foundation for future research and clinical trials aimed at validating these findings. In essence, the integration of advanced analytical techniques and complex datasets represents a significant advancement in biomedical research, enhancing the comprehension and the treatment of diseases. This is the first comprehensive study performing Differential Causal Network (DCN) analysis across the entire GTEx dataset , considering 40 human tissues . Identification of sex-specific causal relationships in gene expression through a novel application of DCNs. Integration of causal discovery techniques with differential expression analysis, ensuring robustness through bootstrap resampling . Identification of structural differences (edges) in gene causal networks between males and females, with implications for precision medicine , by providing a causal understanding of sex-based biological differences.

This study provides insights into the regulatory mechanisms highlighting the complexity and specificity of genetic interactions in determining both biological and pathological states. The analysis consistently identified strong and reproducible interactions between XIST and TSIX genes that play a central role in the X-chromosome inactivation process, as evidenced by pathway enrichment. As shown in Fig. 57 , comparing DCNs between men minus females, only the connection from XIST to TSIX was found in one network. However, females minus males DCNS showed a total of 23 edges linking these two genes, with 22 edges directed from XIST to TSIX and 21 edges directed from TSIX to XIST. Additionally, the presence of bidirectional interactions (20 bidirectional edges reported) suggests the existence of feedback loops or mutual regulation in the inactivation process. As inactivation of the sex chromosome is a specific biological process in women, these findings may highlight the strength of our approach in identifying sex-specific differential causal interactions. Comprehensive results of DCN connections between XIST and TSIX across all analyzed tissues. Alt Text: Furthermore, the refined view of the top 20 genes connected to XIST and TSIX, as illustrated in Fig. 58 , reinforces the notion that these two genes are pivotal nodes within broader regulatory networks. Overview of top 20 genes having incoming/outgoing edges to/from XIST and TSIX in DCNs across all analyzed tissues. In comparison to previous research, our study not only identifies associations but also clarifies the directionality and strength of these relationships, providing a more nuanced view of gene regulation. This methodological advancement is crucial for accurately interpreting genetic data, particularly in the context of complex diseases involving multiple genes and pathways. The CNs identified here serve as fundamental tools for dissecting the genetic basis of such diseases, potentially revolutionizing our approach to diagnosis and treatment. Despite the innovative techniques and significant insights gained, our approach does have limitations. The construction and analysis of CNs in high-dimensional datasets are challenging, requiring careful consideration to avoid overfitting and ensure the findings’ robustness. Nevertheless, these challenges are mitigated by the comprehensive nature and scalability of the developed methodologies, providing a strong foundation for further research and application in the field. Looking ahead, the implications of this research are substantial. By extending the application of CNs to other biological systems, it could be possible to facilitate the understanding of a wide range of diseases, leading to breakthroughs in genetic research and pharmaceutical development. Moreover, integrating multi-omic data with the presented CNs could enhance the precision and accuracy of the proposed models, facilitating the development of predictive tools for disease progression and treatment responses. Furthermore, this study lays the groundwork for transforming clinical practices by applying genetic insights to patient care. The causal models developed here could eventually guide clinical decision-making, enabling the practice of truly personalized medicine. By predicting individual responses to treatments based on genetic profiles, the therapeutic outcomes could be optimized and the adverse effects could be minimized, fundamentally changing the healthcare landscape. This work was designed as an exploratory and systematic application of DCNs across the entire GTEx dataset, with the goal of highlighting both physiological and potentially pathological causal differences between sexes across tissues. A promising direction for future work is to apply this methodology to disease-specific datasets, allowing a more direct link between DCNs and clinical outcomes in precision medicine.

Differential gene expression (GE) analysis is essential for identifying genes whose transcription patterns vary between groups, such as between sexes. Sex-biased genes exhibit differential expression between males and females, characterized as male-biased (higher expression in males) or female-biased (higher expression in females). However, sex-biased expression is not an intrinsic or fixed property of genes; instead, it is observed within specific samples under particular conditions. For instance, sex-biased GE can vary during development, reflecting growth stage specificity [ 1 ]. Sex-based differences in GE are observed across numerous human tissues. These differences influence the onset, prevalence, and severity of diverse diseases and are modulated by factors such as hormones, sex chromosomes, behavioral differences, and environmental exposures. Notably, these differences are driven by tissue-specific effects rather than merely differences in cell type abundances [ 2 ]. In this context differential network analysis (DiNA) allows for comparing biological networks across different conditions, highlighting unique characteristics. DiNA usually focus solely on comparing the network architectures derived from the same set of genes across conditions without an explicit differential expression step [ 3 ]. However, DiNA could integrates differential expression analysis (DEA), identifying genes with significant changes in expression. Lichtblau et al., combined DiNA with DEA to capture complex features associated with changes in biological networks [ 3 ]. Despite such advances, causal mechanisms underlying these differences remain poorly understood. Causal networks (CNs) model causal relationships, offering a potential starting point for identifying complex and divergent interactions among genes from experimental data [ 4 , 5 ]. Consequently, the need to compare different causal mechanisms in two different experimental condition has lead to the introduction of differential causal networks (DCNs) [ 6 ]. Other studies have primarily focused on identifying genes that display sex-specific patterns of expression, both at a large scale across multiple tissues and within individual tissues [ 7 , 8 ]. These efforts have provided important insights into the differential gene regulation associated with sex. In contrast, our work introduces a novel perspective by analyzing the distinct causal mechanisms underlying these sex-biased expression patterns, offering a deeper understanding of the biological processes driving such differences. In addition, other works have investigated sex differences in microRNA expression, which play key regulatory roles in GE, or copy number alteration [ 9 , 10 ]. In contrast, our analysis is restricted to gene transcripts, allowing us to specifically examine sex-related differences at the level of mRNA expression and their underlying causal mechanisms. DCNs framework was introduced in [ 6 ] to uncover not only whether two CNs differ, but also how they differ structurally, characterizing changes in network topology, such as the presence or absence of specific causal edges. This study focuses on understanding differences in GE between the sexes using DCNs on most human tissues using publicly available Open Access Data from the Genotype-Tissue Expression (GTEx) portal. Based on differentially expressed genes, four CNs were constructed for each of the 40 selected tissues: one for males and one for females, further stratified by age groups. The age specific networks were aggregated by sex to derive two comprehensive CNs (one for males and one for females). Subsequently, three types of differential analyses [ 6 ] were performed to generate the DCNs. Figure 1 summarizes the workflow followed in this study. Illustration of the DCN generation process, emphasizing how causal relationships were derived in the study’s analysis.