Prevalent bursty expression contributes to dynamic aRME during iPSC reprogramming

Transcription of many genes occurs in a stochastic manner, where genes undergo sporadic bursting to produce RNA. The kinetics of transcriptional bursting are deduced through the well-known “two-state model” of transcription. According to the “two-state model,” the promoter of a gene switches stochastically from an inactive/OFF state to an active/ON state, and burst kinetics is determined through two parameters: burst frequency and burst size (Fig 1A) (Chubb et al, 2006; Raj & van Oudenaarden, 2008; Larson, 2011). Burst frequency is defined by the number of bursts per unit time, and burst size explains the average amount of mRNA molecules per burst (Fig 1A). To date, burst kinetics have been studied mainly at the gene level, not the allelic level. However, profiling burst kinetics at allelic resolution is important as the kinetics of bursting between alleles often differ (Naik et al, 2021). Therefore, to understand how individual alleles of a gene contribute to the overall gene expression, we have delineated the kinetics of allelic bursting genome-wide across different stages of reprogramming of MEF to iPSC (Fig 1B). To profile, genome-wide allelic transcriptional burst kinetics, we performed SCALE (Single-Cell ALlelic Expression) analysis using scRNA-seq datasets (Fig 1B). SCALE relies on the Empirical Bayes Framework, which first classifies the genes into monoallelic, biallelic, and silent based on the allele-specific read counts across the cells. Subsequently, the biallelic genes are further categorized into biallelic bursty and biallelic constitutive (Fig S1A) (Jiang et al, 2017). Biallelic genes that have zero expression of both alleles in a few cells but not all cells are considered bursty genes, whereas genes that have at least one allele constitutive expression or both allele constitutive expression across cells are considered biallelic constitutive genes (Fig S1A). The MEF cells used in the experiment were derived from a cross of two divergent mouse strains M. Musculus (129S1) and M. Castaneous (CAST), thereby enabling us to perform allele-specific analysis based on strain-specific SNPs (Fig 1B) (Talon et al, 2021). We excluded low-expressed genes from our study to avoid allelic dropout-related technical noise, which could lead to inaccurate estimation of the allelic expression (Kim et al, 2015; Santoni et al, 2017; Zhao et al, 2017).

First, we performed UMAP-based clustering to identify homogeneous cell populations across the different stages of reprogramming. We found that day 0-MEF cells and iPSC cells formed individual tight clusters (Fig 1C). Intermediate cells (day 8, day 9, day 10) also clustered together. However, day 12 cells formed two different sub-clusters (day 12A and day 12B). Therefore, we performed SCALE analysis separately for day 12A and day 12B sub-clusters to reduce the effect of heterogeneity in our burst kinetics analysis (Fig 1C). Through SCALE analysis, we found that most of the biallelic genes (∼90%) have bursty expression across all stages of reprogramming (Fig 1D). Next, we compared the bursty/constitutive pattern of gene expression across different reprogramming stages and found that most genes (n = 756) maintained bursty expression throughout the reprogramming (Fig 1E). Interestingly, these genes were enriched toward key reprogramming-related processes such as ribosome biogenesis, proteasome assembly, blastocyst growth, etc. (Fig 1F). Moreover, genes (n = 60) that remained bursty on all days except iPSCs were also enriched in many reprogramming-related biological processes, including chordate embryonic development (Fig 1F). These results indicated that many genes involved in iPSC reprogramming pathways exhibit bursty expression. However, a few genes (n = 20) that maintained bursty expression in all-day points except MEF did not show enrichment towards major iPSC reprogramming pathways (Fig S1B). Genes (n = 18), which remained constitutive throughout the reprogramming, were not enriched to major reprogramming pathways (Fig S1B). On the other hand, we have previously shown that bursty expression results in dynamic aRME, which creates gene expression heterogeneity among cells (Naik et al, 2021). Therefore, we explored the landscape of allelic expression throughout the different stages of iPSC reprogramming by allele-specific scRNA-seq analysis. We categorized a gene as monoallelic within a cell if at least 95% of the expression originated from one allele. We classified genes based on their allelic expression patterns into four categories: Cat 1: nonrandom monoallelic (same allele expressed across the cells); Cat 2: RME with one allele (genes expressed biallelically in some cells whereas exhibiting monoallelic expression from one allele consistently in other cells); Cat 3: random monoallelic with either allele (genes expressed biallelically in some cells whereas exhibiting monoallelic expression from either of the allele in other cells); and Cat 4: biallelic (genes expressed biallelically across the cells) (Fig 1G). Our analysis revealed that ∼80–90% of genes have dynamic aRME across different stages of reprogramming (Fig 1G). Analysis of allelic expression of X-linked genes in MEF showed monoallelic expression from the CAST allele as 129S1-X is inactivated in these cells, thereby validating our allele-specific analysis pipeline (Fig S1C). X-inactivation is a process through which female mammals compensate for the dosage of X-linked gene expression between sexes (Gayen et al, 2015, 2016; Sarkar et al, 2015; Saiba et al, 2018; Kaur et al, 2020; Samanta et al, 2022). Altogether, our analysis suggested widespread bursty expression of genes during iPSC reprogramming, resulting in dynamic aRME.

Alleles of genes exhibit similar burst kinetics but have different coordination during iPSC reprogramming

Next, we investigated if the two alleles of a gene have similar burst kinetics or not. To explore this, we profiled burst frequency and burst size at the allelic level for the biallelic bursty genes using SCALE. SCALE relies on a “two-state model” of transcription where a gene switches from active to inactive state with activation rate of Kon and deactivation rate of Koff. When a gene is in the active state, the rate of transcription is S and the rate of RNA decay is d (Fig 1A). Burst frequency is determined by the number of bursts per unit time (Kon), and burst size is deduced through average number of mRNA molecules per burst when the gene is in active state (S/Koff) (Fig 1A). We observed a high degree of correlation of both burst frequency (r = 0.64–0.778) and burst size (r = 0.669–0.803) between two alleles across all-day points (Fig 2A and B). Very few genes exhibited significant burst frequency and size differences between the two alleles, as marked by the red triangles (Fig 2A and B). Taken together, our results suggest that alleles of most of the genes have similar burst kinetics. Next, we explored the degree of coordination of bursting between two alleles by plotting the percent of cells expressing neither allele (p₀) versus the percent of cells expressing both alleles (p₂) (Fig 2C). Blue diagonal line represents perfect coordination (p0 + p2 = 1), whereas the red curve signifies independent bursting with shared kinetics (Fig 2C). We categorized genes based on the degree of coordination of allelic bursting into three categories: (1) highly coordinated genes: with p0 + p2 > 0.90 marked by gray asterisk between blue dotted diagonal lines; (2) independent genes: genes near the red curve with a threshold of +0.05, signified by the upper red curve, and −0.05, signified by the lower red curve, marked by rosewood triangles; (3) semi-coordinated genes that lie between the uppermost red curve line and lower blue dotted diagonal line (Persian blue dots) (Fig 2C).

We observed that most of the genes exhibited semi-coordinated allelic bursting at all-day points. However, many genes also showed highly coordinated allelic bursting across the different stages of reprogramming (Fig 2C).

Genes involved in iPSC reprogramming undergo highly coordinated allelic bursting

Next, we performed a gene ontology (GO) enrichment analysis of highly coordinated genes across all stages of reprogramming (Fig 3). Interestingly, we found that highly coordinated genes in iPSCs are enriched in processes linked to stem cell population maintenance, in-utero embryo development, endoderm formation, etc. (Fig 3A, Supplemental Data 1). Furthermore, highly coordinated genes in iPSCs are enriched in cellular respiration, cristae formation, and glutamine metabolism, which are relevant to metabolic remodeling in iPSCs (Fig 3A) (Teslaa & Teitell, 2015; Mathieu & Ruohola-Baker, 2017; Tohyama et al, 2017). Importantly, genes involved in ribosome biogenesis, crucial for stem cell maintenance, were highly enriched into the highly coordinated gene cohort in iPSC (Fig 3A) (Gabut et al, 2020). On the contrary, highly coordinated genes in day 0 MEF cells were not enriched in such related pathways. Interestingly, across the different intermediate stages of reprogramming, highly coordinated genes were enriched towards many reprogramming-related biological processes such as ribosome biogenesis, protein folding, stem cell maintenance, NFkB signaling, etc. (Fig 3B) (Gabut et al, 2020; Yan et al, 2020; Kaltschmidt et al, 2021). Next, we performed a cross-comparison of the three allelic coordination categories across all-day points and found that, whereas many genes maintained a similar degree of allelic coordination across all-day points, some genes did not (Fig S2A). Strikingly, we found that the genes (n = 52) that become highly coordinated in iPSCs are enriched towards ribosome biogenesis, aerobic respiration, nuclear pore formation, which are highly critical and directly relevant for iPSC reprogramming (Fig 3C) (Xu et al, 2013; Teslaa & Teitell, 2015). Moreover, genes that remain highly coordinated (n = 25) on all days showed enrichment towards some iPSC reprogramming-linked functions like protein stability, NFkB pathway, mitochondrion organization, etc. (Fig 3C). Altogether, our analysis revealed that the two alleles of many genes involved in iPSC reprogramming have a high degree of transcriptional bursting coordination (Supplemental Data 1). In parallel, we also observed that genes that remained semi-coordinated (n = 266) on all days or converted from highly coordinated on day 0 to semi-coordinated on other days (n = 33) were enriched towards some iPSC reprogramming-related pathways like proteasome assembly, oxidative stress response, etc. (Fig S2B) (Szutorisz et al, 2006; Schröter & Adjaye, 2014; Hawkins et al, 2016; Mathieu & Ruohola-Baker, 2017). Next, we investigated if higher coordination of transcriptional bursting between alleles of individual genes facilitates the optimal expression dosage of these highly coordinated genes involved in iPSC reprogramming. To explore this, we compared the overall expression level of highly coordinated, semi-coordinated, and independent genes. We found that the indeed highly coordinated genes have higher expression level compared to the semi-coordinated and independent genes (Fig 3D). Taken together, we conclude that highly coordinated allelic bursting helps to fine-tune the optimal expression of genes involved in iPSC reprogramming.

Coordinated allelic bursting is linked to chromatin accessibility

To understand the mechanisms of allelic bursting coordination, we asked whether the degree of coordination of allelic bursting is linked to epigenomic states. To address this, we profiled genome-wide allelic chromatin accessibility across different stages of MEF to iPSC reprogramming through allele-specific analysis of available ATAC-sequencing (ATAC-seq) datasets (Talon et al, 2021). The same hybrid MEFs (129S1X CAST), as described for scRNA-seq, were used for this experiment, allowing us to profile chromatin accessibility at the allelic level. We analyzed ATAC-seq in MEFs (day 0) and across reprogramming stages (SSEA1+ reprogramming intermediates at days 8, 9, 10, 12, and iPSCs), like the burst kinetics analysis (Fig 4A). We first validated our allele-specific ATAC-seq analysis pipeline by quantifying the difference in the enrichment of ATAC-seq reads between active-X (129S1) versus inactive-X (CAST allele) (Fig S3A). In consistence with previous reports, in MEFs and early reprogramming intermediates, active-X (CAST) showed strong enrichment of ATAC-seq reads, whereas the inactive-X (129S1) showed almost no enrichment, and upon reactivation of the inactive-X towards the attainment of iPSCs, there was a gain of chromatin accessibility (Fig S3A). Taken together, enrichment analysis of ATAC-seq reads of X-linked genes validated the accuracy of our method. Next, we compared the enrichment of ATAC-seq reads between two alleles of different categories of genes (highly coordinated, semi-coordinated, and independent) across the gene body and 3 kb upstream of TSS and 3 kb downstream of TES during reprogramming (Fig 4A). Interestingly, our analysis revealed that the two alleles of highly coordinated genes have very similar enrichment at most day points. Whereas enrichment of ATAC-seq reads of alleles of semi-coordinated/independent genes differed in most cases (Fig 4A). As expected, allelic enrichment of ATAC-seq reads considering all autosomal genes was quite similar (Fig 4A).

Notably, genes that maintained semi-coordinated bursting throughout reprogramming maintained differences in allelic accessibility (Fig 4B). Similarly, genes that maintained highly coordinated bursting throughout reprogramming always maintained similar allelic accessibility (Fig 4B). Furthermore, we show that allelic accessibility differences reduce upon becoming highly coordinated in iPSC after maintaining a semi-coordinated bursting throughout the other day points (Fig S3B). A similar trend was found for genes which switched from highly coordinated in MEF to semi-coordinated on other days (data not shown). Together, our analysis suggested a positive correlation between the coordination of allelic bursting and the similarity of allelic chromatin accessibility.

Next, we explored if allelic accessibility differences in semi-coordinated and independent genes are associated with the differential binding of transcription factors (TFs) between alleles of individual genes. To test this, we determined TFs binding scores of individual alleles of a gene using TOBIAS (Transcription factor Occupancy prediction By Investigation of ATAC-seq Signal). Interestingly, we found that many TFs had a significantly different binding score between individual alleles of semi-coordinated and independent genes, which was not the case for highly coordinated genes (Fig 5). Taken together, our analysis suggests that allelic accessibility differences in semi-coordinated or independent genes allow differential binding of certain TFs among alleles of genes, which in turn leads to semi-coordinated or independent transcriptional bursting.

Chromatin accessibility factors contribute to allelic bursting coordination

Next, we investigated if the degree of allelic coordination is dependent on chromatin accessibility factor enrichment. To explore this, we profiled the enrichment of different chromatin accessibility factors in the gene body and across 2 kb upstream of TSS and 2 kb downstream of TES of highly coordinated, semi-coordinated, and independent genes in MEF and iPSCs (Figs 6 and S4). Interestingly, we found that many accessibility factors such as H3.3, H3K36me3, H3K27ac, H3K4me3, H3K4me2, H3K4me1, and RNAPolII-S2P are highly enriched on highly coordinated genes compared to the semi-coordinated/independent genes in both MEF and iPSC (Figs 6 and S4). Moreover, in iPSC, we found that highly coordinated genes are enriched with chromatin remodeler BRD4 (Fig 6). H3K79me2, H3K9ac, and RNA PolII-S5P were highly enriched to the highly coordinated genes in MEF cells (Figs 6 and S4). Next, we tested if switching of coordination pattern of genes in MEF to iPSC is associated with changes in the pattern of chromatin accessibility factor enrichment. We found that highly coordinated genes in MEF having higher enrichment of H3K36me3, H3K27ac, H3.3, and RNA PolII-S2P compared to the semicordinated genes, which flip their pattern of enrichment in iPSC upon conversion of their coordination pattern (Fig S5A and B). Whereas genes that do not switch their coordination pattern in MEF to iPSC maintain the enrichment pattern of these marks (Fig S5A and B). Altogether, our analysis suggested that higher enrichment of these chromatin accessibility factors ensures highly coordinated allelic bursting. However, other chromatin modifications, such as H3K9me3 and CTCF, did not show such differences in enrichment between different categories of genes (Fig S4).

Next, we investigated if the chromatin states of highly coordinated, semi-coordinated, and independent genes are different. To explore this, we identified 10 chromatin states in both MEFs and iPSCs using the ChromHMM based on the seven histone marks (H3K27ac, H3K27me3, H3K36me3, H3K4me1, H3K4me2, H3K4me3, and H3K9me3) using the available ChIP-seq datasets (Fig S6). We annotated each state putatively based on related functional element enrichment. Next, we compared the enrichment of these different chromatin states in the gene body in the TSS/TES neighbourhood of highly coordinated, semi-coordinated, and independent genes. Interestingly, we found that chromatin states enriched around strong enhancers showed higher enrichment in highly coordinated genes than the semi-coordinated or independent genes in both MEF and iPSC (Fig S6A and B). Notably, independent genes were highly enriched with bivalent chromatin state compared to the highly coordinated genes in both MEF and iPSC (Fig S6A and B). On the other hand, enrichment of other chromatin states among highly coordinated, semi-coordinated, and independent genes did not exhibit much difference. Taken together, our analysis suggests that enrichment of enhancer elements in highly coordinated genes might play a role in coordinating the allelic transcriptional bursting.

Allele-specific expression and burst kinetics analysis

To obtain allelic read counts from scRNA-seq data, we performed allele-specific expression analysis following the pipeline as reported previously (Mandal et al, 2020; Naik et al, 2021, 2022). Briefly, we first constructed in silico CAST/EiJ and 129S1-specific parental genome by incorporating CAST/EiJ or 129S1-specific SNPs into the mm10 genome using variant calling file tool (Danecek et al, 2011). Variant calling file was downloaded from the mouse genome project (https://www.sanger.ac.uk/science/data/mouse-genomes-project). Next, we aligned RNA-seq reads into both parental genomes using STAR aligner (STAR-2.7.10a), allowing no multi-mapped reads (--outFilterMultimapNmax 1). We filtered out those genes for allele-specific read counts which had at least two informative SNPs and a minimum of three reads per SNP site. We took an average across SNPs to get gene-level allelic read counts. We normalized the allelic read counts by reads per kilobase million. Because low-expressed genes are dropout prone, we removed low-expressed genes from our analysis to avoid potential dropout effect (Kim et al, 2015; Santoni et al, 2017; Zhao et al, 2017). We considered only those genes, which were expressed in at least ∼27% of cells and had mean expression >10 reads per kilobase million. Allelic ratio was calculated individually for each gene using the formula = (129S1/CAST reads) ÷ (129S1 + CAST reads). A gene was considered monoallelic if at least 95% of the allelic reads came from only one allele. We performed genome-wide allele-specific burst kinetics analysis using SCALE (Jiang et al, 2017). In brief, SCALE relies on Empirical Bayes Framework, which first classifies the genes into monoallelic, biallelic, and silent based on the allele-specific read counts and deduces the allelic burst kinetics based on the two-state model of transcription. It infers different burst kinetic parameters at the allelic level, such as burst frequency (K_on) and burst size (S/K_off), as described in the results. We excluded X-linked genes for our burst kinetics analysis. We performed SCALE analysis on each time point of the reprogramming separately.

Allelic ATAC-seq analysis

For allele-specific ATAC-seq analysis, first we created an “N-masked reference genome mm10” through substituting strain-specific (129S1/SvImJ and CAST/EiJ) SNP position with “N” using SNPsplit_genome_preparation (0.5.0) (Krueger & Andrews, 2016). Strain-specific SNPs were obtained as described above. Next, reads for all samples were mapped to this N-masked genome using Bowtie2 (Langmead & Salzberg, 2012). We removed duplicate reads and mitochondrial reads from our analysis. SNPsplit was then used to create allele-specific BAM files by segregating the aligned reads into two distinct alleles (129S1/SvImJ and CAST/EiJ). Bigwig files were generated from these allelic BAM files using deepTools (version 3.5.1) function bamCoverage (--binSize 100 --smoothLength 500 --normalizeUsing RPGC) (Ramírez et al, 2016). Allelic enrichment metaplots were generated for different gene categories using deepTools functions computeMatrix and plotProfile. Few genes were excluded from our analysis because of the lack of SNP for each day: day 0–15 (of 1,745); day 8–13 (of 1,849); day 9–14 (of 1,906); day 10–19 (of 1,982); day 12A–21 (of 2,521); day 12B–23 (of 2,562); iPSC—14 (of 1,863).

ChIP-seq analysis

For analysis of chromatin mark enrichment in MEF and iPSC in different categories of genes, we analyzed available ChIP-seq data. For ChIP-seq analysis, the reads were mapped to either mm9 or mm10 genome using Bowtie tool (Langmead & Salzberg, 2012). We removed duplicate reads, mitochondrial reads from our analysis. Different enrichment metaplots were generated for different gene categories using deepTools (version 3.5.1) functions computeMatrix and plotProfile.

Transcription factor binding analysis using ATAC-seq footprinting

Transcription factor occupancy between individual alleles of genes was predicted using Transcription factor Occupancy prediction By Investigation of ATAC-seq Signal (TOBIAS, version 0.16.0) (Bentsen et al, 2020). Allelic ATAC-seq BAM files were first corrected for insertion bias of the Tn5 transposase using the command ATACorrect for each gene category (highly coordinated, semi-coordinated, and independent) separately. Gene-wise allelic peak regions were identified through allelic peak calling using macs2 (v 2.2.7.1), and peaks were annotated using uropa (v 4.0.3). Next, footprinting scores for individual alleles were generated using the command ScoreBigwig using the ATACorrect data.

Allelic footprint scores were matched to the curated list of JASPAR motifs using the command BINDetect for both alleles separately. Then the resulting scores of each TF motifs for both alleles were correlated using Pearson correlation method.

Identifications of chromatin states

Chromatin states in both MEFs and iPSCs were profiled using the ChromHMM (v 1.24) based on the seven chromatin marks (H3K27ac, H3K27me3, H3K36me3, H3K4me1, H3K4me2, H3K4me3, and H3K9me3) using the available ChIP-seq datasets (Ernst & Kellis, 2012). In brief, we first constructed binarized bam files of the chromatin marks required as input for ChromHMM using the binarizeBam function of ChromHMM v 1.24. Chromatin state models were generated with the LearnModel function with 10 states. Then we manually curated to assign putative functions to each state according to the combination of epigenomic marks. Finally, enrichment for each category of genes was shown across the states using the OverlapEnrichment function of ChromHMM.

Quantification and statistical analysis

All statistical analysis was performed using the R software (https://www.R-project.org/). Mann–Whitney two-sided U test was used for statistical significance analysis, and P-values < 0.05 was considered significant. For correlation analysis, the Pearson test was used.