A multi-scale model of hepcidin promoter regulation reveals factors controlling systemic iron homeostasis

Abstract

Systemic iron homeostasis involves a negative feedback circuit in which the expression level of the peptide hormone hepcidin depends on and controls the iron blood levels. Hepcidin expression is regulated by the BMP6/SMAD and IL6/STAT signaling cascades. Deregulation of either pathway causes iron-related diseases such as hemochromatosis or anemia of inflammation. We quantitatively analyzed how BMP6 and IL6 control hepcidin expression. Transcription factor (TF) phosphorylation and reporter gene expression were measured under co-stimulation conditions, and the promoter was perturbed by mutagenesis. Using mathematical modeling, we systematically analyzed potential mechanisms of cooperative and competitive promoter regulation by the transcription factors, and experimentally validated the model predictions. Our results reveal that hepcidin cross-regulation primarily occurs by combinatorial transcription factor binding to the promoter, whereas signaling crosstalk is insignificant. We find that the presence of two BMP-responsive elements enhances the steepness of the promoter response towards the iron-sensing BMP signaling axis, which promotes iron homeostasis in vivo. IL6 co-stimulation reduces the promoter sensitivity towards the BMP signal, because the SMAD and STAT transcription factors compete for recruiting RNA polymerase to the transcription start site. This may explain why inflammatory signals disturb iron homeostasis in anemia of inflammation. Taken together, our results reveal why the iron homeostasis circuit is sensitive to perturbations implicated in disease.

Figure 1. Signal integration at the hepcidin promoter.

(A) Schematic representation of two critical signaling pathways controlling hepcidin expression. SMAD and STAT transcription factors are phosphorylated upon BMP and IL6 stimulation, and bind BMP-responsive elements (BRE) and a STAT-binding site (STATBS) in the hepcidin promoter, respectively. The importance of signaling crosstalk is not clear. (TSS: transcription start site; RNAP: RNA polymerase) (B) and (C) Analysis of transcription factor crosstalk at the promoter level by reporter gene assays. Luciferase expression is driven by the wildtype (WT) hepcidin promoter (3 kb upstream of TSS) or promoter mutants lacking one of the transcription factor binding sites (panel B; BRE1m = BRE1 mutated; STATdel = deleted for STATBS). Luciferase activity of each reporter construct (shown on a log10-scale) was measured for increasing doses of IL6 and/or BMP (n = 6). All heatmaps represent the mean of at least four biological replicates (see Methods), and are given in the same arbitrary concentration units (D) Moderate inhibitory signaling crosstalk at the signaling level. Immunoblots against phosphorylated SMAD1/5/8 and STAT3 (Supplemental Fig. S2) were quantified by densitometric analysis. The data points and error bars represent mean and standard deviation of biological replicates (N = 2), respectively (see Methods). Lines are fits of the sigmoidal Hill equation (y = y_basal + y_max * Sⁿ/(Sⁿ + EC₅₀ⁿ), S…stimulus, y_basal…basal signaling activity, y_max…maximal pathway activation, EC50….half-maximal-stimulus, n…Hill coefficient). The fits with and without non-canonical stimulation (blue and green lines, respectively) solely differ in the y_max values.

Figure 2. Mathematical modeling of signaling and promoter crosstalk.

(A) Thermodynamic modeling of promoter states. Depending on the transcription factor concentrations, the hepcidin promoter may be occupied by pSMAD (bound to BRE1 or BRE2), pSTAT (bound to STATBS) and RNAP (bound to the transcription start site), alone or in combination, giving rise to 16 different promoter states. A central presumption of thermodynamic modeling is that all RNAP-bound states are capable of transcription initiation, while RNAP-less states are considered silent. (B) A model selection approach allows for the identification of protein-protein interactions on the promoter. Various model variants were tested for their ability to fit the data in Fig. 1C. The minimal model (model 1) assumes that each transcription factor independently activates RNAP (grey arrows), while more complex variants additionally take cooperativity among transcription factors into account (red arrows). Statistical criteria (Akaike information criterion, likelihood ratio test) indicate that model topology 4 is best suited to describe all data (see Methods, Supplemental Protocol S2). (C) and (D) Integrative crosstalk model simultaneously fits luciferase data and dose-response curves of transcription factor phosphorylation. The thermodynamic promoter model (topology 4 in panel B) was coupled to a simple signaling model describing inhibitory crosstalk between phospho-SMAD and phospo-STAT transcription factors (Supplemental Protocol S2). Solid lines in C represent model trajectories in comparison to experimentally measured data points (shown as mean +/− std). The simulated luciferase activities in D agree well with the corresponding data in Fig. 1C.

Figure 3. Verification of model predictions using double-mutant promoters.

(A) Schematic representation of double-mutant promoters which lack two transcription factor binding sites (cf. Fig. 1B). (B) Systematic analysis of transcription factor binding site deletion effects confirms cooperativity of BRE1 and STATBS. The impact of binding site deletions was calculated by taking the luciferase activity ratios of different promoters (indicated in the legend) and expressed as a log10-fold change (y axis). As expected for a system where both sites cooperatively enhance transcription, the fold-change upon a combined deletion of BRE1 and STATBS (red) is less than the product of the single deletion fold-changes (green and blue; see text). Data points are mean and standard deviation, and model predictions represent the range of measurement-compliant parameter sets, as derived from a parameter identifiability analysis (see Methods, Supplemental Protocol S2). Only BRE1 and STATBS (but not BRE2) contribute to expression upon IL6 stimulation. (C) and (D) Co-stimulation heatmaps of double mutant promoters reveal that BRE1 and BRE2 are functionally similar in the absence of STATBS. (C) Heatmaps of luciferase activity under co-stimulation conditions. (D) Two-dimensional projection of the BRE1mSTATdel and BRE2mSTATdel data in C (averaged over all IL6 concentrations). Data points are mean (panel C, bottom row) or mean +/− std (panel D) (n = 6). Model predictions were formulated as ranges based on a parameter identifiability analysis (see Methods, Supplemental Protocol S2), and show measurement-compliant parameter sets with highest and lowest predicted luciferase expression (top and middle rows in panel C; edges of shaded corridors in panel D). Data and model in D were normalized to basal luciferase expression in the BRE2mSTATdel construct. (E) Systematic analysis of transcription factor binding site deletion effects confirms buffering of BRE1 and BRE2 single deletions. Concepts similar to panel B. The fold-change upon a combined deletion of BRE1 and BRE2 (red bars) is higher than the product of the single deletion fold-changes (green and blue bars; see text). BMP stimulation conditions were considered to ensure that BRE1 and BRE2 both contribute to expression.

Figure 4. Systems properties of hepcidin expression.

(A) The presence of two BREs enhances promoter sensitivity towards BMP stimulation. Hepcidin expression (fold over basal) is shown as a function phospho-SMAD levels for the WT, BRE1m, and BRE2m promoter (phospho-STAT was assumed zero). The dashed lines indicate the maximal steepness of the WT dose-response. The grey corridor indicates range of phospho-SMAD levels in HuH7 cells. (B) and (C) Hepcidin expression is highly sensitive to BMP stimulation, and less sensitive to IL6. The luciferase activity (fold over basal) is plotted as function of the IL6 (blue) or BMP (red) concentration. Panel B shows simulations of the best-fit model, while panel C contains experimental data (n = 3–6) and fits of the Hill equation (solid lines). Dashed lines in C indicate the maximal steepness of the BMP response. (D) Extended mathematical model describing negative feedback regulation of iron blood levels by hepcidin in vivo. Iron blood levels (Fe_b) are controlled by influx and efflux reactions, and the iron influx rate is proportional to the intestinal iron concentration (species Fe_i). Iron blood levels control the BMP signaling pathway, and thus the expression of hepcidin, which in turn lowers the iron influx. Hepcidin expression regulation by IL6 and BMP was modeled using the best-fit crosstalk model (Fig. 2C and D; Supplemental Text S4). (E) Iron homeostasis requires two BMP-responsive elements and is abolished by inflammatory stimulation. Simulations of the extended model (panel D) show how iron blood levels respond to changes in the intestinal iron concentration. The model with a WT hepcidin promoter (blue solid line) shows relatively constant iron blood levels over a broad range of intestinal iron concentrations (‘homeostasis range’). Homeostasis is less efficient and the homeostasis range is narrower in model variants with BRE1m and BRE2m promoters, or if strong IL6 stimulation is assumed (see legend) The mutants are characterized by altered iron blood levels (reflecting iron overload and deficiency, respectively). (F) IL6 stimulation reduces the BMP sensitivity of the promoter. The best-fit model (Fig. 2C and D) was employed to simulate how increasing IL6 stimulation affects the BMP dose-response curve of the promoter. Dashed lines indicate the maximal slope in the absence of IL6. Grey corridor same as in A.