Skip to main content

Sex-biased and parental allele-specific gene regulation by KDM6A

Abstract

Background

KDM6A is a demethylase encoded by a gene with female-biased expression due to escape from X inactivation. Its main role is to facilitate gene expression through removal of the repressive H3K27me3 mark, with evidence of some additional histone demethylase-independent functions. KDM6A mutations have been implicated in congenital disorders such as Kabuki Syndrome, as well as in sex differences in cancer.

Methods

Kdm6a was knocked out using CRISPR/Cas9 gene editing in F1 male and female mouse embryonic stem cells (ES) derived from reciprocal crosses between C57BL6 x Mus castaneus. Diploid and allelic RNA-seq analyses were done to compare gene expression between wild-type and Kdm6a knockout (KO) clones. The effects of Kdm6a KO on sex-biased gene expression were investigated by comparing gene expression between male and female ES cells. Changes in H3K27me3 enrichment and chromatin accessibility at promoter regions of genes with expression changes were characterized by ChIP-seq and ATAC-seq followed by diploid and allelic analyses.

Results

We report that Kdm6a KO in male and female embryonic stem (ES) cells derived from F1 hybrid mice cause extensive gene dysregulation, disruption of sex biases, and specific parental allele effects. Among the dysregulated genes are candidate genes that may explain abnormal developmental features of Kabuki syndrome caused by KDM6A mutations in human. Strikingly, Kdm6a knockouts result in a decrease in sex-biased expression and in preferential downregulation of the maternal alleles of a number of genes. Most promoters of dysregulated genes show concordant epigenetic changes including gain of H3K27me3 and loss of chromatin accessibility, but there was less concordance when considering allelic changes.

Conclusions

Our study reveals new sex-related roles of KDM6A in the regulation of developmental genes, the maintenance of sex-biased gene expression, and the differential expression of parental alleles.

Highlights

  • The majority of KDM6A targets are regulated from one allele or the other with allelic downregulation showing a maternal allele bias.

  • Kdm6a knockout affects genes that have been implicated in phenotypes of Kabuki syndrome caused by KDM6A mutations in human.

  • KDM6A facilitates sex-biased gene expression, particularly female-biased expression.

Introduction

Embryonic development is governed by specific patterns of epigenetic modifications to control gene expression. For example, timely activation of Hox gene expression via resolution of bivalency in early development is facilitated by the histone demethylase KDM6A that removes methylation at lysine 27 of histone H3 [1,2,3,4]. Activation of gene expression by KDM6A is implicated in many developmental pathways such as myogenesis, cardiac development, neural stem cell differentiation, and immune cell functions [5,6,7,8,9]. KDM6A can also enhance gene expression independent of its demethylase activity by association to chromatin remodeling complexes that acetylates H3K27 at active enhancers through recruitment of p300 to the MLL complex [10,11,12]. In addition to its major role in gene activation, evidence suggests that KDM6A can also act as a repressor of gene expression, highlighting the complexity of this protein [12, 13].

The higher level of KDM6A in female versus male mammals due to escape from X inactivation suggests a role in sex-specific gene regulation and in sex differences in health and disease [14,15,16]. Kdm6a modulates sex differences in autoimmunity and cancer, and more recently, it has been implicated in sex differences in Alzheimer’s disease [8, 17, 18]. KDM6A is encoded by a dosage-sensitive X-linked gene and has a Y-encoded homolog (UTY) that has little or no demethylase activity [19,20,21]. Knockout of Kdm6a is homozygous lethal in female mice, while a small number of runt males survive, suggesting that UTY partially compensates for loss of KDM6A possibly via demethylation-independent mechanisms [22].

We hypothesized that increased dosage of Kdm6a in females may lead to gene regulation of sex-biased genes and that differential expression of Kdm6a in male and female germline could potentially lead to parent-of-origin biased gene expression in the next generation. To address these possibilities, we investigated the role of KDM6A in sex-biased autosomal gene regulation as well as expression and epigenetic regulation of maternal and paternal alleles in mouse embryonic stem (ES) cells. Kdm6a knockouts (KO) were obtained by CRISPR/Cas9 approaches in ES cells derived from F1 hybrid mice from reciprocal crosses between C57BL/6 J (BL6) and Mus castaneus (cast) in which alleles can be distinguished based on SNPs. As measured by RNA-seq, sex-biased gene expression was markedly reduced following KO, with a preferential loss of female-biased expression. RNA-seq followed by allele-specific analyses revealed a subset of genes with parental allele-specific expression changes following Kdm6a KO, with more genes downregulated on maternal than paternal alleles.

Methods

Cell lines

BC male (E14) hybrid ES cells as well as CB male hybrid ES cells derived from reciprocal crosses between C57BL/6 J and Mus castaneus were provided by J. Gribnau (Erasmus MC Rotterdam, NL). All ES cells were maintained in a humidified incubator at 37ºC and 5% CO2 in the presence of 1000 U/ml leukemia inhibitory factor (LIF) (Millipore) on a monolayer of inactivated mouse embryonic fibroblasts (MEFs) in ES media containing high glucose DMEM media supplemented with 15% fetal bovine serum (FBS), 1% non-essential amino acids, 10 mg/ml APS, 0.1 mM 2-mercaptoethanol and 25 mM L-glutamine. Medium with serum and LIF was preferred to medium containing MEK1/2 and GSK3β inhibitors (2i), which can lead to DNA methylation changes [23, 24]. Just prior to use, plates were enriched for ES cells by incubation on 0.1% gelatin coated dishes for 1 h to allow MEFs to attach, followed by transfer to fresh gelatin-coated plates for overnight culture. Following expansion, ES cells were split once (1:10) to further reduce potential MEF contamination.

ES cell differentiation

For embryoid body (EB) formation, 2 × 106 wt control and Kdm6aΔE1 (see Fig. 1A and Additional file 1: Fig. S1 for definition of KO clones) cells were cultured on non-adherent bacterial culture dishes without LIF for 6 days, with media changed every two days. For retinoic acid induced neuronal differentiation, 1.2 × 106 wt control and Kdm6aΔE1 cells were grown on gelatin-coated tissue culture plates without LIF and in the presence of 1 μM all-trans retinoic acid for 8 days, with media changed every two days. Pictures were taken at the indicated time points using a Motic AE2000 camera.

Fig. 1
figure 1

Global gene expression changes after Kdm6a KO in BC and CB clones. A Schematic illustrating the Kdm6a CRISPR deletions in BC and CB ES cells. Also included are clone identifiers (see Additional file 8). Exonic (E) deletion results in the loss of about 45 kb in length while the promoter (P) deletion results in a loss of about 4 kb in length. B, C Scatter plots of differential gene expression in three BC male Kdm6aΔ/Y clones (Kdm6aΔE1, Kdm6aΔE3, and Kdm6aΔP4) versus two male BC wt clones (B), and three CB male Kdm6a Δ/Y clones (Kdm6aΔE2.5, Kdm6aΔE2.7, and Kdm6aΔP2.1) versus two CB male wt clones (C). Log2 values of genes with > 1TPM in at least one replicate are shown. Downregulated DEGs are labeled in red, upregulated DEGs in green, and genes lacking differential expression in grey. DESeq2 was employed to identify differentially expressed genes (DEGs) in each cross using an FDR cutoff of < 0.05 and a ≥ 1.5 fold-change. D Venn diagrams of downregulated and upregulated genes as measured by a ≥ 1.25 fold change in log2 TPM in male BC and CB Kdm6aΔ/Y clones versus wt. E GO analysis of dysregulated genes with a ≥ 1.25 fold change in log2 TPM found in common between the BC and CB crosses. Numbers in parentheses represent the fold enrichment over expected number of genes within a given biological process, while those not in parentheses represent the number of genes in each category. F Log2 expression fold changes for genes dysregulated in male Kdm6aΔ/Y clones derived from BC and CB crosses. The genes shown are associated with processes implicated in phenotypes seen in Kabuki syndrome (*p ≤ 0.05; **p ≤ 0.01; Has2 p = .05). See also Additional file 9

CRISPR/Cas9 gene editing

We chose a dual sgRNA approach because of the capacity to create large deletions (Fig. 1A and Additional file 1: Fig. S1) [25]. In addition, two independent sets of sgRNAs were used to account for off target effects. First, plasmids containing sgRNAs located in exons 2 and 4 were used to delete ~ 45 kb, resulting in multiple stop codons in five of the six possible translation frames. Six-Frame translation analysis was used for in silico verification of newly acquired translation stop codons following Kdm6a targeting. For a second editing approach, we targeted the Kdm6a promoter region including exons 1 and 2. sgRNAs designed using the CHOPCHOP online tool were cloned into the CRISPR/Cas9 plasmid pX330 as described [26,27,28].

ES cells were transfected with sgRNAs containing CRISPR/Cas9 constructs and a plasmid carrying puromycin resistance (pPGKpuro) using UltraCruz® transfection reagent at a 3:1 CRISPR to pPGKpuro ratio in media with no antibiotics. Two days later, cells were selected for 48-72 h with 1 μg/ml puromycin in complete ES media. Cells were allowed to recover in media with antibiotics. Following recovery, cells were cloned into 96 well plates using serial dilutions. Clones were expanded and screened for deletions using PCR and Sanger sequencing. PCR analysis confirmed correct gene editing of exons 2–4 in BC male ES cell clones Kdm6aΔE1 and Kdm6aΔE3 and CB male ES cell clones Kdm6aΔE2.5 and Kdm6aΔE2.7 (Additional file 1: Fig. S1). Note that BC clones Kdm6aΔE1 and Kdm6aΔE3 were derived from a single CRISPR editing event. Indeed, Sanger sequencing showed that the deletion breakpoints in Kdm6aΔE1 and Kdm6aΔE3 were identical, suggesting that the clones were derived from a single targeting event. PCR amplification and Sanger sequencing confirmed the predicted ~ 4 kb deletion of the promoter region in BC male ES cell clones Kdm6aΔP4 and CB ES clone Kdm6aΔP2.1 (Additional file 1: Fig. S1).

Kdm6a expression in Kdm6aΔE clones was assayed by RT-PCR using primers specific for regions inside the deletion (spanning exons 3–6) and regions downstream of the deletion (spanning exons 23–26). RT-PCR and RNA-seq analyses showed low-level expression of exons downstream of the deletion, but Western blot analysis showed complete ablation of the KDM6A protein, consistent with premature stop codons and published Kdm6a exon KO models (Additional file 1: Fig. S1) [21, 28]. RNA-seq and RT-PCR expression analyses verified the lack of Kdm6a expression in clones with a promoter deletion. Note that the small alternative transcript annotated in UCSC (http://genome.ucsc.edu/) was absent in Kdm6aΔP clones as shown by the lack of reads in this region.

Reverse-transcription quantitative PCR

RNA was extracted from cells using Qiagen RNeasy mini kit followed by cDNA synthesis with the Promega GoScript reverse transcription kit. Relative transcript levels were determined using SYBR Green PCR master mix on the ABI 7900HT machine. All qRT-PCR assays were conducted in triplicates and normalized to Actb; the comparative CT method was used to analyze the data.

Western blotting

Nuclear protein lysates were prepared using the NE-PER protein extraction kit according to manufacturer’s instructions. Protein immunoblots were done to confirm protein ablation in clones Kdm6aΔE1 and Kdm6aΔE3. Concentrations were determined using the BCA protein assay kit. Equal concentrations of protein lysates were loaded onto SDS-PAGE gels, Tris–glycine 8%. Proteins were transferred to nitrocellulose membranes in Tris–glycine plus methanol buffer overnight at 4 °C. After blocking with Tris-buffered saline (TBS), 0.1% Tween-20, 5% nonfat dry milk for 1 h at room temperature, membranes were incubated at 4 °C overnight with the primary antibody: rabbit polyclonal antibody against KDM6A (1:3000) (Bethyl Laboratories). After primary antibody incubation, immunoblots were incubated with HRP-conjugated secondary antibodies (1:10,000) and visualized by chemiluminescence. Ponceau S staining of total protein served as a loading control.

Sanger sequencing

All plasmid, gDNA and cDNA lysates used for Sanger sequencing were made using Qiagen PCR purification columns. All reactions were prepared and submitted according to the protocols provided by Eurofin Genomics (www.operon.com).

RNA-seq with diploid analysis

Bulk RNA-seq indexed libraries were prepared using Illumina TruSeq RNA sample preparation kit V2. All libraries were generated from 1 µg of total RNA starting material prior to mRNA isolation. For BC derived ES cells libraries were prepared from two biological replicates from wt male controls and three Kdm6aΔ/Y male clones. For CB derived ES cells libraries were prepared from two biological replicates from wt controls and three Kdm6aΔ/Y clones. cDNA libraries were constructed for analysis on a NextSeq sequencer and 75 bp single-end reads were generated. Diploid gene expression was estimated using Tophat/v2.0.14 [29] with default parameters and gene-level expression was normalized using TPM (transcripts per kilobase of exon length per million mapped reads). Differentially expressed genes were determined using DESeq2 analysis [30] with a false discovery rate (FDR) threshold of 0.05 and a fold-change cutoff of 1.5. Gene ontology (GO) analysis was done using Panther Gene Ontology or DAVID GO [31,32,33]. Sex-biased gene expression in BC cells was determined by comparing female and male wt cells. Genes with ≥ twofold TPM expression differences (p ≤ 0.05; student’s t-test) were considered sex-biased. A gain of sex-biased expression was defined as a ≥ twofold TPM expression differences and a p ≤ 0.05 (student’s t-test) following KO. A summary of RNA-seq mapping metrics is shown in Additional file 14: Table S14. For analysis of sex-biased expression in CB cells, data from GSE90516 was re-analyzed to determine sex-biased genes in CB wt female and male cells [34]. Comparisons to CB Kdm6a KO cells were made following a quantile normalization.

RNA-seq with allele-specific analysis

To estimate allele-specific gene expression in hybrid cells, a “pseudo-castaneus” genome was first assembled by substituting known heterozygous SNPs of cast (obtained from Sanger Institute 2015/May, version 5 [35]) into the BL6 mm10 reference genome. RNA-seq reads were mapped using bowtie/v2.2.5 [36] to the genome and transcriptome of both the BL6 and the pseudo-cast genomes. Only those reads that mapped uniquely and with a high-quality mapping score (MAPQ ≥ 30) to either the BL6 or pseudo-cast genomes were kept for further analyses. As previously described, all uniquely mapped and high-quality (MAPQ ≥ 30) reads were segregated into three categories: (1) BL6-SNP reads containing only BL6-specific SNP(s); (2) cast-SNP reads containing only cast-specific SNP(s); (3) allele-uncertain reads, that is, reads that do not contain valid SNPs (Additional file 14: Table S14) [37]. This resulted in ~ 11% and ~ 9% of reads mapped uniquely to the BL6 allele and the cast allele, respectively (Additional file 14: Table S14). For each gene \(i\), denote the number of exonic BL6-SNP reads and cast-SNP reads be \({n}_{i0}\) and \({n}_{i1}\), respectively. The BL6-allele expression proportion for gene \(i\) was then calculated as \({p}_{i0}={n}_{i0}/\left({n}_{i0}+ {n}_{i1}\right)\). To account for the mapping biases between the BL6 and the cast alleles, the BL6-allele expression proportion was further adjusted using the average BL6-allele mapping ratio of all autosomal genes \({r}_{m}={N}_{A0}/{N}_{A1}\), where \({N}_{A0}\) and \({N}_{A1}\) are the number of allele-specific autosomal reads in the BL6 genome and the cast genome, respectively. Specifically, the adjusted BL6-allele expression proportion is \(\widehat{{p}_{i0}}={n}_{i0}/\left({n}_{i0}+ {{r}_{m}n}_{i1}\right)\). We assume that the diploid TPM value of gene \(i\) is the sum of haploid TPM values from the BL6-allele and the cast-allele. Thus, the BL6-allele TPM is estimated as \({\mathrm{TPM}}_{i0}= \widehat{{p}_{i0}}{\mathrm{TPM}}_{i}\), while the cast-allele TPM is estimated as \({\mathrm{TPM}}_{i1}= {\left(1-\widehat{{p}_{i0}}\right)\mathrm{TPM}}_{i}\). Samples derived from the CB cross were analyzed in a similar manner with the BL6- and pseudo-cast reference genomes switched between the maternal and paternal alleles.

Allele-specific differential gene expression was detected using DESeq2 with similar parameters as described in the above diploid DESeq2 analysis [30]. Post filtering of allele-specific expression levels was done to remove genes with estimated allelic expression < 1TPM. Allele-specific DEGs were first identified on the maternal and paternal alleles separately, using the corresponding allelic RNA-seq read counts. Next, allelic downregulated DEGs were categorized into three groups: genes downregulated specifically on the maternal allele (Group A), genes downregulated specifically on the paternal allele (group B), and genes downregulated from both alleles (Group C). Allelic upregulated DEGs were similarly categorized into Groups D, E, and F as described in the main text. Note that our allele-specific RNA-seq analyses were performed independently from the diploid RNA-seq analyses, therefore the resulting allelic DEG list does not completely overlap with the diploid DEG list. For example, one gene could be upregulated from the maternal allele and downregulated from the paternal allele and therefore not identified as a diploid DEG.

ATAC-seq with allele-specific analysis

ATAC-seq was done on wt male ES cells from reciprocal crosses and on clones Kdm6aΔE1 and Kdm6aΔE2.5 as previously described [38]. Data were generated from libraries using a NextSeq sequencer to obtain 75 bp paired-end reads. ATAC-seq reads were mapped to the BL6 genome and the pseudo-cast genome, separately, using bwa/ 0.7.12 [39]. Uniquely mapped reads with a high-quality mapping score (MAPQ ≥ 30) to either the BL6 or the pseudo-cast genome were kept for the following analyses. To assist the allele-specific analysis, ATAC-seq reads were first classified as BL6-SNP reads, cast-SNP reads, or allele-uncertain reads, as described in the RNA-seq analysis. As long as one-end of the paired-end read was allele-certain, then the pair was assigned to the same allele. Read pairs with two ends assigned to different alleles were discarded from further analysis. In the following text, we refer to paired-end reads simply as “reads”. Duplicate reads were removed using the MarkDuplicates function in Picard/v2.14.1. To evaluate DNA accessibility of genes transcribed on the maternal- and paternal-alleles, allelic ATAC-seq read coverage were calculated at promoter regions (2 kb upstream and downstream of TSSs). Quantile normalization was used to minimize allele biases and batch differences between cell lines derived from reciprocal crosses by assuming that the promoter read coverage has a similar distribution between alleles and cell lines.

ChIP-seq with allelic-specific analyses

H3K27me ChIP-seq was performed on wt male ES cells from reciprocal crosses and on clones Kdm6aΔE1 and Kdm6aΔE2.5 as previously described [40]. ChIP libraries were prepared using the TruSeq V2 DNA preparation kit and matching input controls were sequenced by a NextSeq sequencer to produce 75 bp paired-end reads. Similar to ATAC-seq analysis (see above), ChIP-seq reads were mapped to the BL6 genome and the pseudo-cast genome separately, then uniquely mapped reads with high-quality mapping scores were assigned to the correct allele. Allelic promoter (± 2 kb of the TSS) coverage was scaled according to the sequencing depth and then normalized by the input control. To account for potential sequencing depth biases and allele differences, a quantile normalization method was applied (see above).

Results

Kdm6a KO in mouse ES cells results in expression changes in developmentally important genes

Kdm6a was edited using CRISPR/Cas9 in F1 male ES cells derived either from a cross between BL6 x cast (BC) or from the reciprocal cross cast x B6 (CB), in which alleles can be distinguished using SNPs (Fig. 1A and Additional file 1: Fig. S1) [41]. We derived stable Kdm6a-targeted ES cell clones from both BC and CB crosses by inducing either a deletion between exons 2 and 4 (ΔE) or a deletion of exons 1 and 2 in the promoter region (ΔP) (Fig. 1A and Additional file 1: Fig. S1A-D and Additional file 8: Table S8A). We verified editing and loss of protein in seven Kdm6a KO clones, including three male hemizygous clones from the BC cross, hereafter called BC Kdm6aΔ/Y (if all three clones Kdm6aΔE1, Kdm6aΔE3, Kdm6aΔP4 are included), and three male hemizygous clones from the CB cross, hereafter called CB Kdm6aΔ/Y (if all three clones Kdm6aΔE2.5, Kdm6aΔE2.7, Kdm6aΔP2.1 are included) (Fig. 1A and Additional file 1: Fig. S1A-E and Additional file 8: Table S8A). An off-target 314 kb deletion of three adjacent genes (Sh3kbp1, Map3k15, Pdha1) was identified in clones Kdm6aΔE1 and Kdm6aΔE3, consistent with a common origin of these lines (Additional file 1: Fig. S1F). These genes were excluded in DEG analysis when considering expression changes between BC wt clones and BC Kdm6aΔ/Y clones. We also isolated control male ES cell clones (hereafter called CRISPR-controls) derived from the BC cross, which were subject to the CRISPR/Cas9 treatment but did not exhibit a deletion of Kdm6a (Additional file 8: Table S8A).

RNA-seq analysis was done to compare gene expression between wild-type (wt) and Kdm6a KO clones derived from the BC and CB crosses. Using principal component analysis and hierarchal clustering based on expression values for all transcribed genes, wt clones segregated from KO clones (Additional file 2: Fig. S2A, B). The CB lines clustered away from the BC lines, likely due to parental strain-specific differences including the presence of an actin-GFP transgene in the wt CB line available to us. DESeq2 was employed to identify differentially expressed genes (DEGs) in each cross using an FDR cutoff of < 0.05 and a ≥ 1.5 fold-change [30]. Kdm6a KO did not lead to reduced expression of pluripotency and self-renewal genes, confirming that KDM6A is not necessary to maintain pluripotency, nor did the editing process induce ES cell differentiation (Additional file 2: Fig. S2C) [42]. Consistent with previous reports, differentiation kinetics for clone Kdm6aΔE1 compared to wt showed slower formation of embryoid bodies and delayed appearance of neuronal cells after retinoic acid treatment (Additional file 2: Fig. S2D) [43, 44].

In the BC cross, 195 downregulated and 334 upregulated DEGs were identified in the three male KO clones (Kdm6aΔ/Y) compared to the two male wt clones (Fig. 1B; Additional file 9: Table S9A). As expected, downregulated DEGs included known KDM6A targets, which were confirmed by qPCR and RT-PCR (Additional file 2: Fig. S2E, F). There did not appear to be a preference for autosomal genes (~ 4%) or sex-linked genes (~ 5%). After applying our expression fold change threshold and FDR cutoff to a previously published dataset derived from a Kdm6a KO in a male mouse ES line derived from an independent BC cross, we found good concordance with our results for 40% of downregulated and 31% of upregulated DEGs (Additional file 9: Table S9B) [12].

In the CB cross, we identified 334 downregulated and 213 upregulated DEGs in the three male KO Kdm6aΔ/Y clones compared to two wt clones (Fig. 1C; Additional file 9: Table S9C). Again, there did not appear to be a preference for autosomal genes (~ 4%) or sex-linked genes (~ 4%). Surprisingly, the lists of DEGs differed between the CB and BC KO clones, with only 1% downregulated and 6% upregulated DEGs in common (Additional file 2: Fig. S2G and Additional file 9: Table S9A, C). Indeed, of the eight known KDM6A targets confirmed in the BC cross, only Bcar3 and Hsd17b11 were called as DEGs in the CB cross. This low degree of overlap was also observed when published data from a KO line from an independent BC cross was compared to the CB KO clones (11% downregulated and 8% upregulated DEGs) [12]. The discrepancy in DEGs may be due to parental strain differences. The low concordance between the CB and BC KO clones prompted us to use a less stringent approach and simply compare the lists of genes with a ≥ 1.25 fold change in TPM value between male KO and wt clones. Using this approach we identified 345 downregulated and 226 upregulated genes shared between the reciprocal crosses (Fig. 1D; Additional file 9: Table S9D, E). These genes include six known KDM6A targets including Wt1, Bcar3, Foxh1, Dnmt3a, Sox3, and Hsd17b11.

GO terms for downregulated genes common to BC and CB KO clones (≥ 1.25 TPM fold change) are enriched in cardiac ventricle development (e.g. Foxh1 and Adamts6), consistent with heart malformations observed in utero in Kdm6a KO mice (Fig. 1E) [21, 43]. Interestingly, the GO terms bone development (e.g. Has2) and hard palate development (e.g. Cbfb) are also highlighted, consistent with high-arched palate and other distinctive bone and facial features observed in Kabuki type 2 syndrome caused by KDM6A mutations in human (Fig. 1D, E) [45,46,47,48]. There is also decreased expression of Nes, a KDM6A target gene identified in neural crest cells that contribute to patterning and formation of all anterior facial bone and cartilage (Additional file 9: Table S9D) [10]. Some of the downregulated genes in KO clones are associated with regulation of T-cell activation (e.g. Ephb6) and neuron projection development (e.g. Tnik), consistent with autoimmune disorders and neurological issues in Kabuki syndrome (Fig. 1E, F) [48]. Upregulated DEGs in BC and CB KO clones include a subset of genes known to be repressed by KDM6A (e.g. Fam114a1, H1f0, Naglu) (Additional file 9: Table S9E) [13]. Note, we did not observe increased expression of Kdm6a’s Y-linked paralogue, Uty, in any of the Kdm6a KO clones.

Taken together, our findings support a role for KDM6A in the regulation of a subset of key developmentally critical genes in both male and female mouse ES cells. Despite cell line and clonal variability, downregulated DEGs affirm the role of KDM6A as an enhancer of gene expression, whereas upregulated DEGs indicate a repressive function, which may reflect indirect effects on gene expression. Interestingly, some of the dysregulated genes common to the reciprocal crosses are candidates for a role in clinical features of Kabuki type 2 syndrome in human.

KDM6A facilitates sex-biased autosomal gene expression in mouse ES cells

Next, we investigated the role of KDM6A in regulating sex-biased autosomal gene expression. By comparing wt female and male BC ES cells we identified a 2.4 fold greater number of female-biased (2048) versus male-biased genes (868) (≥ twofold TPM difference cut off; p < 0.05) (Fig. 2A; Table 1; Additional file: Fig. S3A and Additional file 10: Table S10A). Interestingly, there was a marked loss of sex-biased expression in BC Kdm6aΔ/Y cells (< twofold TPM difference between female wt and Kdm6aΔ/Y), which amounted to 35% of female-biased genes and 25% of male-biased genes (Fig. 2; Additional file 3: Fig. S3B and Additional file 10: Table S10B). To confirm these results we re-analyzed a published dataset of expression in BC wt cells and found concordance with our data for 767 sex-biased genes (Additional file 10: Table S10C) [34]. Thirty-two percent of these concordant genes showed loss of sex-biased expression following Kdm6a KO, which again affected more female-biased genes than male biased (Additional file 10: Table S10D). GO analysis of genes that lost female-biased expression following KO revealed several biological processes involved in lipid metabolism consistent with the role of KDM6A in adipocytes and metabolic functions (Fig. 2B) [49]. Biological processes for genes with loss of male-biased expression include bone morphogenesis, again consistent with facial features observed in Kabuki type 2 syndrome (Fig. 2B). We also observed some new sex-biased expression after Kdm6a KO, but at a much reduced level, with only 5% genes gaining female-biased expression and 3% genes, male-biased expression (Additional file 3: Fig. S3C, D; Table 1).

Fig. 2
figure 2

Sex-biased gene expression changes in BC Kdm6a KO cells. (A) Number of genes with sex-biased expression (≥ twofold TPM difference cut off; p < 0.05) in BC wt and Kdm6aΔ/Y clones. Overall, there was a greater loss of female-biased genes (pink) than male-biased genes (blue). See also Table 1 and Additional file 10A, B. (B) GO analysis of genes with loss of sex-biased expression in Kdm6aΔ/Y. Numbers in parentheses represent the fold enrichment over expected number of genes within a given biological process, while those not in parentheses represent the number of genes in each category

Table 1 Sex-biased gene expression in wt and Kdm6a KO clones

In the CB cross loss of sex-biased expression after Kdm6a KO was examined using published data for CB male and female ES lines (Table 1; Additional file 3: Fig. S3E and Additional file 10: Fig. S10E) [34]. We found fewer sex-biased genes in wt CB cells and there was a lesser effect of Kdm6a KO on sex-biased gene expression in CB versus BC clones. Nonetheless, Kdm6a KO in CB cells again caused a greater decrease in female- (32%) than male-biased (24%) genes (Table 1; Additional file 3: Fig. S3F and Additional file 10: Fig. S10E). There was little overlap between genes that lost sex-biased expression in BC and CB clones, suggesting strain or clonal effects.

Taken together, these results suggest that KDM6A, either directly or indirectly, plays a regulatory role in perpetuating autosomal sex-biased expression in ES cells with an emphasis on maintenance of sex-biased, especially female-biased gene expression.

The majority of gene expression changes in Kdm6a KO ES cells are allele-specific

To measure allelic gene expression, we assigned RNA-seq reads to parental alleles using SNPs between B6 and cast in the hybrid F1 mouse ES cells from the reciprocal BC and CB crosses. Most genes (94%) had informative SNPs allowing them to be classified as having either bi-allelic expression or parent-of-origin expression bias (see Methods). Allele-specific differential gene expression changes between Kdm6a KO and wt clones were identified using DESeq2 with the same parameters as described above (≥ 1.5 TPM fold change, FDR < 0.05). Overall, scatter plots of differential maternal and paternal expression were similar to those generated using diploid data (Additional file 4: Fig. S4A, B). However, the number of DEGs detected by allelic analysis may differ from those found by diploid analysis (see Methods).

Surprisingly, a majority of downregulated and upregulated DEGs displayed changes limited to one allele. In the BC cross, allelic expression analysis of Kdm6aΔ/Y clones versus wt showed 167/194 (86%) downregulated and 139/189 (74%) upregulated DEGs with allele-specific changes (Table 2). Such changes were equally frequent in Kdm6aΔ/Y clones from the CB cross and affected 237/269 (88%) downregulated and 171/198 (86%) upregulated genes (Table 2). However, similar to the diploid analysis above, there was little overlap between allelic DEGs in CB and BC crosses (Table 2; Additional file 11: Table S11). Downregulated DEGs in KO clones were categorized into three groups: group A includes genes downregulated on the maternal allele, group B, on the paternal allele, and group C, on both alleles (Fig. 3A, C; Table 2). Upregulated DEGs in KO clones were categorized into three groups: group D includes genes upregulated on the maternal allele, group E, on the paternal allele, and group F, on both alleles (Fig. 3B, D; Table 2).

Table 2 Number of allelic DEGs after Kdm6a KO
Fig. 3
figure 3

Allele-specific gene expression changes after Kdm6a KO in male BC and CB clones. A, B Heatmaps of allelic gene expression normalized to the mean of each gene across alleles show the extent of fold changes for maternal (purple) and paternal (dark blue) alleles of downregulated DEGs (DESeq2—≥ 1.5 fold change and FDR < 0.05) (A) and upregulated DEGs (B) in two male wt controls, three male BC KO clones Kdm6aΔ/Y including two clones Kdm6aΔ/E with an exon 2–4 deletion (Kdm6aΔE1 Kdm6aΔE3) and one clone with a promoter deletion (Kdm6aΔP4). The genotypes of the clones are color-coded at the top. X-linked gene expression from the paternal allele was masked as zero (white boxes). Z scores (color-coded) representing deviations from the mean for each row are shown for each group of genes. C, D Same analysis as in (A, B) for but for CB clones. See also Table 2 and Additional file 11: Table S11. Group A includes genes downregulated on the maternal allele, group B, on the paternal allele, and group C, on both alleles. Group D includes genes upregulated on the maternal allele, group E, on the paternal allele, and group F, on both alleles

Remarkably, group A genes show a preferential downregulation of maternal alleles in both the BC and CB crosses. Considering the BC cross, group A comprises 109 genes, whereas group B consists of only of 58 genes (Fig. 3A; Table 2; Additional file 11: Table S11). When only considering clones with exonic Kdm6a deletion (Kdm6aΔE1 and Kdm6aΔE3) overlapping genes continue to show preferential downregulation of maternal alleles (38 for group A versus 20 for group B). Further, comparisons between individual BC KO clones show highly similar expression changes for all allelicly regulated groups (Additional file 4: Fig. S4C). In CB KO clones, there were again more group A (133) than group B genes (104). Similar to BC cells, comparisons between individual CB KO clones show highly similar expression changes for all allelicly regulated groups (Additional file 4: Fig. S4D). Despite differences between BC and CB crosses in A and B genes identified as allelic DEGs, additional analyses of changes in TPM values reveal that 33% of group A and 50% of group B genes identified in male BC clones showed a trend of decreased expression in male CB clones (Fig. 3C; Table 2; Additional file 11: Table S11 and Additional file 12: Table S12). Among downregulated genes in BC clones we found a subset of maternally (but not paternally) expressed imprinted genes that include the Meg3 cluster, Xlr cluster, and Phlda2 (Additional file 5: Fig. S5A-E and Additional file 11: Table S11A). Another maternally expressed imprinted gene H19 that regulates Igf2 showed a significant decrease in clones Kdm6aΔE1 and Kdm6aΔE3 (p = 0.01) and a trend for a decrease in clone Kdm6aΔP4, as well as a corresponding increase at Igf2 (Additional file 5: Fig. S5F, G). Changes in maternally expressed imprinted genes were not verified in the CB cross where there were few changes in imprinted gene expression, suggesting clonal or parental strain effects (Additional file 5: Fig. S5H). Further comparisons between BC and CB crosses were made to address allelic expression bias of group A and B genes in wt cells. We defined allelic bias as ≥ twofold difference between alleles. We found that a majority of group A genes (maternally downregulated genes) identified in male BC clones (65%) and CB clones (57%) exhibited maternally biased expression in wt clones (Fig. 3A, C; Additional file 5: Fig. S5I, J). This maternal allele biased expression of group A genes in wt cells was also reported in an independently derived CB ES cell line [50]. Group B genes (paternally downregulated) showed a lower percentage of paternally biased expression in wt (BC—36%, CB—15%). Similar to our observations using diploid analyses, there was little overlap in the identity of genes with allelic expression biases in BC and CB wt cells, suggesting parental strain effects (Additional file 5: Fig. S5I, J).

In contrast to downregulated genes, upregulated genes did not display an allelic bias, as reflected by a similar number of genes in groups D (68) and E (71) in the male BC clones, but did display a maternal allele bias when comparing groups D (109) and E (69) in the male CB clones (Fig. 3B, D; Table 2; Additional file 11: Table S11G, H).

We conclude that KDM6A preferentially facilitates expression of maternal alleles of a subset of mouse genes. While this preference is seen in Kdm6a KO cells from both BC and CB crosses, the lack of overlap between the lists of genes affected in the two crosses suggests that KDM6A regulation of maternal alleles is not gene-specific. Interestingly, a majority of these genes exhibit maternal allele-biased expression in wt ES cells, suggesting that differential Kdm6a expression in male and female germ cells may influence parental allele expression.

Epigenetic changes induced by Kdm6a KO

Epigenetic changes and chromatin accessibility were examined at promoter regions (± 2 kb from the TSS) of genes with expression changes after Kdm6a KO, using H3K27me3 ChIP-seq and ATAC-seq to compare male BC Kdm6aΔE1 and CB Kdm6aΔE2.5 clones with their respective wt controls. As an initial control, we verified the expected higher H3K27me3 level at the downregulated Hoxb cluster in Kdm6a KO clones (Fig. 4A) [4, 10, 22]. Consistent with diploid gene expression changes in BC and CB clones, non-allelic analyses show most promoters of downregulated genes present with increased levels of H3K27me3 and decreased chromatin accessibility, while promoters of upregulated genes show decreased H3K27me3 and increased chromatin accessibility (Fig. 4B-D). To rule out indirect effects we verified that expression levels of genes encoding PRC2 complex-associated proteins that mediate H3K27 methylation (Ezh1/2, Suz12, Eed, and Rbbp7) were unchanged in KO clones (Additional file 6: Fig. S6A). Similarly, expression levels of genes (Ep300, Smarca4, Eomes) encoding chromatin remodeling proteins known to mediate H3K27 acetylation by association with KDM6A in the MLL complex were unchanged (Additional file 6: Fig. S6A). However, this analysis does not exclude changes in recruitment of these proteins.

Fig. 4
figure 4

Diploid changes in H3K27me3 and chromatin accessibility after Kdm6a KO in BC and CB clones. A UCSC Genome browser (GRCm38/mm10) view of H3K27me3 profiles in wt and Kdm6aΔE1 confirms a general increase across the entire Hoxb locus following Kdm6a KO. wt (black); KO cells (orange). B Box plots of average ratios of H3K27me3 read coverage in BC (Kdm6aΔE1) and CB KO clones (Kdm6aΔE2.5) versus wt at promoters (± 2 kb of the TSS) of genes with similar changes between the reciprocal crosses. Downregulated and upregulated genes exhibit expected increases and decreases in H3K27me3 (***p ≤ 0.0001). C Same analysis as in (B), but for ATAC-seq average ratios of read coverage. Downregulated and upregulated genes exhibit the expected decreases and increases in chromatin accessibility (***p ≤ 0.0001; **p = 0.011). D Example profiles of H3K27me3 read coverage in wt (black) and KO cells (orange) show an increase at promoters of downregulated genes (Tnik, Has2, Ephb6) in BC and CB KO clones. (UCSC Genome browser GRCm38/mm10)

Next, we performed allele-specific ChIP-seq and ATAC-seq to detect H3K27me3 and chromatin accessibility changes on each parental allele. In BC KO cells a significant increase in H3K27me3 at the promoters of maternally downregulated genes (group A) was accompanied by a significant decrease in chromatin accessibility, while maternally upregulated group D and E genes showed the expected increase in chromatin accessibility (Fig. 5A-C; Additional file 7: Fig. S7A and Additional file 13: Table S13). In contrast, paternally downregulated alleles (group B) showed an increase in H3K27me3 but no significant decrease in chromatin accessibility, while paternally upregulated (group E) genes showed the expected increase in chromatin accessibility (Fig. 5D-F; Additional file 7: Fig. S7A and Additional file 13: Table S13). Importantly, in clones with exonic Kdm6a deletion (Kdm6aΔE1 and Kdm6aΔE3), promoters of overlapping genes show similar changes in H3K27me3 and chromatin accessibility (Figs. 3A, B; Additional file 7: Fig. S7B, D, F, H). In the CB cross similar findings for H3K27me3 were obtained for group A genes, while upregulated groups D and E genes showed only a modest increase in chromatin accessibility (Fig. 5G, H; Additional file 13: Table S13). Of note, differences in chromatin accessibility were only observed at promoters of KDM6A target genes, and there were no overall differences in accessibility at other promoters or enhancers (Additional file 6: Fig. S6B, C). Additionally, H3K27me3 showed no overall differences at enhancers regardless of parental allele (Additional file 6: Fig. S6).

Fig. 5
figure 5

Allelic changes in H3K27me3 and chromatin accessibility after Kdm6a KO in BC and CB clones. A, B Box plots show average ratios of H3K27me3 (A) and ATAC-seq (B) read coverage between clone Kdm6aΔE1 and wt cells at promoters (± 2 kb of the TSS) on maternal alleles. Both downregulated groups (A and B) show a significant increase in H3K27me3 levels on maternal alleles (A **p = 0.001; B ***p = 5.33e−06). Maternally downregulated genes (group A) show significant loss of chromatin accessibility (***p = 5.9e−11) on the maternal allele in KO cells, while paternally downregulated genes (group B) do not. Upregulated groups (D and E) show no or slight changes in H3K27me3 (E *p = 0.020), and show the expected increase in chromatin accessibility (D ***p = 2.4e−07; E *p = 0.03). Values were normalized to the input and converted to log scale. C UCSC Genome browser (GRCm38/mm10) view of maternal H3K27me3 and chromatin accessibility profiles in wt (−) and Kdm6aΔE1 ( +) at the maternally downregulated gene Exoc8 shows an increase in H3K27me3 and a corresponding decrease in ATAC-seq around the promoter. D, E Same analysis as in (A, B) but for paternal alleles. Both downregulated groups (A and B) show an increase in H3K27me3 levels (A ***p = 4.32e−07; B ***p = 1.73e−06), while upregulated groups (D and E) show a slight but significant increase (D *p = 0.03; E *p = 0.04). Only maternally downregulated group A and paternally upregulated group E show significant changes in ATAC-seq (A ***p = 0.0001; B ***p = 0.001) at paternal alleles. F UCSC Genome browser (GRCm38/mm10) view of paternal H3K27me3 and chromatin accessibility profiles in wt (−) and Kdm6aΔE1 ( +) at the paternally downregulated Neurog3 gene shows increased H3K27me3 but no corresponding decrease in ATAC-seq around the promoter. G Same analysis as in (A, B) but in the CB Kdm6aΔE2.5 clone. Only maternally downregulated genes (group A) show increased H3K27me3. There were no significant ATAC-seq changes. H Same analysis as in (D, E) but in the CB Kdm6aΔE2.5 clone. Maternally downregulated group A show increased H3K27me3, while ATAC-seq changes were only seen on paternal alleles. Group A includes genes downregulated on the maternal allele and group B on the paternal allele. Group D includes genes upregulated on the maternal allele and group E, on the paternal allele

Taken together our results show the expected changes in H3K27me3 and chromatin accessibility at upregulated and downregulated DEGs. Allelic analysis demonstrates parental allele-specific changes, notably for downregulated group A genes, which display a consistent increase in H3K27me3 on both parental alleles in both crosses. However, chromatin accessibility changes were partly inconsistent with changes in H3K27me3 in BC and CB KO clones, a phenomenon observed in other studies [10, 51, 52].

Discussion

Here, we show that Kdm6a KO in ES cells derived from reciprocal crosses between C57BL/6 J and Mus castaneus results in significant dysregulation of genes involved in development and reproduction. Depletion of KDM6A leads to a loss in sex-biased gene expression, particularly female-biased expression, suggesting a new role for the histone demethylase. We also found a strong downregulation of maternal alleles in Kdm6a KO clones from both reciprocal crosses. Furthermore, a large subset of the target genes exhibits maternal allele expression bias in wt ES cells from both crosses. The higher expression of Kdm6a in females versus males due to escape from X inactivation may explain the preferential regulation of maternal allele expression.

Gene regulation by KDM6A mainly depends on its catalytic function. Indeed, depletion or inhibition of KDM6A results in decreased gene expression due to higher levels of the repressive histone mark H3K27me3 at promoters and/or enhancers of target genes, which readily explains our findings of H3K27me3 enrichment and downregulation of genes after Kdm6a KO [14, 42, 53]. Conversely, higher KDM6A expression in human cancer causes a decrease of H3K27me3 enriched heterochromatin [54]. Aside from its catalytic function, KDM6A can also enhance gene expression through other mechanisms that involve recruitment of the MLL complex to increase histone acetylation [10,11,12]. Finally, KDM6A may act as a repressor of gene expression, underscoring the complexity of this protein [10, 11, 13, 55]. Indirect effects are also expected, which could explain some of our findings of upregulated genes after Kdm6a KO. In particular, the more than twofold reduction in Dnmt3a expression in KO male BC clones could cause DNA hypomethylation and upregulation of DNMT3A target genes, such as Celsr1 and Fah [56]. Similarly, the decrease in Tet1 expression in KO female ES cells could explain increased expression of TET1 target genes such as Casp8.

A number of KO dysregulated genes were found to differ between crosses, which could be due to differences between the wt ES cell lines and/or to strain differences. Clonal variability is a limitation of the study, which we addressed by using both a stringent DESeq2 approach and a less stringent fold change approach to identify gene expression changes in BC and CB KO cells. Importantly, among the KDM6A targets common to both reciprocal crosses we identified a subset of genes associated with multiple phenotypes observed in individuals with Kabuki syndrome due to mutations in KDM6A in human. For example, Has2 and Cbfb have been implicated in hard palate development, malformations of which are often seen in Kabuki syndrome [48, 57, 58]. Ephb6 and Tnik were also downregulated in Kdm6a KO cells derived from reciprocal crosses consistent with immunopathological phenotypes and neurological disorders seen in the syndrome [48]. Our identification of KO dysregulated genes involved in development, immune cell functions, neuronal cell functions, and metabolism may help further identify KDM6A targets implicated in Kabuki syndrome, a disorder with more severe phenotypes in males [48, 59, 60].

We report on a new role for KDM6A, which is to maintain sex-biased, especially female-biased gene expression. KDM6A appears to maintain female-biased expression via repression of genes in males, which is surprising since KDM6A’s main function is to facilitate gene expression through removal of H3K27me3. However, this may be due to indirect effects of KDM6A depletion as mentioned above. Indeed, there appears to be little commonality in genes that gain sex bias following KO.

A surprising finding from our allelic expression analyses is that a majority of gene dysregulation after Kdm6a KO occurs on one allele or the other, with a strong maternal allele bias. This novel observation was made in both reciprocal crosses, even though the sets of dysregulated alleles differ. Female-biased expression of Kdm6a due to escape from X inactivation is consistent with a role in differential regulation of the maternal and paternal genomes [61]. How this allelic regulation is transmitted through generations will need to be further investigated. H3K27me3 is known to play a critical role in the repression of the maternal allele at non-canonical imprinted sites, loss of which leads to loss of imprinting [62]. Thus, manipulating the dosage of Kdm6a in differentiating hybrid primordial germ cells where genomic imprints are established may yield insights into its role in regulating non-canonical imprinted genes. Further, utilizing independent reciprocal crosses with Kdm6a KO will help to identify potential non-canonical imprinted gene targets. Future allelic expression studies using differentiating ES cells following KO would be interesting since KDM6A target genes may also have roles in germ layer development and ES cell differentiation [14, 17]. We found changes in H3K27me3 and chromatin accessibility at dysregulated alleles but these were not always concordant with expression changes. This was not entirely unexpected as H3K27me3 levels are not always correlated with expression [10, 52]. Furthermore, many of allelicly repressed genes were more highly expressed from downregulated alleles in wt suggesting a higher sensitivity to changes in H3K27me3.

Perspectives and significance

Results from our study provide important findings implicating KDM6A, a female-biased chromatin modifier, in allelic gene regulation with a preference for maternal alleles. These data emphasize that inclusion of allelic analyses of gene regulation will enhance our understanding of parent-of-origin effects and sex-specific differences. Our study also provides evidence of a role of KDM6A in the regulation of sex-biased gene expression, which could have broad implications in understanding sex differences. Our CRISPR deletions encompass known pathogenic mutations seen in Kabuki syndrome, including loss of exons 1–2 [63]. Thus, results from this study may reveal therapeutic targets for individuals with Kabuki syndrome.

Considering KDM6A’s role in embryonic development, it will be of interest to determine whether sex-biased and parental allele effects of depletion persist during differentiation. Kdm6a KO mice have a sex-specific phenotype, with male mice runt but able to survive, while female mice die during embryogenesis [21]. This has been attributed to a demethylase-independent mechanism of compensation by UTY, which has little or no demethylase activity [21]. A number of escape genes have a Y-linked paralogues that may have synergistic roles in gene regulation. Therefore, it will be important to fully tease out the possible compensatory role of Uty on gene expression by generating ES cell clones with Kdm6a and Uty specific KO. Future analyses considering additional epigenetic marks such as H3K4me3, H3K9ac and DNA methylation will help to ascertain the complex network of gene regulation by KDM6A. Finally, extending our studies of KDM6A’s role in gene regulation to human cells will help determine whether KDM6A’s functions are conserved between mice and humans at target genes and regulatory elements.

Conclusions

Here, we demonstrate for the first time a role in the regulation of sex-biased genes and a maternal allele bias in gene regulation by a gene that escapes X inactivation, a uniquely female phenomenon. Our results implicate KDM6A as a master regulator of fetal growth and development through the control of a subset of maternally biased genes. While our studies are performed in ES cells where developmental pathways are not fully represented, our data provide information on pathways and lineages that may be impacted in syndromes due to KDM6A mutations such as Kabuki Syndrome.

Availability of data and materials

The datasets generated and/or analyzed during the current study are available in the NCBI GEO database. Oligonucleotide sequences used in this study are available in Additional file 1.

References

  1. De Kumar B, Parrish ME, Slaughter BD, Unruh JR, Gogol M, Seidel C, et al. Analysis of dynamic changes in retinoid-induced transcription and epigenetic profiles of murine Hox clusters in ES cells. Genome Res. 2015;25(8):1229–43. https://doi.org/10.1101/gr.184978.114.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Dhar SS, Lee SH, Chen K, Zhu G, Oh W, Allton K, et al. An essential role for UTX in resolution and activation of bivalent promoters. Nucleic Acids Res. 2016;44(8):3659–74. https://doi.org/10.1093/nar/gkv1516.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Agger K, Cloos PA, Christensen J, Pasini D, Rose S, Rappsilber J, et al. UTX and JMJD3 are histone H3K27 demethylases involved in HOX gene regulation and development. Nature. 2007;449(7163):731–4. https://doi.org/10.1038/nature06145.

    Article  CAS  PubMed  Google Scholar 

  4. Yu J, Wang L, Pei P, Li X, Wu J, Qiu Z, et al. Reduced H3K27me3 leads to abnormal Hox gene expression in neural tube defects. Epigenet Chromatin. 2019;12(1):76. https://doi.org/10.1186/s13072-019-0318-1.

    Article  CAS  Google Scholar 

  5. Seenundun S, Rampalli S, Liu QC, Aziz A, Palii C, Hong S, et al. UTX mediates demethylation of H3K27me3 at muscle-specific genes during myogenesis. EMBO J. 2010;29(8):1401–11. https://doi.org/10.1038/emboj.2010.37.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Bosselut R. Pleiotropic functions of H3K27Me3 demethylases in immune cell differentiation. Trends Immunol. 2016;37(2):102–13. https://doi.org/10.1016/j.it.2015.12.004.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Lei X, Jiao J. UTX affects neural stem cell proliferation and differentiation through PTEN signaling. Stem Cell Rep. 2018;10(4):1193–207. https://doi.org/10.1016/j.stemcr.2018.02.008.

    Article  CAS  Google Scholar 

  8. Itoh Y, Golden LC, Itoh N, Matsukawa MA, Ren E, Tse V, et al. The X-linked histone demethylase Kdm6a in CD4+ T lymphocytes modulates autoimmunity. J Clin Invest. 2019;130:3852–63. https://doi.org/10.1172/JCI126250.

    Article  Google Scholar 

  9. Shan Y, Zhang Y, Zhao Y, Wang T, Zhang J, Yao J, et al. JMJD3 and UTX determine fidelity and lineage specification of human neural progenitor cells. Nat Commun. 2020;11(1):382. https://doi.org/10.1038/s41467-019-14028-x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Shpargel KB, Starmer J, Wang C, Ge K, Magnuson T. UTX-guided neural crest function underlies craniofacial features of Kabuki syndrome. Proc Natl Acad Sci U S A. 2017;114(43):E9046–55. https://doi.org/10.1073/pnas.1705011114.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Miller SA, Mohn SE, Weinmann AS. Jmjd3 and UTX play a demethylase-independent role in chromatin remodeling to regulate T-box family member-dependent gene expression. Mol Cell. 2010;40(4):594–605. https://doi.org/10.1016/j.molcel.2010.10.028.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Wang SP, Tang Z, Chen CW, Shimada M, Koche RP, Wang LH, et al. A UTX-MLL4-p300 transcriptional regulatory network coordinately shapes active enhancer landscapes for eliciting transcription. Mol Cell. 2017;67(2):308-21 e6. https://doi.org/10.1016/j.molcel.2017.06.028.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Gozdecka M, Meduri E, Mazan M, Tzelepis K, Dudek M, Knights AJ, et al. UTX-mediated enhancer and chromatin remodeling suppresses myeloid leukemogenesis through noncatalytic inverse regulation of ETS and GATA programs. Nat Genet. 2018;50(6):883–94. https://doi.org/10.1038/s41588-018-0114-z.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Berletch JB, Deng X, Nguyen DK, Disteche CM. Female bias in Rhox6 and 9 regulation by the histone demethylase KDM6A. PLoS Genet. 2013;9(5): e1003489. https://doi.org/10.1371/journal.pgen.1003489.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Greenfield A, Carrel L, Pennisi D, Philippe C, Quaderi N, Siggers P, et al. The UTX gene escapes X inactivation in mice and humans. Hum Mol Genet. 1998;7(4):737–42. https://doi.org/10.1093/hmg/7.4.737.

    Article  CAS  PubMed  Google Scholar 

  16. Tukiainen T, Villani AC, Yen A, Rivas MA, Marshall JL, Satija R, et al. Landscape of X chromosome inactivation across human tissues. Nature. 2017;550(7675):244–8. https://doi.org/10.1038/nature24265.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Tran N, Broun A, Ge K. Lysine demethylase KDM6A in differentiation, development, and cancer. Mol Cell Biol. 2020;40:20. https://doi.org/10.1128/MCB.00341-20.

    Article  Google Scholar 

  18. Davis EJ, Broestl L, Abdulai-Saiku S, Worden K, Bonham LW, Minones-Moyano E, et al. A second X chromosome contributes to resilience in a mouse model of Alzheimer’s disease. Sci Transl Med. 2020;12:558. https://doi.org/10.1126/scitranslmed.aaz5677.

    Article  CAS  Google Scholar 

  19. Bellott DW, Hughes JF, Skaletsky H, Brown LG, Pyntikova T, Cho TJ, et al. Mammalian Y chromosomes retain widely expressed dosage-sensitive regulators. Nature. 2014;508(7497):494–9. https://doi.org/10.1038/nature13206.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Cortez D, Marin R, Toledo-Flores D, Froidevaux L, Liechti A, Waters PD, et al. Origins and functional evolution of Y chromosomes across mammals. Nature. 2014;508(7497):488–93. https://doi.org/10.1038/nature13151.

    Article  CAS  PubMed  Google Scholar 

  21. Shpargel KB, Sengoku T, Yokoyama S, Magnuson T. UTX and UTY demonstrate histone demethylase-independent function in mouse embryonic development. PLoS Genet. 2012;8(9): e1002964. https://doi.org/10.1371/journal.pgen.1002964.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Shpargel KB, Starmer J, Yee D, Pohlers M, Magnuson T. KDM6 demethylase independent loss of histone H3 lysine 27 trimethylation during early embryonic development. PLoS Genet. 2014;10(8): e1004507. https://doi.org/10.1371/journal.pgen.1004507.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. von Meyenn F, Iurlaro M, Habibi E, Liu NQ, Salehzadeh-Yazdi A, Santos F, et al. Impairment of DNA methylation maintenance is the main cause of global demethylation in naive embryonic stem cells. Mol Cell. 2016;62(6):848–61. https://doi.org/10.1016/j.molcel.2016.04.025.

    Article  CAS  Google Scholar 

  24. Yagi M, Kishigami S, Tanaka A, Semi K, Mizutani E, Wakayama S, et al. Derivation of ground-state female ES cells maintaining gamete-derived DNA methylation. Nature. 2017;548(7666):224–7. https://doi.org/10.1038/nature23286.

    Article  CAS  PubMed  Google Scholar 

  25. Byrne SM, Ortiz L, Mali P, Aach J, Church GM. Multi-kilobase homozygous targeted gene replacement in human induced pluripotent stem cells. Nucleic Acids Res. 2015;43(3): e21. https://doi.org/10.1093/nar/gku1246.

    Article  CAS  PubMed  Google Scholar 

  26. Montague TG, Cruz JM, Gagnon JA, Church GM, Valen E. CHOPCHOP: a CRISPR/Cas9 and TALEN web tool for genome editing. Nucleic Acids Res. 2014;42(Web Server issue):W401–7. https://doi.org/10.1093/nar/gku410.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Labun K, Montague TG, Gagnon JA, Thyme SB, Valen E. CHOPCHOP v2: a web tool for the next generation of CRISPR genome engineering. Nucleic Acids Res. 2016;44(W1):W272–6. https://doi.org/10.1093/nar/gkw398.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Cong L, Ran FA, Cox D, Lin S, Barretto R, Habib N, et al. Multiplex genome engineering using CRISPR/Cas systems. Science. 2013;339(6121):819–23. https://doi.org/10.1126/science.1231143.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36. https://doi.org/10.1186/gb-2013-14-4-r36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550. https://doi.org/10.1186/s13059-014-0550-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Mi H, Muruganujan A, Ebert D, Huang X, Thomas PD. PANTHER version 14: more genomes, a new PANTHER GO-slim and improvements in enrichment analysis tools. Nucleic Acids Res. 2019;47(D1):D419–26. https://doi.org/10.1093/nar/gky1038.

    Article  CAS  PubMed  Google Scholar 

  32. da Huang W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57. https://doi.org/10.1038/nprot.2008.211.

    Article  CAS  Google Scholar 

  33. da Huang W, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37(1):1–13. https://doi.org/10.1093/nar/gkn923.

    Article  CAS  Google Scholar 

  34. Werner RJ, Schultz BM, Huhn JM, Jelinek J, Madzo J, Engel N. Sex chromosomes drive gene expression and regulatory dimorphisms in mouse embryonic stem cells. Biol Sex Differ. 2017;8(1):28. https://doi.org/10.1186/s13293-017-0150-x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Keane TM, Goodstadt L, Danecek P, White MA, Wong K, Yalcin B, et al. Mouse genomic variation and its effect on phenotypes and gene regulation. Nature. 2011;477(7364):289–94. https://doi.org/10.1038/nature10413.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9(4):357–9. https://doi.org/10.1038/nmeth.1923.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Berletch JB, Ma W, Yang F, Shendure J, Noble WS, Disteche CM, et al. Escape from x inactivation varies in mouse tissues. PLoS Genet. 2015;11(3): e1005079. https://doi.org/10.1371/journal.pgen.1005079.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Buenrostro JD, Giresi PG, Zaba LC, Chang HY, Greenleaf WJ. Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position. Nat Methods. 2013;10(12):1213–8. https://doi.org/10.1038/nmeth.2688.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25(14):1754–60. https://doi.org/10.1093/bioinformatics/btp324.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Yang F, Deng X, Ma W, Berletch JB, Rabaia N, Wei G, et al. The lncRNA Firre anchors the inactive X chromosome to the nucleolus by binding CTCF and maintains H3K27me3 methylation. Genome Biol. 2015;16:52. https://doi.org/10.1186/s13059-015-0618-0.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Barakat TS, Rentmeester E, Sleutels F, Grootegoed JA, Gribnau J. Precise BAC targeting of genetically polymorphic mouse ES cells. Nucleic Acids Res. 2011;39(18): e121. https://doi.org/10.1093/nar/gkr550.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Mansour AA, Gafni O, Weinberger L, Zviran A, Ayyash M, Rais Y, et al. The H3K27 demethylase Utx regulates somatic and germ cell epigenetic reprogramming. Nature. 2012;488(7411):409–13. https://doi.org/10.1038/nature11272.

    Article  CAS  PubMed  Google Scholar 

  43. Wang C, Lee JE, Cho YW, Xiao Y, Jin Q, Liu C, et al. UTX regulates mesoderm differentiation of embryonic stem cells independent of H3K27 demethylase activity. Proc Natl Acad Sci U S A. 2012;109(38):15324–9. https://doi.org/10.1073/pnas.1204166109.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Morales Torres C, Laugesen A, Helin K. Utx is required for proper induction of ectoderm and mesoderm during differentiation of embryonic stem cells. PLoS ONE. 2013;8(4): e60020. https://doi.org/10.1371/journal.pone.0060020.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Guo Z, Liu F, Li HJ. Novel KDM6A splice-site mutation in kabuki syndrome with congenital hydrocephalus: a case report. BMC Med Genet. 2018;19(1):206. https://doi.org/10.1186/s12881-018-0724-4.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Adam MP, Banka S, Bjornsson HT, Bodamer O, Chudley AE, Harris J, et al. Kabuki syndrome: international consensus diagnostic criteria. J Med Genet. 2019;56(2):89–95. https://doi.org/10.1136/jmedgenet-2018-105625.

    Article  PubMed  Google Scholar 

  47. Yap CS, Shekhar Jamuar S, Lai AHM, Tan ES, Ng I, Ting TW, et al. Identification of KMT2D and KDM6A variants by targeted sequencing from patients with Kabuki syndrome and other congenital disorders. Gene. 2020. https://doi.org/10.1016/j.gene.2020.144360.

    Article  PubMed  Google Scholar 

  48. Boniel S, Szymanska K, Smigiel R, Szczaluba K. Kabuki syndrome—clinical review with molecular aspects. Genes (Basel). 2021;12:4. https://doi.org/10.3390/genes12040468.

    Article  CAS  Google Scholar 

  49. Li F, Wang S, Cui X, Jing J, Yu L, Xue B, et al. Adipocyte Utx deficiency promotes high-fat diet-induced metabolic dysfunction in mice. Cells. 2022;11:2. https://doi.org/10.3390/cells11020181.

    Article  CAS  Google Scholar 

  50. Marks H, Kerstens HH, Barakat TS, Splinter E, Dirks RA, van Mierlo G, et al. Dynamics of gene silencing during X inactivation using allele-specific RNA-seq. Genome Biol. 2015;16:149. https://doi.org/10.1186/s13059-015-0698-x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Gentile C, Berlivet S, Mayran A, Paquette D, Guerard-Millet F, Bajon E, et al. PRC2-associated chromatin contacts in the developing limb reveal a possible mechanism for the atypical role of PRC2 in HoxA gene expression. Dev Cell. 2019;50(2):184-96 e4. https://doi.org/10.1016/j.devcel.2019.05.021.

    Article  CAS  PubMed  Google Scholar 

  52. Young MD, Willson TA, Wakefield MJ, Trounson E, Hilton DJ, Blewitt ME, et al. ChIP-seq analysis reveals distinct H3K27me3 profiles that correlate with transcriptional activity. Nucleic Acids Res. 2011;39(17):7415–27. https://doi.org/10.1093/nar/gkr416.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Faralli H, Wang C, Nakka K, Benyoucef A, Sebastian S, Zhuang L, et al. UTX demethylase activity is required for satellite cell-mediated muscle regeneration. J Clin Invest. 2016;126(4):1555–65. https://doi.org/10.1172/JCI83239.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Gurrion C, Uriostegui M, Zurita M. Heterochromatin reduction correlates with the increase of the KDM4B and KDM6A demethylases and the expression of pericentromeric DNA during the acquisition of a transformed phenotype. J Cancer. 2017;8(14):2866–75. https://doi.org/10.7150/jca.19477.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Lee S, Lee JW, Lee SK. UTX, a histone H3-lysine 27 demethylase, acts as a critical switch to activate the cardiac developmental program. Dev Cell. 2012;22(1):25–37. https://doi.org/10.1016/j.devcel.2011.11.009.

    Article  CAS  PubMed  Google Scholar 

  56. Haney SL, Upchurch GM, Opavska J, Klinkebiel D, Hlady RA, Roy S, et al. Dnmt3a is a haploinsufficient tumor suppressor in CD8+ peripheral T cell lymphoma. PLoS Genet. 2016;12(9): e1006334. https://doi.org/10.1371/journal.pgen.1006334.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Lan Y, Qin C, Jiang R. Requirement of hyaluronan synthase-2 in craniofacial and palate development. J Dent Res. 2019;98(12):1367–75. https://doi.org/10.1177/0022034519872478.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Sarper SE, Inubushi T, Kurosaka H, Ono Minagi H, Murata Y, Kuremoto KI, et al. Anterior cleft palate due to Cbfb deficiency and its rescue by folic acid. Dis Model Mech. 2019;12:6. https://doi.org/10.1242/dmm.038851.

    Article  CAS  Google Scholar 

  59. Lintas C, Persico AM. Unraveling molecular pathways shared by Kabuki and Kabuki-like syndromes. Clin Genet. 2018;94(3–4):283–95. https://doi.org/10.1111/cge.12983.

    Article  CAS  PubMed  Google Scholar 

  60. Miyake N, Mizuno S, Okamoto N, Ohashi H, Shiina M, Ogata K, et al. KDM6A point mutations cause Kabuki syndrome. Hum Mutat. 2013;34(1):108–10. https://doi.org/10.1002/humu.22229.

    Article  CAS  PubMed  Google Scholar 

  61. Xu J, Deng X, Watkins R, Disteche CM. Sex-specific differences in expression of histone demethylases Utx and Uty in mouse brain and neurons. J Neurosci. 2008;28(17):4521–7. https://doi.org/10.1523/JNEUROSCI.5382-07.2008.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Inoue A, Chen Z, Yin Q, Zhang Y. Maternal Eed knockout causes loss of H3K27me3 imprinting and random X inactivation in the extraembryonic cells. Genes Dev. 2018;32(23–24):1525–36. https://doi.org/10.1101/gad.318675.118.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Faundes V, Goh S, Akilapa R, Bezuidenhout H, Bjornsson HT, Bradley L, et al. Clinical delineation, sex differences, and genotype-phenotype correlation in pathogenic KDM6A variants causing X-linked Kabuki syndrome type 2. Genet Med. 2021;23(7):1202–10. https://doi.org/10.1038/s41436-021-01119-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank Drs. Di Kim Nguyen and Xinxian Deng (University of Washington) for consultation on this study and for critical reading of the manuscript and Charlie Lee (University of Washington) for his help with next-generation sequencing. Hybrid mouse ES cells used in this study were graciously donated by Joost Gribnau (Erasmus MC, University Medical Center).

Funding

This work was supported by NIH grants GM131745 (CMD), MH105768 (JB), and NSF grant DBI-1751317 (WM). The funders had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information

Authors and Affiliations

Authors

Contributions

Conceptualization, J.B. and C.M.D.; Investigation, J.B., H.F., G.N.F, N.P.; Formal analysis, W.M.; Data Curation, W.F.; Writing—Original Draft, J.B. and C.M.D.; Writing—Review and Editing, J.B., W.M., G.N.F., and C.M.D.; Visualization, J.B., W.M., and C.M.D.; Supervision, J.B. and C.M.D.; Funding Acquisition, J.B. and C.M.D. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Wenxiu Ma or Joel B. Berletch.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

All authors have given their consent for publication of the manuscript in Biology of Sex Differences.

Competing interests

The authors declare no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Fig. S1.

CRISPR/Cas9 editing of Kdm6a in BC and CB mouse ES cells. (A) Schematic shows location of the exonic deletion (Kdm6aΔE) that removed exons 2–4 of Kdm6a in male BC and CB ES cells. Exons are shown as vertical bars with location of the PCR and RT-PCR primers (color-coded arrows) used to confirm the deletion and measure expression. Above, images of gels after electrophoresis of PCR products (gDNA), RT-PCR products (cDNA), and Western blot (protein) confirm KDM6A KO in male BC clones (Kdm6aΔE1, Kdm6aΔE3) and CB clones (Kdm6aΔE2.5, Kdm6aΔE2.7) compared to wt (color-coding refers to primers shown on schematic). Actb was run as a control for PCR and RT-PCR and Ponceau S staining was used as loading control for the Western blot. Sanger sequencing shown below the schematic was done to verify deletions in all male Kdm6aΔE clones compared to wt. (B) Schematic shows location of the promoter region deletion (Kdm6aΔP) made in male BC and CB ES cells. Exons are shown as vertical bars with location of the PCR and RT-PCR primers (color-coded arrows) used to confirm the deletion and measure expression. Above, images of gels after electrophoresis of PCR products (gDNA) and RT-PCR products (cDNA) confirm KDM6A KO in male BC clone (Kdm6aΔP4) and CB clone (Kdm6aΔP2.1) compared to wt. Actb was run as a control for PCR and RT-PCR. Sanger sequencing shown below the schematic was done to verify deletions in all Kdm6aΔP male BC and CB clones. (C-D) UCSC genome browser (GRCm38/mm10) view of RNA-seq profiles for (C) BC male and female wt, BC male clones Kdm6aΔE and Kdm6aΔP, and (D) CB male wt, and CB male Kdm6aΔE and Kdm6aΔP clones. Note that a low level of 3’ end reads are present in clones with deletion of exons 2–4, although there is no evidence of protein (see A above). In clones with deletion of the promoter, there are no reads over Kdm6a, including no reads overlapping the small alternative Kdm6a transcript, confirming absence of expression and suggesting the absence of a cryptic promoter for the smaller annotated transcript. Actb was used as a PCR control. Primer sequences and sgRNAs are listed in Additional file 8. (E) Evidence of an off-target event in male BC clones Kdm6aΔE1 and Kdm6aΔE3. RNA-seq genome browser profiles in BC derived wt and Kdm6a KO male clones show absence of Pdha1 reads in Kdm6aΔE but not Kdm6aΔP clones. The inset shows PCR analysis for regions A, B and C amplified by primers indicated by colored arrows, which confirms an off-target deletion of part of Sh3kbp1 and all of Map3k15 and Pdha1 in BC Kdm6aΔE clones. No other KO clones show this off-target deletion. Actb was run as a control.

Additional file 2: Fig. S2.

Diploid gene expression changes in Kdm6a KO BC and CB clones. (A) Principal component analysis (PCA) based on diploid expression values for all transcribed genes from RNA-seq in wt and Kdm6a KO clones derived from the BC and CB crosses. BC clones include two male wt, two female wt, three male KO (Kdm6aΔ/Y), while CB clones include two male wt and three male KO (Kdm6aΔ/Y). Clone identifiers are included in the plot. (B) Hierarchal clustering of the clones described in (A). The color scale represents the sample-to-sample distance. (C) Expression fold changes between BC derived Kdm6aΔ/Y clones (Kdm6aΔE1, Kdm6aΔE3, and Kdm6aΔP4) and wt measured by RNA-seq show a significant decrease in Kdm6a expression, but no significant decrease in expression of pluripotent genes (Nanog, Pou5f1, Sox2, Cd9, Stat3, Fut4) (**p < 0.01 using a student’s t-test). (D) Deficiencies in ES cell differentiation potential following Kdm6a KO were tested by removal of LIF and by all-trans retinoic acid (RA) treatment. Smaller and less dense embryoid bodies were observed 6 days after removal of LIF in Kdm6aΔE1. In the presence of RA, wt cells show morphological signs of differentiation after 8 days, while Kdm6aΔE1 cells remain similar in morphology to day 0 controls. (E) Expression fold changes between male BC Kdm6aΔ/Y clones (Kdm6aΔE1, Kdm6aΔE3, and Kdm6aΔP4) and wt measured by RNA-seq confirm decreased expression of known KDM6A target genes Wt1, Bcar3, Foxh1, Dnmt3a, Sox3, Hsd17b11. Expression is based on diploid analysis at **FDR < 0.001. (F) Expression fold changes between male BC Kdm6aΔ/Y clones (Kdm6aΔE1, Kdm6aΔE3, and Kdm6aΔP4) and wt measured by quantitative RT-PCR analysis confirm downregulation of Vrtn, Bcar3, and Sall2 (*p < 0.01 using a student’s t-test) (Additional file 9A). Expression is normalized to Actb. (G) Venn diagrams to compare the number of downregulated and upregulated DEGs in male BC Kdm6aΔ/Y clones (Kdm6aΔE1, Kdm6aΔE3, and Kdm6aΔP4) and male CB Kdm6aΔ/Y clones (Kdm6aΔE2.5, Kdm6aΔE2.7, Kdm6aΔP2.1). DESeq2 was employed to identify differentially expressed genes (DEGs) in each cross using an FDR cutoff of < 0.05 and a ≥ 1.5 fold-change.

Additional file 3: Fig. S3.

Changes in sex-biased gene expression in Kdm6a KO BC and CB clones. (A) Scatter plots of log2 gene expression between male and female lines illustrating the degree of sex-biased expression in wt cells. Sex-biased genes were identified by comparing two BC male wt clones versus two BC female male wt clones. A gene was classified as sex-biased if its expression is ≥ 2 TPM fold expression higher in either female wt or male wt cells (p ≤ 0.05). Female-biased genes are in light orange and male-biased genes in light blue. (B) Scatter plot showing loss of sex-biased expression in three BC male Kdm6aΔ/Y clones (Kdm6aΔE1, Kdm6aΔE3, and Kdm6aΔP4) versus two BC female wt clones. Genes that lost sex-biased expression in KO cells are in dark orange if female-biased in wt, and in dark blue if male-biased in wt. (C) Scatter plot of genes without sex-biased expression (< 2 TPM fold expression difference) in wt BC cells. (D) Gain of sex biased expression in in three BC male Kdm6aΔ/Y clones (Kdm6aΔE1, Kdm6aΔE3, and Kdm6aΔP4). (E, F). Scatter plots of log2 gene expression of sex-biased genes in (E) three CB male wt clones versus three CB female wt clones, (F) three CB male Kdm6aΔ/Y (Kdm6aΔE2.5, Kdm6aΔE2.7, and Kdm6aΔP2.1) clones versus three CB female wt clones. Same analysis as in A, B. All wt expression values are from re-analysis of published RNA-seq data in CB ES cells, as we did not have access to a female wt CB line [34]. See also Table 1 and Additional file 10: Table S10.

Additional file 4: Fig. S4.

Correlation of gene expression (diploid, maternal and paternal alleles) between wt and KO clones. Scatter plots of average gene expression based on TPM values show a high correlation between male BC wt and (A) Kdm6aΔ/Y clones (Kdm6aΔE1, Kdm6aΔE3, and Kdm6aΔP4), and (B) male CB Kdm6aΔ/Y clones (Kdm6aΔE2.5, Kdm6aΔE2.7, and Kdm6aΔP2.1) for diploid and allele-specific expression. Correlation coefficients are all ≥ 0.96. (C, D) Scatter plots of diploid TPM expression values comparing individual BC KO clones (C) and CB KO clones (D). Red dots represent all allelicly regulated genes (combined groups A-F). Black dots represent all other genes. High correlation coefficients (c.c) indicate similar expression changes following KO in all three clones.

Additional file 5: Fig. S5.

A subset of maternally expressed imprinted genes are downregulated after Kdm6a KO in BC but not in CB clones. (A) Scatter plot of average expression changes for canonical imprinted genes, either maternally (purple) or paternally (blue) expressed, in Kdm6aΔ/Y (Kdm6aΔE1, Kdm6aΔE3, and Kdm6aΔP4) versus wt (average from 2 clones) based on diploid RNA-seq analysis in male BC derived clones. Only imprinted genes with > 1TPM in at least one sample are shown. (B) RT-PCR analysis confirms decreased expression of Meg3 and Rian in all male BC KO clones (Kdm6aΔE1, Kdm6aΔE3, Kdm6aΔP4). CRISPR + /– labels indicate presence/absence of transfection with CRISPR/Cas9 and Kdm6a sgRNAs, while a ± for Kdm6aΔP and Kdm6aΔE clones indicates editing result. Non-edited CRISPR control clones (Kdm6aE14, Kdm6aP2) maintain expression of these genes. (C) Sanger sequencing in wt male BC clones confirmed the presence of SNPs (gDNA) and mono-allelic expression (cDNA) of Meg3 and Rian. (D, E) UCSC Genome browser (GRCm38/mm10) views of allele-specific RNA-seq profiles at the maternally expressed Dlk1/Mirg polycistron region and at a subset of imprinted Xlr genes. Allelic read profiles on the maternal chromosome (purple) and the paternal chromosome (blue) are based on SNP analysis in male BC wt clone (KO-), compared to male BC Kdm6aΔE1 clone (KO +). (F) UCSC Genome browser (GRCm38/mm10) view of allele-specific RNA-seq, ATAC-seq and H3K27me3 ChIP-seq profiles generated in male wt and Kdm6aΔE1 at the imprinted genes H19 (maternally expressed) and Igf2 (paternally expressed). The DMR (differentially methylated region) is denoted by a green bar. (G) TPM expression fold change measured by RNA-seq for H19 and Igf2 in Kdm6AΔE1 versus wt shows H19 downregulation concordant to the expected upregulation of Igf2. (H) Same analysis as in (A) for male CB Kdm6a Δ/Y clones (Kdm6aΔE2.5, Kdm6aΔE2.7, and Kdm6aΔP2.1). (I, J) Venn diagrams to compare the number of allele-biased genes in BC male wt, BC female wt, and CB male wt. Allelic bias was defined as a ≥ twofold difference between alleles. BC derived cells show more commonality in allele-biased genes than with CB derived cells.

Additional file 6: Fig. S6.

Chromatin modifier expression and chromatin accessibility changes at enhancers and promoters after Kdm6a KO. (A) Expression fold changes between Kdm6aΔ/Y clones (Kdm6aΔE1, Kdm6aΔE3, and Kdm6aΔP4) and wt measured by RNA-seq for genes encoding known chromatin modifying proteins in the PRC2 complex (Ezh1/2, Suz12, Eed, Rbbp7) and the MLL complex (Ep300, Smarca4, Eomes). (B) Box plots of allelic ATAC-seq read coverage at gene enhancers in male BC and CB clones (Kdm6aΔE1 and Kdm6aΔE2.5) compared to their respective wt controls. Enhancers (n = 25,346) were defined as regions enriched in H3K4me3 and H3K27ac in male ES cells as described [12]. (C) Same analysis as in (B) but for promoters (n = 20,745). ATAC-seq coverage was calculated using ± 2 kb surrounding the transcription start sites determined using GENCODE. No significant chromatin accessibility differences were seen, nor was there a significant difference between alleles in either cross. (D) Same analysis as is (B) but for H3K27me3.

Additional file 7: Fig. S7.

Allelic epigenetic features of genes regulated by KDM6A. (A) Scatter plots of H3K27me3 read coverage in log2 scale on maternal and paternal alleles in male BC Kdm6aΔE1 cells versus wt. Allele-specific ChIP-seq reads around promoters (± 2 kb of the TSS) were used to calculate promoter coverage normalized by sequencing depth and allele differences (see also Fig. 5A, D). (B) Same analysis as in (A) but for DEGs that overlap between exonic KO clones only (Kdm6aΔE1 and Kdm6aΔE3) and Kdm6aΔY clones (Kdm6aΔE1, Kdm6aΔE3 and Kdm6aΔP4). (C) Scatter plots of ATAC-seq read coverage in log2 scale on maternal and paternal alleles in male BC Kdm6aΔE1 cells versus wt. Allele-specific reads around promoters (± 2 kb of the TSS) were used to calculate promoter coverage normalized by sequencing depth and allele differences (see also Fig. 5B, E). (D) Same analysis as in (C) but for DEGs that overlap between exonic KO clones only (Kdm6aΔE1 and Kdm6aΔE3) and Kdm6aΔY clones (Kdm6aΔE1, Kdm6aΔE3 and Kdm6aΔP4). (E, F) Same analysis as in (A) and (B), but for upregulated DEGs. (G, H) same analysis as in (C) and (D), but for upregulated DEGs. (I, J) Scatter plots of ratios of ATAC-seq read coverage versus ratios of H3K27me3 ChIP-seq read coverage between CB clone Kdm6aΔE2.5 and wt CB cells at promoter regions (± 2 kb of the TSS) of downregulated groups on maternal and paternal alleles. (K, L) Same analysis as in (I) and (J) but for upregulated groups. See also Fig. 5G, H. Unedited image figure legend. (A) Raw images for electrophoresis gels shown in Additional file 1A for BC and CB clones. Left panel shows screening of CRISPR exon deletion clones by PCR. The top right gel shows screening of CRISPR exon deletion clones in CB cells by PCR (gDNA). The associated control Actb lanes are in (C) lower right gel (gDNA). For all panels, primer names correspond to those in Additional file 1A. (B) Unedited images of western blots and controls for KDM6A. Top two are film images after 30 s and 5 min exposure (30 s image is shown in Additional file 1A). Lower image shows a picture of the Ponceau S stain used for loading control. For all images the clones are labelled at the top. (C) Raw gels for Kdm6a expression changes in exon deletion clones. Left panel shows Kdm6a expression changes in BC clones. The center panel shows a repeat electrophoresis run of the Actb controls which is shown in Additional file 1A. For CB clones, Kdm6a expression in controls and KO cells are shown in the lower right gel in (A) (cDNA). The right panel shows the Actb control (cDNA). Primer names correspond to those in Additional file 1A. (D) Raw images for electrophoresis gels shown in Additional file 1B for BC and CB clones. Left panel shows promoter deletion screening by PCR in BC clones (gDNA). The center and right panels show the same for the CB promoter clone (gDNA). For all panels, primer names correspond to those in Additional file 1B. (E) Raw gels for Kdm6a expression changes in control and promoter deletion clones. Top gel shows expression changes in BC KO clones and the lower gel shows expression changes in the CB KO clone. Primer names correspond to those in Additional file 1B. (F) Gel image of PCR screening for the off-target effect seen in CRISPR KO BC exon clones ΔE1, and ΔE3. Primer names correspond to those in Additional file 1F. (G) Gel images showing expression changes of Meg3 (left) and Rian (center) in controls and BC KO clones. The right panel shows the Actb controls.

Additional file 8: Table S8.

A: Summary of clones, B: Oligos used for CRISPR and PCR.

Additional file 9: Table S9.

A: Expression changes in BC Kdm6aΔY clones (see Figure 1B); B: GSE97702 comparison; C: Expression changes in CB Kdm6aΔY KO clones (see Figure 1C); D: Downregulated genes common between BC and CB KO clones (see Figure 1D-F); E: Upregulated genes common between BC and CB KO clones (see Figure 1D-F).

Additional file 10: Table S10.

A: Sex-biased genes in wt ES cells (see Figures 2A and Additional file 3A); B: Loss of sex-biased expression in Kdm6aΔY ES cells (see Figures 2 and Additional file 3B); C: Concordant sex-biased genes (Werner et al. 2017); D: Loss of sex-biased expression in Kdm6aΔY for concordant genes; E: Loss of sex-biased expression in CB Kdm6aΔY clones (see Additional file 3F).

Additional file 11: Table S11.

A: BC group A (see Figures 3A, B); B: BC group B (see Figures 3A, B); C: BC group D (see Figures 3A, B); D: BC group E (see Figures 3A, B); E: CB group A (see Figures 3C, D); F: CB group B (see Figures 3C, D); G: CB group D (see Figures 3C, D); H: CB group E (see Figures 3C, D).

Additional file 12: Fig. S12.

A: Expression of BC Group A genes in CB clones; B: Expression of BC Group B genes in CB clones.

Additional file 13: Fig. S13.

A: BC H3K27me3 ChIP-seq changes (see Figures 5 and Additional file 7); B: BC ATAC-seq changes (see Figures 5 and Additional file 7); C: CB H3K27me3 changes (see Figures 5 and Additional file 7); D: CB ATAC-seq changes (see Figures 5 and Additional file 7).

Additional file 14: Fig. S14.

Summary of RNA-seq mapping.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Ma, W., Fang, H., Pease, N. et al. Sex-biased and parental allele-specific gene regulation by KDM6A. Biol Sex Differ 13, 40 (2022). https://doi.org/10.1186/s13293-022-00452-0

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13293-022-00452-0

Keywords

  • Sex biases
  • Parent-of-origin
  • Allele-specific
  • Imprinting
  • Gene regulation
  • Development
  • Epigenetics
  • Histone methylation
  • Chromatin