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

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 KOFig. 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

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

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].