Ethical compliance

The research carried out in this study has been approved by the Swedish Board of Agriculture, Jordbruksverket: N343/12.

Cell culture

Mouse primary fibroblasts were derived from adult (>10 weeks old) CAST/EiJ × C57BL/6J or C57BL/6J × CAST/EiJ mice by skinning, mincing and culturing tail explants (for at least 10 d) in DMEM high glucose, 10% embryonic stem cell FBS, 1% penicillin/streptomycin, 1% nonessential amino acids (NEAAs), 1% sodium pyruvate, 0.1 mM 2-mercaptoethanol (Sigma-Aldrich) in culture dishes coated with 0.2% gelatin (Sigma-Aldrich). NIH/3T3 cells were maintained in DMEM supplemented with 10% FBS and 1% penicillin/streptomycin. All supplements were purchased from Thermo Fisher Scientific (unless stated otherwise).

Generation of Smart-seq2 libraries

Smart-seq2 libraries were prepared as described earlier18 using the following parameters: (1) 20 cycles of PCR for preamplification; (2) a ratio of 0.8:1 for bead:sample purification of preamplified complementary DNA (in-house-produced 22% polyethylene glycol (PEG) beads); (3) tagmentation of approximately 1 ng bead-purified cDNA (in-house-generated Tn5 (ref. 46)); (4) 10 cycles of PCR for library amplification of the tagmented samples using Nextera XT Index primers; and (5) a ratio of 1:1 for bead purification of DNA sequencing libraries (in-house-produced 22% PEG beads). Sequencing was carried out on an Illumina HiSeq 2000 generating 43 base pair (bp) single-end reads. The libraries related to Figs. 1, 3 and 6 were derived from one tail explant (F1 offspring of C57 × CAST mouse female adult) and combined with previously published Smart-seq2 data20. The additional Smart-seq2 data generated for Fig. 7 were derived from one additional tail explant (F1 offspring of CAST × C57 mouse female adult).

Generation of Smart-seq3 libraries

Smart-seq3 libraries were generated according to a previously published protocol6. Cells were stained with propidium iodide (PI) before being sorted (BD FACSMelody 100 μM nozzle; BD Biosciences) into 384 well plates containing 3 μl of Smart-seq3 lysis buffer (5% PEG (Sigma-Alrich), 0.10% Triton X-100 (Sigma-Aldrich), 0.5 U μl−1 of recombinant RNase inhibitor (Takara), 0.5 μM Smart-seq3 oligo(dT) primer (5′-biotin-ACGAGCATCAGCAGCATACGA T30VN-3′; Integrated DNA Technologies), 0.5 mM deoxynucleoside triphosphate (Thermo Fisher Scientific)) and stored at −80 °C. From this point, the standard protocol for Smart-seq3 was applied: (1) 20 cycles of PCR for preamplification of cDNA; (2) a ratio of 0.6:1 for bead:sample purification of preamplified cDNA (in-house-produced 22% PEG beads); (3) tagmentation of 150 ng bead-purified cDNA using 0.1 μl of Amplicon Tagment Mix; and (4) 12 cycles of PCR for library amplification of the tagmented samples using custom-designed Nextera Index primers containing 10-bp indexes. Samples were pooled, bead-purified at a ratio of 0.7:1 (in-house-produced 22% PEG beads) and prepared for sequencing on a DNBSEQ-G400RS (MGI) generating 100-bp paired-end reads. The data related to Fig. 2 were obtained from one tail explant (F1 offspring of C57 × CAST female adult mouse) and is also part of a previous study47. The libraries with siRNA-perturbed lncRNAs (related to Fig. 8) were derived from one tail explant (F1 offspring of C57 × CAST female adult mouse).

Processing of RNA-seq data

A subset of primary fibroblasts analyzed in this study (sequenced by Smart-seq2) are part of previously published studies and were reanalyzed for consistency7,20 (NCBI Sequence Read Archive ID SRP066963). The zUMIs v.2.7.1b pipeline48 was used for alignment (mm10 assembly), gene quantification (Ensembl, GRCm38.91) and allelic calling for primary fibroblast data. To pass quality control, cells were required to have (1) ≥500,000 reads, (2) 4,000 genes expressed at ≥ 5 read counts, (3) distribution of allelic counts within −0.10 < allelic SNPs < 0.10 on autosomes (imprinted and genes on the X chromosome excluded) and (4) no more than 20% of allelic counts mapped to the imprinted X chromosome (escapee genes excluded). Genes with at least five read counts in two cells were kept for downstream analysis (unless stated otherwise).

Smart-seq3 libraries of HEK293 cells had previously been generated by Hagemann-Jensen et al.6 (ArrayExpress ID E-MTAB-8735). The zUMIs v.2.7.0a pipeline48 was used for alignment (hg38 assembly) and quantification of gene expression (Ensembl, GRCh38.95). Cells were required to have (1) ≥500,000 read counts mapped to exons and (2) ≥7,500 genes (≥1 read count). Genes with at least one read count in three cells were considered for downstream analysis. Gene types were annotated according to BioMart release 91.

The Smart-seq2 libraries of mouse embryonic stem cells had previously been generated by Ziegenhain et al.19 (Gene Expression Omnibus (GEO) ID GSE75790). The zUMIs v.2.7.2a pipeline was used for alignment (mm10 assembly) and quantification of gene expression (Ensembl, GRCm38.91). Cells were required to have ≥ 400,000 read counts mapped to exons and ≥8,000 genes (≥5 read counts). Genes with at least five read counts in two cells were considered for downstream analysis and gene types were annotated according to Supplementary Table 2 (downloaded from https://m.ensembl.org/biomart/martview/; gene list also available at https://github.com/sandberg-lab/lncRNAs_bursting).

For the Smart-seq3 libraries of primary fibroblasts treated with siRNAs, the zUMIs v.2.9.4b pipeline48 was used for alignment (mm10 assembly) and quantification of gene expression (Ensembl, GRCh38.95). Cells were required to have (1) ≥100,000 read counts mapped to exons, (2) ≥50,000 unique UMI counts and (3) ≥5,000 genes (≥1 UMI count). Genes with at least one UMI count in three cells were considered for downstream analysis.

Annotation of lncRNAs

The Ensembl BioMart annotation (GRCm38.p6; Supplementary Table 2) was used to assign lncRNAs. Genes were first filtered (above) and lncRNAs categorized as: (1) divergent (no gene–gene overlap and TSS not separated by more than 500 bp); (2) convergent (gene–gene overlap and TSS not separated by more than 2 kb); (3) intergenic (no gene–gene overlap and at least 4 kb from any other expressed gene); and (4) separated transcriptional units (TSS separated with at least 4,000 bp from any other expressed gene). The threshold of 4 kb was established by manual inspection of Extended Data Fig. 1d where mean expression had been measured (median of sliding window size = 51) against the distance between the 2 most closely located TSSs (only genes passing quality control were considered for these analysis).

Permutation test for CV2

For the analysis of cell-to-cell variability, only genes meeting the following criteria were considered: (1) not imprinted; (2) not encoded on the X chromosome; and (3) being classified as separated transcriptional units (Extended Data Fig. 1d).


For each lncRNA meeting the criteria, ten separated transcribed protein-coding genes having the most similar mean expression (min(mean(RPKMlncRNA) − mean(RPKMmRNA))) were selected. The matching allowed for the same protein-coding gene to be selected multiple times (sample replacement). For the permutation test (n = 10,000), 1 expression-matched protein-coding gene was randomly sampled for each lncRNA and the expected CV2 (median) was calculated for each permutation. The P value represents the frequency of median(CV2sampled) > median(CV2lncRNA).

To estimate the number of lncRNAs needed to detect median(CV2lncRNA) > median(CV2mRNA) (Extended Data Fig. 2c), the permutation test was repeat 100 times for each subsampling size (between 10 and 200 lncRNAs) of the frequency where 50% and 95% of the permutations reached median(CV2lncRNA) > median(CV2sampled) was assessed.

Transcriptional bursting kinetics inference

Transcriptional bursting kinetics were inferred from homogenous sets of cells using the two-state model of transcription, based on previous methodology7. In detail, we first computed the UMI expression values from the Smart-seq3 libraries6 and the fraction of allele-sensitive reads were used to assign the UMI counts to the CAST or C57 allele, respectively. Cells having UMIs but lacking allelic read counts for individual genes were assigned as missing values for the inference whereas cells lacking UMIs and allelic information were considered as ‘true’ zeros and included in the analysis. The allelic expression level per cell was provided as input to the maximum likelihood inference (https://github.com/sandberg-lab/txburst); instead of using profile likelihood to estimate CIs, we performed 1,000 bootstraps per gene and allele and collected the inferred burst frequency and size of each sampled input, and importantly, each new bootstrap used a random initialization of kinetic parameter to ensure proper sampling of kinetic space. We continued with 95% CIs based on the bootstrapped parameters. For the downstream analyses we required that each gene had: (1) ≥1 UMI count in ≥5 cells; (2) burst size within 0.2 < size < 50; (3) burst frequency 0.01 < frequency < 30; (4) UMI expression 0.01 < UMImean < 100; and (5) width of CIs (CIHigh/CILow) below 101.5 (for burst size and frequency). Finally, only non-imprinted autosomal genes, identified as independent transcriptional units, were considered for downstream analysis.

Permutation test of bursting kinetics for lncRNAs with highly ranked CV2

The CV2 for each lncRNA was ranked to 100 mRNAs of similar mean expression (using allele-distributed UMIs, equally distributed with 50 mRNAs with higher or lower expression). The top 50 ranked lncRNAs, for each individual allele, were used for downstream analysis of bursting kinetics where each lncRNA was matched with 10 mRNAs of similar expression followed by subsampling 1 expression-matched mRNA for each lncRNA (similar as for Fig. 1f). The P values represent the frequency where lncRNAs (median) was higher (for burst frequencies) or lower (for burst sizes) than the burst parameters for sampled mRNAs (median).

Identification of cell cycle stage of individual primary fibroblasts

The most variable genes were identified using the R package Seurat27 v.4.0.5. Genes were first filtered for being expressed in ≥5 cells (≥5 read counts). Counts were normalized using LogNormalize (setting scale.factor = 10,000) and the most variable genes were identified using the vst method of FindVariableFeatures. We next extracted the cell cycle-related genes reported by Whitfield et al.28 (Supplementary Table 3) and used the top 50 ranked genes with the highest variability for PCA. The cell cycle phase of individual cells was identified using the first three principal components as input for the R package princurve v.2.1.6 and the Lambda factor used to align cells to the cell cycle. Expression of individual genes was illustrated using a rolling mean of 15 cells (using the R package zoo v.1.8.9). The assignment of cells to cell cycle phase was performed based on the expression levels of known cell cycle regulators (Gas1, Ccne2, Ccnb1 and Ccnd1) using the rolling mean of Seurat-normalized read counts.

Differential expression of lncRNAs in the cell cycle

Differential expression analysis between cell cycle phases (G0, G1, G1/S and G2/M) was performed using a one-way ANOVA (Benjamini–Hochberg-adjusted, P < 0.01) with normalized read counts (log-normalized, Seurat).

Correlation of cell cycle genes

Genes were first filtered for being expressed in ≥2 cells (≥5 read counts). Seurat was used to log-normalize the read counts and the normalized counts were used to calculate the Spearman correlation of cell cycle genes28. For each pairwise comparison, cells lacking expression of both genes were excluded from the analysis.

Cell cycle analysis

NIH/3T3 cells were washed twice in PBS and treated either with 0.1% FBS, 2 mM thymidine or 800 nM nocodazole for 16–24 h. Cells were collected using TrypLE Express, washed in PBS, resuspended in 70% ethanol and stored at −20 °C. For analysis, cells were washed in PBS and resuspended in 500 µl staining buffer (PBS containing 40 µg ml−1 PI, 100 µg ml−1 RNase A, 0.1% Triton X-100), incubated on ice for approximately 1 h and analyzed by flow cytometry. The same conditions were used for analysis on RT–qPCR.

Identification of apoptosis-related lncRNAs

Cells assigned to the G1 cell cycle phase were extracted; fitting to the squared coefficients of variations against the means of normalized gene expressions (reads per kilobase million (RPKM)) was performed using the R function glmgam.fit() (similar to the method presented by Brennecke et al.49). The cell-to-cell variability of genes was ranked and the top 75 apoptosis-related genes (GO:0043065) were used for PCA. Cell clusters were identified using the pam function of the R package cluster v.2.1.2.


RNA was extracted (QIAGEN RNeasy Mini Kit) followed by DNase treatment (Ambion DNA-free DNA Removal Kit). Equal amounts of DNase-treated RNA was used to prepare cDNA (SuperScript II or Maxima H Minus RT; Thermo Fisher Scientific) and oligo(dT)18 primer according to the manufacturer’s recommendations. Quantification was carried out with Power SYBR Green Master Mix (Thermo Fisher Scientific) on a StepOnePlus or ViiA 7 Real-Time PCR System (Applied Biosystems). The Delta-Delta Ct method was used to quantify relative expression levels (normalized to siControl/ASOControl treatments and Beta-actin unless stated otherwise). Sequences for oligonucleotides are provided in Supplementary Table 5a. Samples were required to have similar RNA content (on DNase treatment) and similar Ct values of the Beta-actin internal control (on RT–qPCR) to be included in the analysis.

Cloning and generation of lentiviral U6 expressed shRNAs

Single-stranded oligonucleotides with Nhe1/Pac1 overhangs (synthesized by Integrated DNA Technologies; Supplementary Table 5) were phosphorylated (T4 Polynucleotide Kinase; New England Biolabs), linearized (95 °C for 3 min on a PCR cycler) and annealed by slowly decreasing the temperature on the PCR cycler. The previously generated pHIV7-IMPDH2-U6 construct50 was digested by Nhe1/Pac1 restriction enzymes, dephosphorylated (Antarctic Phosphatase; New England Biolabs) and gel-purified (QIAquick Gel Extraction Kit). The annealed oligonucleotides were ligated into the Nhe1/Pac1 and the digested pHIV7-IMPDH2-U6 construct (T4 DNA Ligase; Thermo Fisher Scientific); integration of shRNAs was verified by colony PCR and Sanger sequencing (Eurofins Genomics).

Lentiviral stable cell lines

HEK293FT cells were transfected with pCHGP-2, pCMV-G pCMV-rev and pHIV7-IMPDH2-U6 (refs. 50,51) at a 1:0.5:0.25:1.5 ratio using Lipofectamine 2000 and PLUS Reagent (Thermo Fisher Scientific) in serum-depleted DMEM medium. Medium was changed approximately 6 h post-transfection to DMEM containing 10% FBS, 1% penicillin/streptomycin, 1% NEAA, 1% sodium pyruvate, 2 mM L-glutamine, 0.37% sodium bicarbonate (supplements purchased from Thermo Fisher Scientific) and 1× Viral Boost Reagent (Alstem Cell Advancements). The viral supernatant was collected approximately 48 h post-transfection, passed through a 0.45-µm filter (Sarstedt) and concentrated with PEG-it (System Biosciences) according to the manufacturer’s recommendations. NIH/3T3 cells were transduced using a low titer of lentiviral particles (<10% of transduced cells) and green fluorescent protein+ cells sorted at the CMB Core Facility (Karolinska Institutet).

Colony formation assay

For stable NIH/3T3 cell lines, cells were seeded at 500 cells per well (6-well plates). After 10–14 d, cells were washed in PBS, stained for 20 min with 0.5% Crystal Violet, washed in water and left to dry. For quantification, stained cells were resolubilized in 10% acetic acid solution and then the absorbance was measured.

For siRNAs, NIH/3T3 cells were seeded at 1,000–5,000 cells per well in 6-well plates. Transfection was carried out 24 h after seeding and the procedure described above was repeated.

siRNA and ASO knockdown

NIH/3T3 and primary cells were transfected using Lipofectamine RNAiMAX Reagent (Thermo Fisher Scientific) according to the manufacturer’s protocol. A final concentration of 10 nM siRNA and 10 nM ASO was used. Cells were transfected the day after seeding and sorted (for Smart-seq3) or RNA-extracted (for RT–qPCR) 72 h after transfection. Sequences, company names and catalog numbers for siRNAs and ASOs are provided in Supplementary Table 5b.

PI-annexin V staining

PI-annexin V staining was carried out using the Annexin-V-FLUOS Staining Kit (catalog no. 11858777001; Roche) according to the manufacturer’s protocol. MMC treatment was initiated 24 h after siRNA transfection and samples were analyzed on a BD FACSMelody Cell Sorter 48 h later.

Functional prediction of lncRNAs using allelic imbalance

Genes were first filtered for (1) ≥3 allelic read counts in ≥20 cells, (2) not imprinted, (3) not encoded on the X chromosome and (4) having one of the following Ensembl BioMart annotations (GRCm38.p6, Supplementary Table 2): protein_coding; lncRNA; pseudogenes; transcribed_processed_pseudogene; transcribed_unitary_pseudogene; unitary_pseudogene; unprocessed_pseudogene; and transcribed_unprocessed_pseudogene.

Allelic imbalance of gene expression was measured as defined previously: (CASTallelicCounts / (CASTallelicCounts + C57allelicCounts) – 0.5). The allelic score (allelicImbalancelncRNA + allelicImbalancemRNA – diff(allelicImbalancelncRNA, allelicImbalancemRNA)) was calculated for each lncRNA-mRNA gene pair within 500 kb of the lncRNA TSS. The allelic score of the lncRNA-mRNA gene pairs was compared to a permutation test where each lncRNA (n = 542) was moved to 1,000 randomly selected mRNA gene positions. (The 1,000 genomic loci were kept the same for all lncRNAs and required to have at least 2 other genes in proximity.) The allelic score was computed for each lncRNA-mRNA gene pair over the randomly selected genomic loci (within ±500 kb pairs (kbp)) and P values were calculated as: allelicScorelncRNA:mRNA:random ≥ allelicScorelncRNA:mRNA:real / nrandomGeneInteractions.

Functional prediction of lncRNAs using allele-resolved RNA expression

Coordinated allelic expression of lncRNA-mRNA gene pairs (at the single-cell level) was addressed for all lncRNA-mRNA gene pairs within ±500kb of the lncRNA TSS (n = 542 lncRNAs). The expression pattern for each gene pair (≥3 allelic read counts) was evaluated using Fisher’s exact test (PReal, Benjamini–Hochberg-adjusted). To estimate the background, each lncRNA was translocated to 1,000 randomly selected gene locations and a Fisher’s exact test applied for all randomly generated gene pairs (PRandom, Benjamini–Hochberg-adjusted). lncRNA-mRNA gene pairs were considered significant if PReal < 0.01 where PRandom < PReal occurred in less than 1% of the permutated gene interactions.

Estimation of RNA half-lives and decay rates

Primary mouse tail fibroblast explants (F1 offspring from one adult female CAST × C57 and one adult female C57BL6, both in technical duplicates) were treated with actinomycin D (catalog no. SBR00013-1ml; Sigma-Aldrich) at a final concentration of 5 µG ml−1 in quadruplicate. RNA was extracted and global levels of RNA measured by poly(A)+ RNA-seq. Briefly, approximately 60 ng of DNase-treated RNA was prepared for sequencing using the Smart-seq2 protocol (modified for bulk RNA-seq) and sequenced on an Illumina NextSeq 500 (High-Output Kit v.2.5, 75 cycles). Data were processed using the zUMIs v.2.9.3e pipeline and genes filtered for ≥10 read counts in all 4 samples in the untreated condition (t0). Using RPKMs, gene expression was first normalized to the untreated condition (setting t0 = 1) for each individual sample. To normalize expression over the actinomycin D-treated time points, we took advantage of previous estimates of RNA half-lives in mouse embryonic stem cells22. We identified a subset of control genes with half-life estimates 1 h < t1/2 < 8 h with ≥50 read counts at t0 in all 4 actinomycin D-treated samples. The expected expression level of the control genes was calculated (y = 1 × exp(−kcontrol: × t)) and used to compute a ‘normalization factor’ (by taking the median) for each time point and sample, to which all genes were normalized to reach the final relative expression levels. Genes with shorter half-lives than 2 h were excluded from the 7 h and 10 h time points when calculating the ‘normalization factor’.

To estimate the half-lives, the normalized expression was fitted to an exponential decay curve (y = a × exp(−kx)) using the R package drc v.3.0.1. The decay rate (λ) was calculated using the formula: t1/2 = ln(2) / λ. Genes with half-lives <10 h and burst duration <72 h were considered for downstream analysis.

Statistical test for burst inference

To test the hypothesis regarding changes in burst kinetics, we used the likelihood ratio test. The test statistic for this test is essentially the difference between the likelihood of the null hypothesis (no change) and the likelihood of the observed change. Expressed as a formula, it is:

$$lambda _{mathrm{LR}} = – 2left[ {lleft( {theta _0} right) – lleft( {hat theta } right)} right]$$

Where λLR is the likelihood ratio test statistic, l(θ0) is the maximal log-likelihood where the null hypothesis is true, and (lleft( {hat theta } right)) is the log-likelihood of the maximized likelihood function (that is, the observed change). According to Wilk’s theorem, λLR converges asymptotically to the chi-squared distribution under the null hypothesis. This enables hypothesis testing of burst kinetics by comparing λLR to the chi-squared distribution with 1 d.f. At α = 0.05, the critical value is 3.84 for a one-sided test and 7.68 for a two-sided test.

In the context of burst kinetics, we focused on the log-ratio between, for example, burst frequency in the two samples. We set the null hypothesis θ0 = 0 and the alternative hypothesis (hat theta = log _2frac{{k_{on_2}}}{{k_{on_1}}}) where kon1and kon2 are the maximum likelihood estimates for both samples, respectively.

Simulations of burst inference

Simulations of burst inference were used to estimate the spread in inferred kinetics to be expected, given that the observed changes in expression were only caused by changed burst frequency or size, respectively. To evaluate the spread of changed burst frequency, we first modified the burst frequency by the observed change in mean RNA expression (assuming it is 100% explained by frequency); then, we simulated RNA count observations from the Beta-Poisson model (that is, the two-state model) with the same number of cells as present in the experiment. Then, we inferred the kinetic parameters; the densities of inferred parameters were shown as clouds in the ‘burst kinetics parameter space’. The rationale is that an alteration exclusively caused by any of the parameters would be expected to occur in these subsets of space, to guide intuition and further support the hypothesis testing performed above.

Statistics and reproducibility

No statistical method was used to predetermine sample size. The experiments were not randomized. The investigators were not blinded to allocation during the experiments and outcome assessment.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.


Source link