Investigating transcriptome-wide sex dimorphism by multi-level analysis of single-cell RNA sequencing data in ten mouse cell types

Abstract

Background: It is a long established fact that sex is an important factor that influences the transcriptional regulatory processes of an organism. However, understanding sex-based differences in gene expression has been limited because existing studies typically sequence and analyze bulk tissue from female or male individuals. Such analyses average cell-specific gene expression levels where cell-to-cell variation can easily be concealed. We therefore sought to utilize data generated by the rapidly developing single cell RNA sequencing (scRNA-seq) technology to explore sex dimorphism and its functional consequences at the single cell level.

Methods: Our study included scRNA-seq data of ten well-defined cell types from the brain and heart of female and male young adult mice in the publicly available tissue atlas dataset, Tabula Muris. We combined standard differential expression analysis with the identification of differential distributions in single cell transcriptomes to test for sex-based gene expression differences in each cell type. The marker genes that had sex-specific inter-cellular changes in gene expression formed the basis for further characterization of the cellular functions that were differentially regulated between the female and male cells. We also inferred activities of transcription factor-driven gene regulatory networks by leveraging knowledge of multidimensional protein-to-genome and protein-to-protein interactions and analyzed pathways that were potential modulators of sex differentiation and dimorphism.

Results: For each cell type in this study, we identified marker genes with significantly different mean expression levels or inter-cellular distribution characteristics between female and male cells. These marker genes were enriched in pathways that were closely related to the biological functions of each cell type. We also identified sub-cell types that possibly carry out distinct biological functions that displayed discrepancies between female and male cells. Additionally, we found that while genes under differential transcriptional regulation exhibited strong cell type specificity, six core transcription factor families responsible for most sex-dimorphic transcriptional regulation activities were conserved across the cell types, including ASCL2, EGR, GABPA, KLF/SP, RXRα, and ZF.

Conclusions: We explored novel gene expression-based biomarkers, functional cell group compositions, and transcriptional regulatory networks associated with sex dimorphism with a novel computational pipeline. Our findings indicated that sex dimorphism might be widespread across the transcriptomes of cell types, cell type-specific, and impactful for regulating cellular activities.

Keywords: Cell cluster; Differential distribution; Differential expression; Pathway analysis; Sex dimorphism; Single cell RNA sequencing; Transcription regulatory network.

Fig. 1

Schematic overview of processing cell type-specific scRNA-seq data. Expression data were cleaned, imputed, and validated by clustering. Differential distribution analysis using scDD identified differentially expressed genes and various forms of differential distributions. Differentially distributed genes were fed into GSVA to illustrate differential representation of pathways in each individual cell. Gene co-expression information was combined with information of TF-gene motifs and TF-TF interactions in PANDA network inference. Differentially activated TF-gene couplings (edges), differentially active TFs, and genes under differential intensity of TF-regulation were studied and compared among different cell types. TFs are represented by diamonds. Genes are represented by ovals

Fig. 2

Cell type-specific differential gene expression. Volcano plots and GO annotations of differentially expressed genes in a oligodendrocytes and b heart endothelial cells exhibited cell type-specific gene expression and related pathways. Distribution of DE genes showed disequilibrium in females (red) and males (blue). In oligodendrocytes, more genes had female-specific expression, whereas in heart endothelial cells more had male-specific expression. For each cell type, only genes identified by scDD as DE genes having an FDR among the smallest 3% were retained. Only ten enriched GO terms of the smallest p values were presented. GO terms were categorized into biological process (bp), cellular component (cc), and molecular function (mf). Venn diagrams of c brain cells and d heart cells showed imbalanced number and low conservation of DEGs across cell types. Bar graphs record number of DEG in each cell type. In brain cells, 505 DEGs were specific to one cell type, while 65, 22, and 6 DEGs were shared by two, three, and four cell types; in heart cells, 297 DEGs were specific to one cell type, while 41, 7, and 2 DEGs were shared by two, three, and four cell types. “Dendrocyte” in a represents oligodendrocyte, “Cardiac” in b represents cardiac muscle cell, and “Smooth” in b represents smooth muscle cell

Fig. 3

Sub-classification of fibroblast cells identified sub-cell types. a Five sub-cell types of heart fibroblast cells were identified by unsupervised clustering. Most female cells formed cluster 0 and 1 while most male cells formed cluster 2 and 3. b Distribution of male/female cells in each cluster. Cluster 0, 1, 2, and 3 exhibited dominance of one sex, while cluster 4 was a mixture of both sexes. c Expression of top ten marker genes distinguishing each cluster from the other clusters indicated strong sub-cell type specificity. Nonetheless, expression patterns were perceptually more similar between cluster 0 and 2 as well as those between cluster 1 and 3, than other pairs of clusters. Cluster 4 stood out as almost all corresponding marker genes were hardly expressed in other four clusters. d Top 20 marker genes distinguishing cluster 0 from 1 (left) and cluster 2 from 3 (right) showed a strong overlap. Twelve marker genes were conserved in top 20 most significant marker genes for both pairs of comparisons

Fig. 4

Differential edges influenced functional pathways in a cell type-specific manner. Volcano plots of edges in a cardiac muscle cells and b endocardial cells visualized differential edges with an FDR < 5 × 10⁻⁵ and a mean edge weight difference > 0.25 between females (red) and males (blue). Top five significantly enriched GO terms with smallest p values in female-specific edges (red letters) and male-specific edges (blue letters) were selected, respectively. Terms were ordered with regard to their statistical significance, female-specific term with the smallest p value at the top and male-specific term with the smallest p value at the bottom. c Network containing sex-specific edges in leukocytes was illustrated as it displayed moderate complexity suitable for visualization. Labeled nodes are TFs and unlabeled nodes are genes. Female-specific edges are colored red and male-specific edges are colored blue

Fig. 5

Core TF families were conserved in sex-specific regulatory TF-gene pairs. One row represents one TF. Light gray bins indicate that the corresponding TF was not involved in any sex-specific edge in the corresponding cell type. In each cell type, TFs were ranked based on the number of sex-specific edges from which they extend. The more sex-specific edges that a TF contributed to the PANDA network, the higher its ranking is. TFs were arranged in alphabetic order for illustration of TF families. TFs families containing member(s) ranked top 5 (colored red) in any one type of cell were viewed in zoomed-in windows. These TFs suggested seven core TF families showing remarkable property of conservation in all types of cells

Fig. 6

Significantly differential targeting TFs in astrocytes (left) and cardiac muscle cells (right). Each row represents one TF and each column represents one network. In each plot, 100 randomly generated female-specific networks were placed on the left half and 100 randomly generated male-specific networks were placed on the right without clustering. Relative out-degree values were obtained by subtracting median summed edge weight of each TF in all 200 networks of astrocytes/cardiac muscle cells from each individual summed edge weight. Rows were hierarchically clustered with distance measured by Pearson correlation. Bars mark TFs that were significantly more active in either females (red) or males (blue)

Fig. 7

Differentially targeted genes in astrocytes were involved in both life-sustaining and disease-relevant pathways. a Bubble plot shows significantly enriched KEGG pathways (having at least five annotated genes) with an FDR < 0.1. Gene ratio stands for the proportion of differentially targeted genes in corresponding pathways. KEGG pathways of b thermogenesis and c Huntington disease are illustrated. Differentially targeted genes are colored red (female-specific) or blue (male-specific). Saturation is commensurate to the difference in summed edge weight