Kinetics of mRNA nuclear export regulate innate immune response gene expression

Abstract

The abundance and stimulus-responsiveness of mature mRNA is thought to be determined by nuclear synthesis, processing, and cytoplasmic decay. However, the rate and efficiency of moving mRNA to the cytoplasm almost certainly contributes, but has rarely been measured. Here, we investigated mRNA export rates for innate immune genes. We generated high spatio-temporal resolution RNA-seq data from endotoxin-stimulated macrophages and parameterized a mathematical model to infer kinetic parameters with confidence intervals. We find that the effective chromatin-to-cytoplasm export rate is gene-specific, varying 100-fold: for some genes, less than 5% of synthesized transcripts arrive in the cytoplasm as mature mRNAs, while others show high export efficiency. Interestingly, effective export rates do not determine temporal gene responsiveness, but complement the wide range of mRNA decay rates; this ensures similar abundances of short- and long-lived mRNAs, which form successive innate immune response expression waves.

Fig. 1. Determining the kinetics of post-transcriptional events.

A Following mRNA from transcription to degradation. First, the mRNA is attached to the chromatin (caRNA), then it is released to the nucleoplasm (npRNA) upon 3′-cleavage, finally it is transported to the cytoplasm where it is degraded. B Schematic of the experimental setup. Bone-marrow-derived macrophages (BMDMs) are stimulated with LPA. Cells are harvested at different time after stimulation and fractionated into subcellular fractions, and RNA is extracted. C Example of tracks for the Tnf gene. Many intronic reads are found in the chromatin-associated fraction, fewer in the nucleoplasmic fraction, and the cytoplasmic fraction RNA is fully spliced. D Gene selection workflow. Lowly expressed genes were filtered out, and inducible genes at the chromatin level were selected. E Heatmap of gene expression of the selected genes in (D) for each fraction. F Correlation analysis of mRNA abundance in each fraction. Note a stronger correlation between chromatin and nucleoplasmic fraction for nearby timepoints than between nucleoplasmic and cytoplasmic fractions or between cytoplasmic and chromatin fractions. Also, note a time shift, that timepoints on the chromatin correlate better with later timepoints in the cytoplasmic fraction. G Example genes: Top, Cd74 and Btg2 have similar chromatin expression profiles but exhibit different profiles in the nucleoplasmic and cytoplasmic fraction. Bottom, Arl5b and Ccr3 show similar cytoplasmic levels but Arl5b expression is few fold higher than Ccr3 in both chromatin and nucleoplasmic fractions.

Fig. 2. Fitting a kinetic model to the stimulus-response data.

A Model linking mRNA measurements in different subcellular fractions, with associated equations. B Schematic of the fitting workflow that takes the expression of caRNA as input, simulates the nucleoplasmic and cytoplasmic mRNA abundances, and iterates to identify the optimal parameter set for each gene. C Heatmap of fitting results alongside the data. Most genes fit the data really well, as indicated in the Fit Quality bar on the left-hand side. Source data are provided as a Source Data file. D Example of fitting results for Egr1 and Malt1 genes for the three replicates. The error bars represent the 95% prediction interval of the model given the data library size.

Fig. 3. Parameter inference, identifiability, and reproducibility.

A Number of genes for which parameters were identifiable, partially identifiable, or non-identifiable (left panel). Distribution of 95% confidence interval range for identifiable parameters. Also see the Supplementary Data file. Source data are provided as a Source Data file. B Distribution of the parameter values for the different replicates. The distributions are similar for the different replicates and the bulk of the distributions span a 30 to 100-fold difference between the genes depending on the parameter. C Reproducibility of each parameter. The color of the point corresponds to the fit quality metric for the worst replicate, the line corresponds to 95% confidence interval, with the color corresponding to the fit quality of that replicate. Parameters k₁’ and k₂ that were identifiable for both replicates are also quite reproducible, while k₂’ and k_{cyto_deg} are very reproducible and well-defined. The dashed red line indicates the 2-fold reproducibility window. D Reproducibility of composite parameter k₁’/k₂. Even though k₁’ and k₂ are well not defined for some genes, their ratio may be highly reproducible and well-defined. The dashed red line indicates the twofold reproducibility window. E Comparison of model-inferred half-life with half-life values determined with the actinomycin-D method. The dashed lines link replicates of the model-inferred half-life. The green line indicates a 1:1 relationship and the dashed red line a 2-fold range. We notice a high correlation between the two half-life estimates (Spearman rank correlation of 0.8) but that the model seems to capture a larger range in half-lives than the actinomycin-D method.

Fig. 4. Export rates vary widely across immune response genes.

A Left panel shows the histogram of nucleoplasmic export efficiencies (k₂’/k₂) for immune response genes. The distributions are similar for the different replicates and span a 30-fold difference between genes. Middle panel shows the reproducibility of these export efficiency estimates. Right panel shows the number of genes for which this quantity is identifiable. B Left panel shows the histogram of effective transport rates (k₁’k₂’/k₂) for immune response genes. The distributions are similar for the different replicates and span a 100-fold difference between genes. Middle panel shows the reproducibility of these effective transport rate estimates. Right panel shows the number of genes for which this quantity is identifiable. C Correlation of the nucleoplasmic export efficiency and the effective transport rate for the genes for which both parameters are identifiable. D The distributions of effective transport rates of genes involved in various biological processes. E Example expression for two genes having effective transport rates (k₁’k₂’/k₂) on the extremes of the distribution. Malt1 has a low effective transport rate and Egr1 has a high rate. From the line graph (top), we noticed that even though Egr1 is less expressed than Malt1 in the chromatin fraction, its expression ends up being higher in the cytoplasmic fraction.

Fig. 5. Potential determinants of the effective mRNA transport rate.

A Effective transport rates (k₁’k₂’/k₂) show a significant negative correlation with gene length (Pearson’s correlation of −0.47 and −0.42 with associated P value of 3e-12 and 9e-10 for replicate 1 and 2, respectively). B Effective transport rates show a significant negative correlation with the number of introns (Pearson’s correlation of −0.54 and −0.53 with associated P value of 3e-16 and 8e-16 for replicate 1 and 2 respectively). C Effective transport rates show a significant positive correlation with splicing probability, averaged over all introns (Pearson’s correlation of 0.63 and 0.55 with associated P value <2.2e-16 and of 1e-13 for replicate 1 and 2, respectively). This correlation coefficient is higher than when only the most retained intron is considered (Supplementary Fig. S5A). D Sum of ChIP-seq signals of indicated histone mark associated with the gene do not show a correlation with the effective transport rate (alternatively, windows different sizes along the gene, described in Methods, were tried but yielded no better correlation). E Machine-learning models reveal little predictive power in histone modification ChIP-seq signals. Top, plot of R² values that indicates the predictive power of machine-learning models that consider indicated features. Error bars indicate the mean +/− standard deviation of R² value of the cross-validation sets. A two-sided t test was used with */**/*** indicating a P value of <0.05, <0.01, <0.001. Bottom, heatmap of the features’ importance, defined by the gain in accuracy brought by each feature normalized by the total gain. Numbers in the heatmap correspond to the number of ChIP-seq bins selected by the model.

Fig. 6. Prior LPS exposure diminishes transcriptional initiation but has little effect on export rates.

A Comparison of observed basal expression (top) and max fold change (bottom) at the chromatin level. Tolerization increases basal expression and decreases induced expression for most genes. B Heatmap of fold change in the naive condition and tolerized condition. We observe in all fractions that the induction is dimmed in the tolerized condition. C Distribution of fit quality for tolerized and naive macrophages. The fit quality for the tolerized macrophage is not as good as for the naive (heavier right tail). Source data are provided as a Source Data file. D Comparison of post-transcriptional parameters in naive and tolerized cells. We observe that k₁’ and k₂ are not well-defined, as in naive macrophages. However, k₂’ and k_{cyto_deg} are well-defined and very similar between to those in naive condition. E Comparison of the composite parameters, effective transport rate (k₁’k₂’/k₂), and export efficiency (k₂’/k₂), between naive and tolerized conditions. These parameters have similar values for most genes, even if less well-defined in the tolerized condition.

Fig. 7. Understanding the link between responsiveness, transport, and mRNA half-life.

A Simulations of gene induction with different parameter values for mRNA decay (with indicated half-lives, HL) and effective nuclear transport rates ranging from 0.03 to 0.3 cytoFPKM/caFPKM min⁻¹ as observed for the fitted genes. Responsiveness is measured as time to half induction. B Heatmap of responsiveness for various k₂’ and k_{cyto_deg} parameter values. There are three main regimes, in the top left corner, the responsiveness is determined solely by half-life, in the bottom right corner the responsiveness is solely defined by k₂, and the diagonal may depend on both. The parameter valued for the immune response genes are overlaid on top. Most genes fall in the region where responsiveness is solely determined by the cytoplasmic decay rate. C Scatterplots to show the correlation of gene responsiveness to a step function with half-life (left panel) and transport (right panel). A strong relationship between half-life and responsiveness is observed. There is a weaker, nonlinear relationship with the transport rate. D Scatterplot to show the correlation between the chromatin-to-cytoplasm transport rate and the cytoplasmic decay rates. E Heatmap of relative cytoplasmic mRNA abundance as a function the effective transport and cytoplasmic degradation rate. F Scatterplots between the peak cytoplasmic expression level of each gene with either the effective transport rate or the half-life of its mRNA. This confirms that neither quantity is correlated with the expression level, supporting the model that the two balance each other to render the level of expression controlled by other mechanisms, such as transcriptional initiation.