Protein–DNA binding dynamics predict transcriptional response to nutrients in archaea

13 Pages·2013·2.19 MB·English

8546–8558 Nucleic Acids Research, 2013, Vol. 41, No. 18 Published online 26 July 2013
doi:10.1093/nar/gkt659
Protein–DNA binding dynamics predict
transcriptional response to nutrients in archaea
Horia Todor1, Kriti Sharma1, Adrianne M. C. Pittman1 and Amy K. Schmid1,2,*
1Department of Biology, Duke University, Durham, NC 27708, USA and 2Center for Systems Biology, Institute
for Genome Science and Policy, Duke University, Durham, NC 27708, USA
Received May 24, 2013; Revised July 2, 2013; Accepted July 4, 2013 D
o
w
n
lo
a
ABSTRACT andmetabolicregulatorynetworksattheenzymatic,tran- de
d
Organisms across all three domains of life use gene sspcroipntdiionngaltaondflpuoctsut-attrianngscrciopntidointiaolnlseve(1ls).inFoorrganexisammsprlee-, from
regulatory networks (GRNs) to integrate varied transient nutrient changes may result in microsecond- h
stimuliintocoherenttranscriptionalresponsestoen- scale regulation of enzyme activity, whereas prolonged ttps
vironmental pressures. However, inferring GRN exposure during the course of minutes or hours may ://a
c
topology and regulatory causality remains a central trigger changes in the levels of enzyme-coding transcripts a
d
challenge in systems biology. Previous work (2).Becauseofthepotentialforthebuildupoftoxicinter- em
characterized TrmB as a global metabolic transcrip- mediates and futile cycles, temporal dynamics of meta- ic.o
tion factor in archaeal extremophiles. However, it bolic enzyme-coding transcripts can be as important as up
remains unclear how TrmB dynamically regulates overall levels (3). Over evolutionary time scales, the .co
m
its (cid:2)100 metabolic enzyme-coding gene targets. constant presence of a given nutrient may lead to gene /n
loss in competing metabolic pathways. Such streamlining a
UdastiengthaedtyonpaomloigcypoerftuthrbeatTiormnBapmpreotaabcho,licweGReNluciin- olofntgh-etegrmenoamdaepitsatthioonusghint tnoetewnoabrklesftarsutcetrurreepolircattoiopnolaongdy r/article
the model archaeon Halobacterium salinarum. (4). Discerning the dynamic function of the underlying -a
b
Clustering of dynamic gene expression patterns network resulting from the interaction of many different stra
reveals that TrmB functions alone to regulate levels of regulation remains a central challenge. c
central metabolic enzyme-coding genesbutcooper- Inarchaea,evidenceismountingthattranscriptionmay t/41
ates with various regulators to control peripheral be an important mechanism for regulating metabolism. /18
metabolic pathways. Using a dynamical model, we Unlike eukaryotes and bacteria, in which regulation of /85
4
predictgeneexpressionpatternsforsomeTrmB-de- flux through central carbon metabolism and other 6/1
pathways appears to be allosteric (5), archaea seem to 0
pendent promoters and infer secondary regulators 3
forothers.Ourdatasuggestfeed-forwardgeneregu- lack many classic allosteric regulatory control points 554
(6).For example, in hypersaline-adapted archaea, glutam- 2
latory topology for cobalamin biosynthesis. In b
ate dehydrogenase was found to be unresponsive to ADP y
contrast, purine biosynthesis appears to require g
and GDP (7), whereas D-Lactate dehydrogenase was not u
e
TTrrmmBB-iisndanepimenpdoerntatntrcegoumlaptoonres.ntWfoermceodniactliundgemethtaa-t rpeygruulvaatetedkbinyafsreufcrtoomse-H1,a6l-obbiaspctheorsiupmhastaeli(n8a)r.uAmdwdiatsiofnoaulnlyd, st on
1
bolic modularity, integrating nutrient status and to possess only weak allosteric regulation (9). Other po- 6
N
regulating gene expression dynamics alone and in tential regulatory mechanisms such as protein phosphor- o
v
concert with secondary regulators. ylation have been proposed in two different archaeal em
species (10,11); however, in the hypersaline-adapted b
e
archaeal model system H. salinarum, much of the r 2
INTRODUCTION 0
missingregulationseemstotakeplaceattheleveloftran- 1
8
Diverse metabolic processes must be differentially scription(12–15).Thispropertyofhaloarchaeaprovidesa
regulated to maintain homeostasis and optimize growth simplified model system for understanding the underlying
in changing environmental and intracellular conditions. logic of the integrated transcriptional and metabolic
Environmental fluctuations occur at several temporal network. The largely transcriptional nature of archaeal
scales. This is mirrored in the integrated transcriptional responses has already been used to infer regulatory
*To whom correspondence should be addressed. Tel:[phone]; Fax:[phone]; Email: [email protected]
Present address:
Amy K. Schmid, Department of Biology, Duke University, Durham, NC 27708, USA.
(cid:2)TheAuthor(s)2013.PublishedbyOxfordUniversityPress.
ThisisanOpenAccessarticledistributedunderthetermsoftheCreativeCommonsAttributionLicense(http://creativecommons.org/licenses/by/3.0/),which
permitsunrestrictedreuse,distribution,andreproductioninanymedium,providedtheoriginalworkisproperlycited.
NucleicAcidsResearch,2013,Vol.41,No.18 8547
networks with a remarkable degree of accuracy (14). MATERIALS AND METHODS
H. salinarum survives in an extreme and fluctuating
Strains and growth conditions
environment, where daily and seasonal changes in
salinity,oxygenandnutrientsrequireconstantadjustment H. salinarum NRC-1 (ATCC strain 700922) was used as
of metabolism (16). Current evidence suggests that the wild-type strain background for all studies
many of these metabolic adjustments are regulated (Supplementary Table S1). Gene expression was assayed
transcriptionally. in a previously constructed strain containing an in-frame
For example, in H. salinarum, TrmB has been deletion of VNG1451C [(cid:2)ura3DtrmB; (17)] and its
characterized as the central transcriptional regulator of isogenic parent strain ((cid:2)ura3). A previously constructed
carbon metabolic pathways under aerobic conditions strain containing trmB::c-myc fusion on a low copy
(17). Conserved throughout the archaea, TrmB is a number plasmid was used for ChIP-qPCR to determine
D
winged helix-turn-helix transcription factor (TF) that bindingsiteoccupancy(17).Cellsforgeneexpressionwere o
w
binds a palindromic inverted repeat cis-regulatory growninCompleteDefinedMediumcontaining19amino n
sequence in Pyrococcus furiosus (18), Thermococcus acids [modified from (17); Supplementary Table S2] sup- loa
d
litoralis(19)andH.salinarum(17).TrmBactsasarepres- plemented with 50mg/ml uracil to complement the ura3 e
d
soratsomepromotersandasanactivatoratothers.Inthe deletion. Strains carrying the trmB::c-myc construct for fro
obligately anaerobic hyperthermophillic archaea, TrmB ChIP-qPCR were grown in Complete Defined Medium m
h
has been characterized as a regulator specifically supplemented with 20mg/ml mevinolin for plasmid main- ttp
i(n2v0o).lvIendcoinntroalsitg,oTsarcmcBhairnidHe .trsaanlisnpaorurtmahnads bceaetnabfooluisnmd t4e2n(cid:3)aCncuen.dCerullotuwreasmwbeierentglrigohwtn. at 225 r.p.m. shaking at s://ac
a
to bind 113 promoters in the genome to regulate genes d
e
involvedindiverseprocessesincludingcentralcarbonme- mRNA preparation for gene expression time course m
ic
tabolism, TCA cycle, amino acid metabolism and co- .o
H. salinarum Dura3 (parent) and (cid:2)ura3DtrmB (knockout u
factor metabolism (17). p
mutant) strains were grown to early logarithmic phase .c
Yet even in these simplified systems, where metabolic o
(OD600 &0.3). For gene expression time courses, three m
flux should largely be predictable from transcriptional /n
4ml of aliquots were removed from the continuously a
data alone, regulatory causality is complicated by meta- shaking cultures before the addition of 5% (w/v) glucose r/a
bolic feedback. For example, TrmB appears to be ((cid:4)240, (cid:4)60, 0min time points) and seven afterwards (5, rtic
regulated at the level of TF activity (14,17), binding le
10,20,45,90,180,360mintimepoints).Cellswereimme- -a
sugars(glucose,trehalose,maltose)withvaryingaffinities, b
which in turn decrease TrmB affinity for DNA (17–19). diatelypelletedbycentrifugationandsnapfrozeninliquid stra
nitrogen. RNA was prepared using the Absolutely-RNA c
However,byde-activatingtranscriptionofppsA(encoding miniprep kit (Stratagene/Agilent) according to the manu- t/4
phosphoenolpyruvate synthase) and other gluconeogenic facturer’s instructions. RNA quality was assessed using 1/1
genes while de-repressing transcription of pykA the Agilent 2100 BioAnalyzer RNA-Nano Chip and the 8/8
(encoding pyruvate kinase) and other glycolytic genes, 54
Prokaryotic Nano RNA protocol (Agilent Technologies, 6
the amount of glucose is decreased, causing increased /1
Santa Clara, CA), and samples were verified to be free of 0
TrmB-promoter binding and reversing the effect (17). In 3
DNAcontaminationby25cyclesofPCRamplificationon 5
bacteria, the enzymes involved in glycolysis and atleast200ngofRNAbeforecDNAconversioninreverse 542
gluconeogenesis are regulated allosterically (5) in transcriptase quantitative PCR (RT-qPCR). by
response to transient changes, and transcriptionally g
u
through cAMP Responsive protein (Crp in gram e
RT-qPCR s
negative bacteria) in response to more permanent t o
n
changes(21).Incontrast,wehypothesizethatTrmBregu- RNAextractedfromcellswasquantifiedusingRT-qPCR. 1
6
lates metabolic enzyme-coding genes rapidly in response Briefly, RNA was quantified using the Power SYBR N
to stimuli while also adapting the equilibrium levels of RNA-to-CT 1-step kit (Applied Biosciences) in an ove
these genes to longer-term conditions. Eppendorf Mastercycler ep Realplex thermocycler using m
b
Regulation of diverse pathways in response to environ- Realplex software. Reaction size was 20ml, and reactions e
mental fluctuations of varying temporal scales is also were prepared according to manufacturer’s instructions. r 20
1
likely to be mediated by several TFs working together in Plate setup was performed robotically using a Corbitt 8
various combinations to produce appropriate transcrip- Life Sciences liquid handling system. C threshold was
Q
tional dynamics for each pathway (14,22). Such network determined automatically by the software. Primers were
topology is difficult to infer on the basis of steady state crosschecked against the H. salinarum genome using
measurements of gene expression and promoter occu- BLAST to ensure specificity. Amplification efficiencies
pancy. To unravel cause and effect relationships in the for each pair of primers (ppsA – VNG0330G, pykA –
gene regulatory network (GRN) controlling metabolism, VNG0324G, eif1a2 - VNG1756G) were determined by
wemeasuredgeneexpressionusingNanoString(cid:3) (23)and running three serial 10-fold dilutions of the same sample
TrmB binding dynamics using ChIP-qPCR over time in (see Supplementary Table S3 for primers). These
responsetoaglucosestimulus.Integrationofthesedatain efficiencies were used to calculate enrichment relative to
the context of a dynamical model enabled prediction of thereferencegene(VNG1756G)inRNAfromtimecourse
how TrmB and its partners enact several transcriptional sampling from the measured C values as previously
Q
programs in different ways across metabolic pathways. described (24). Each experiment consisted of two to
8548 NucleicAcidsResearch,2013,Vol.41,No.18
three biological replicates (separate cultures) and three point, and for each gene, fit an exponential decay curve
technical replicates. The significance of change in gene through the points that were annotated as valid in the
expression before and after nutrient addition was microarray file using the nls() function in the statistical
assessed using the Welch one-sided t-test. Using the first packageR.Standarderrorwasalsocalculatedusingnls().
three and last two points of each biological replicate time
course,wedeterminedP-valuesforthedifferenceinmeans
Modeling approach and fitting
equivalent to a 1.5-fold up- or downregulation.
Synthesis rate determination. To better understand how
Measurement of gene expression using NanoString binding events at the promoter influence gene expression,
we calculated thederivative with respectto time forevery
Quantification of additional genes was performed on the
gene at each time point by the weighted average of the
same RNA samples using NanoString technology (23).
slope of the segments before and after each time point. D
One hundred micrograms of each RNA sample from WeusedtheK valuesfrom(29)tocalculatethesynthe- ow
g(Sluecaottslee, tWimAe),cowuhrseerse swaemrepledseliwveerreedhytobridNizaendoSttoringa sis rate accordidnegg to Formula 1. nload
e
customprobeset(SupplementaryTableS4)encompassing d½mRNA(cid:5) d
100 genes of interest and counted using an nCounter dt +Kdeg½mRNA(cid:5)¼Ksynthesis ð1Þ from
machine. The same mRNA samples used for RT-qPCR h
were used for NanoString experiments. Data were Gene expression prediction. To determine the import- ttp
s
normalized by total counts per sample across strains. All ance of TrmB regulation to specific genes, we used an ://a
raw and normalized data are included in Supplementary ordinary differential equation (ODE) model to predict ca
Table S5. For clustering, each gene was further mean gene expression dynamics in both WT and (cid:2)ura3DtrmB de
m
centered and normalized to a standard deviation of 1. strains as a function of ChIP-qPCR enrichment at the ic
Gene expression profiles were clustered using k-means ppsA promoter. The first three and last three points of .ou
with a Pearson correlation distance function on the the Dura3 gene expression and TrmB-binding time p.c
o
normalized gene expression data with k=8. The same courses were used to determine the values for Kbasal and m
clustering parameters were used for both (cid:2)ura3 parent Keff by setting the right hand side of Formula 2 to 0. The /na
and (cid:2)ura3DtrmB mutant data sets. Each cluster was Hill coefficient (n) was set to either negative or positive r/a
analyzed for enrichment of Clusters of Orthologous one to indicate repression or activation, respectively. The rtic
le
Groups (COG) categories (25,26) using the level of gene expression during the nutrient stimulus was -a
b
hypergeometric test. predicted by calculating the change in gene expression s
every minute according to Formula 2. tra
c
ChIP-qPCR protocol t/4
d½mRNA(cid:5) 1
H. salinarum cells harboring the trmB::c-myc construct dt ¼(cid:4)Kdeg½mRNA(cid:5)+Kbasal+Keff:TrmBnenrichment /18/8
were grown to early logarithmic phase (OD600 &0.3) as 5
ð2Þ 4
described earlier in the text. DNA–TF complexes were 6
/1
cross-linked and immunoprecipitated using the c-myc Calculating model residuals. To elucidate the topology 03
epitope tag as previously described (17). Primers were of nutritional control of central and peripheral metabol- 55
4
designed according to the criteria presented in (27). The ism, we compared our predicted and actual gene expres- 2 b
qPCR thermocycling reaction and thermocycling condi- sion data. Each gene was normalized to a maximum of 1. y g
tions were as described previously (28), except that the The normalized values were used as input to the TrmB- ue
SsoAdvanced SYBR Green Supermix (Bio-Rad) was enrichment-based predictive ODE model of gene expres- st o
used. Each of five biological replicates was run in three sion (earlier in the text). For each gene in each of the n 1
technical replicates. Enrichment of TrmB binding at the (cid:2)ura3DtrmB and (cid:2)ura3 strain backgrounds at each 6 N
ppsA promoter was calculated as relative enrichment of mRNA time point, we calculated the difference between o
v
the immunoprecipitated sample versus the input. The thepredictedsynthesisrateandtheactualsynthesisrateas em
ratios given in the figures compare enrichment at the calculated earlier in the text (see ‘Synthesis Rate be
binding peak to the 30 end of the gene of interest (28) Determination’ section). We clustered resulting traces on r 2
0
(Supplementary Table S3). the basis of the (cid:2)ura3 data using k-means clustering with 18
k=5 and Euclidean distance.
Degradation constant calculation
Feed-forward loop logic approximation. To assess
To assess the validity of our data and improve the whether a feed-forward loop (FFL) might be responsible
accuracy of our modeling, we decided to incorporate forthedynamicsweobservedinthecobalaminbiosynthe-
genome wide mRNA degradation parameters. Although siscluster,wesimulatedthenetworkusingalogicapproxi-
theseparametershavebeendeterminedforarelatedstrain mation (Figure 8A) (30). Genes involved in cobalamin
[H.salinarumR1ATCC29341/DSM671;(29)],theactual biosynthesis were normalized so that the lowest value
half-liveswerenotpublished.Wethereforereanalyzedthe was 0 and the highest value was 1. A regulator was
microarray data (ArrayExpress accession number E- presumed active when its value was >0.5. We used the
MEXP-1088, http://www.ebi.ac.uk/arrayexpress/) to de- TrmB ChIP-qPCR binding data as an input and fit the
termine these values (Supplementary Table S6). In short, degradation parameter (K & K ) for
deg FFL Effector deg.cob/cbi
we normalized each gene on each array to the t=0 time both of the other genes to the average scaled expression
NucleicAcidsResearch,2013,Vol.41,No.18 8549
profile of the cob/cbi cluster using least squares. The fit point in glycolysis/gluconeogenesis (6), state-change
was evaluated at 5, 10, 20 and 45min following glucose dynamics in response to glucose and glycerol were
addition. This is the period in which feed-forward expected. However, TrmB regulates genes involved in
dynamics would be expected to be most apparent. processes across metabolism (17). To assess the impact
of TrmB on gene expression in response to glucose, we
assayed mRNA levels over time following glucose
RESULTS
stimulus for 100 genes involved in central and secondary
Gene expression dynamics in response to nutrients are metabolism using NanoString (23). NanoString was used
dependent on TrmB to assay gene expression because it has been shown to be
sensitive, preciseandrequire minimalprocessing (33).We
Our previous work demonstrated that TrmB binds the
selected96genesinvolvedincentralmetabolism,including
promoter of metabolic enzyme-coding genes throughout D
both direct and indirect TrmB targets and TrmB itself o
thegenomeinresponsetocarbonsource availability(17). w
(Supplementary Table S4). Using data from previous n
However,thesestudieswereconductedatsteadystate.To lo
microarray experiments (13,14,34,35), we also selected a
investigate the dynamic expression response of TrmB- three genes (VNG1670C, mrp, srp19) expressed at high, de
d
regulated genes to nutrients (glucose, glycerol and medium and low levels with minimal variation across fro
sucrose), we used RT-qPCR to measure repressed (pykA) many conditions (Supplementary Table S4). The same m
and activated (ppsA) gene levels in response to nutrients h
othveesretgimenees(Fbiygbuirnedsi1ngantodt2h)e.pTrrommBotiesrtehitohuegrhttotaoctrievgautelaoter g(AFlusigcNousraeeno1S)retrwsipneorgenhsseaasmntpiomlteedbeienpnotihpnertesvNioauunsseoldSytruifsnoegrdefxoRprTemr-iqmRPeNCnRAt. ttps://ac
torepressexpressionintheabsenceofglucose.Additionof a
quantification in H. salinarum, we validated the data by d
e
glucose to the medium results in TrmB dissociation from comparing measurements of ppsA and pykA gene expres- m
the promoter and de-activation or de-repression of the sion in the same RNA time course samples between the ic.o
targetgene(17).Briefly,H.salinarumcellsweregrownon NanoString and RT-qPCR platforms. We found a linear up
aminoacidsasacarbonandenergysourcetomid-logarith- relationship (Pearson correlation=0.962) across several .co
m
micphase.Cellsweresampledthricebeforeandseventimes orders of magnitude (Figure 3), suggesting that /n
after the addition of nutrients (‘Materials and Methods’ NanoString is a robust alternative to RT-qPCR for the ar/a
section).Weconsideredachangerelevantwhena1.5-fold quantification of mRNA in H. salinarum. rtic
or greater up- or downregulation of the target gene was TodeterminethepatternofTrmB-dependentregulation le
significant (P<0.05). As these genes are a key control -a
across central and secondary metabolic pathways, b
s
point in glycolysis (pykA) and gluconeogenesis (ppsA), NanoString gene expression profiles were clustered separ- tra
their levels are highly informative of the regulation of c
ately in the (cid:2)ura3 and (cid:2)ura3DtrmB strain backgrounds t/4
thatpathway(31,32). usingk-means.Thedistributionofclusterswasintegrated 1/1
WeobservedthatbothppsAandpykAmRNAexhibita with the TrmB-centered metabolic regulatory network re- 8/8
characteristic steady-state expression value during mid- construction(17)(Figure4andSupplementaryTableS5). 54
logarithmic phase growth before glucose addition that 6
ClustersofgenesfrombothDura3and(cid:2)ura3DtrmBexpres- /1
differs between the (cid:2)ura3 parent and (cid:2)ura3DtrmB 0
sionprofileswereanalyzedforsignificant(P<0.01)enrich- 3
mutant strain (Figure 1). On addition of glucose, ppsA is ment of membership in COGs (Table 1) (25, 26). 554
de-activated(Figure1A,3.6-folddecrease)andpykAisde- 2
Surprisingly, we found that temporal gene expression b
repressed (Figure 1B, 12.7-fold increase) by the first time y
patterns clustered according to metabolic pathway g
point (5min), and both genes reached a new equilibrium u
modules in the Dura3 strain and to a lesser extent in the e
level by 45min. Hereafter, we refer to such monotonic (cid:2)ura3DtrmBmutant,suggestingthatTrmBisatleastpar- st o
changes in gene expression as state-change dynamics. n
tiallyresponsibleforthemaintenanceofmetabolicmodu- 1
These dynamics were greatly attenuated in the 6
larity (Figure 4C and F and Supplementary Table S5 and N
(cid:2)ura3DtrmB mutant background in response to glucose Table1). ov
(Figure 1) and in the (cid:2)ura3 parent strain in response to The gluconeogenic enzyme-coding genes exhibited the em
sucrose (Figure 2). Gene expression dynamics similar to b
same state-change downregulation dynamics as ppsA in e
those in response to glucose were observed with glycerol the (cid:2)ura3 strain, with the exception of enolase (eno), r 20
(Figure 2), although ultimately, a different final equilib- 1
which is not a direct TrmB target (17). These dynamics 8
rium level was reached. These data confirm that gene ex-
were not observed in the (cid:2)ura3DtrmB strain (Figure 4A
pression dynamics of ppsA and pykA are specific to
and Table 1). Because state-change dynamics in
glucose and glycerol and dependent on TrmB.
gluconeogenic enzyme-coding genes were greatly dimin-
Furthermore, theyconfirmtheroleofTrmBinregulating
ished in the (cid:2)ura3DtrmB strain, it is likely that these
genes coding for enzymes in glycolysis/gluconeogenesis at
genes are predominantly dependent on TrmB. As
both short and intermediate temporal scales.
expected, the expression profiles of enzyme-coding genes
ingluconeogenesisclusteredinthe(cid:2)ura3parentstrainbut
NanoString measurement of temporal gene expression
not in the (cid:2)ura3DtrmB strain (Table 1).
reveals metabolic modularity governed by TrmB
Incontrasttothestate-changedynamics(i.e.monotonic
In our initial experiments, both ppsA and pykA gene increasesordecreasesandastatisticallysignificantchange
expression exhibited state-change TrmB-dependent in equilibrium gene expression) observed in the glycolytic
dynamics. Because these genes are an important control and gluconeogenic pathways, impulse-like dynamics
8550 NucleicAcidsResearch,2013,Vol.41,No.18
(A) (A)
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
(B) (B) ://a
c
a
d
e
m
ic
.o
u
p
.c
o
m
/n
a
r/a
rtic
le
-a
b
s
tra
c
t/4
1
/1
8
/8
5
4
6
/1
0
3
5
5
4
2
Figure 1. ppsA and pykA gene expression exhibits state-change b
Figure 2. ppsA (A) and pykA (B) do not respond significantly to 5% y
dynamics in response to glucose perturbation when TrmB is present. sucrose in the (cid:2)ura3 strain (black lines). ppsA and pykA respond to gu
GeneexpressionofppsA(A)andpykA(B)toa5%glucosestimulusin 0.167%glycerolstimulusinthe(cid:2)ura3strain(graylines).Geneexpres- es
the (cid:2)ura3 (black lines) and (cid:2)ura3DtrmB (gray lines) strains were sionwasmeasuredbyRT-qPCRandisshownplottedonalogarithmic t o
measured by RT-qPCR and are shown plotted on a logarithmic axis. n
Errorbarsrepresentthestandarderrorfromtheaverageofatleasttwo abxetisw.eAenstethrieskbseignindniciantgeasnigdnitfihceanecnedoofftthheedtiifmfeerecnocuersine;e*xpsrigesnsiifiocnanletvaetl 16
biological replicate experiments. Asterisks indicate significance of the N
P<0.05; ** significant at P<0.01; *** significant at P<0.001. o
differenceinexpressionlevelbetweenthebeginningandtheendofthe v
e
time course; * significant at P<0.05; ** significant at P<0.01; *** m
significant at P<0.001. be
TrmB-dependent,thedifferenceinoverallmRNAconcen- r 2
0
trationcomparedwiththeparentstrainatthestartofthe 18
(transient increases or decreases in gene expression) were timecoursebutnotattheendsuggeststhatTrmBisnone-
observed in the genes coding for enzymes in purine, co- theless important in controlling these genes.
balamin and amino acid biosynthesis in the Dura3 strain. Expression patterns in genes coding for enzymes
The variety of impulse-like dynamics made metabolic involved in cobalamin (vitamin B-12) biosynthesis also
modularity especially apparent in these clusters. The differed between the (cid:2)ura3DtrmB mutant and (cid:2)ura3
pattern of gene expression in the purine biosynthesis parent strains. In the (cid:2)ura3 parent strain, these genes
cluster was distinct in both the (cid:2)ura3 and (cid:2)ura3DtrmB were downregulated as soon as glucose was added. This
strains(Figure4BandGandTable1).Expressionofthese was followed by a transient upregulation centered 45min
genes displayed a transient increase in mRNA immedi- after the addition of glucose (Figure 4E). These distinct
ately on the addition of glucose in both backgrounds expression patterns were tightly correlated and clustered
(Figure 4B and G). Although the impulse-like dy- togetherintheDura3strain(Figure4E).Incontrast,inthe
namics observed in purine biosynthesis genes are not (cid:2)ura3DtrmB background, these genes were slowly
NucleicAcidsResearch,2013,Vol.41,No.18 8551
damped oscillations before reaching a new steady state
level (Figure 5A). To quantitatively assess the role
played by TrmB in the gene expression profile of ppsA,
we calculated a derived mRNA synthesis rate and
compared it with promoter occupancy. We reasoned
that a derived mRNA synthesis rate is representative of
actualpromoteractivityanddeconvolvestheeffectofdif-
fering mRNA degradation rates in different genes.
Derivation of an mRNA synthesis rate from mRNA
data requires the degradation constant (K ) for the
deg
gene of interest. For this, we used genome-wide mRNA
D
degradation rates that have been experimentally o
w
determined previously in H. salinarum (29). The derived n
lo
mRNA synthesis rate of ppsA inferred from the gene ex- a
d
pression data (Figure 5B) correlated strongly with the ed
ChIP-qPCR measurement of actual enrichment (Pearson fro
m
correlation=0.94). Taken together, these calculations h
suggest that TrmB is the primary regulator of ppsA ttp
s
mRNA synthesis. ://a
c
a
d
Prediction of gene expression based on TrmB-promoter e
Figure 3. NanoString gene expression measurements are tightly m
enrichment ic
correlated with qPCR measurements. mRNA levels for ppsA (black .o
points) and pykA (gray points) were measured by RT-qPCR and TrmB dynamics at the ppsA promoter explained the up
NanoString. Both RT-qPCR and NanoString are shown plotted on a temporal gene expression pattern observed for ppsA .co
logarithmic axis. Error bars represent standard error from two repli- m
cates of NanoString and at least two replicates of RT-qPCR. mRNA (Figure 5). This led us to ask to what extent /n
a
TrmB-promoter binding dynamics alone could explain r/a
and constantly upregulated following the addition of the temporal expression patterns of other genes in the rtic
TrmB regulon. To address this question, we modeled le
glucose without any rapid changes in expression level expression across the 100 genes in the NanoString data -ab
(SupplementaryFigureS1).Therewasnosignificantclus- s
tkenroincgkoouft (tThaesbelee1x)p.rTesosgioetnhepr,rothfielesse dinatathseug(cid:2)geusrtat3hDattrmcoB- spertomaostear fouvnecrtiotinmeo.fATsrmTBrmBenrbicinhdmsenat caits-rtehgeulaptposrAy tract/4
motif that is conserved throughout its regulon (TACT- 1
balamin biosynthesis is predominantly TrmB regulated /1
N -GAGTA) and its biochemical mechanism involves 8
but that other factors may be involved. (7–8) /8
reduced affinity for DNA when sugar is present in the 5
dyInnamsicumgmenaerye,xptrheessiosnurppraitstienrgnsdsiuvgegresisttys tohaft oTbrsmerBveids growth medium, we expect that TrmB binding will be 46/1
qualitatively similar at all regulated promoters. This is 03
required for temporal coordination of gene expression 5
supported by previous genome-wide ChIP-chip studies 5
across metabolism in response to glucose. These patterns 4
implicateotherunknownregulatoryfactorsinthisprocess. showing that TrmB bound its regulon only in the 2 b
The expression of genes coding for enzymes at the core of aabsssuemnceed offorglcuocmospeutaotriongalylcemroolde(l1in7g). pWurepotsheesretfhoaret y gue
central metabolism appears to be predominantly TrmB s
regulated, while branching cofactor pathways seem to TrmBdynamicsattheppsApromoterwererepresentative t on
require additional, as of yet unidentified regulators. of those at other promoters in the genome. ppsA 16
promoter-binding dynamics (Figure 5) were used as an N
input to an ODE model of gene expression with two par- ov
TrmB promoter occupancy can explain gene expression e
ameters: a basal synthesis rate (K ) and a scaling term m
dynamics for ppsA basal b
(K )(Methods). For 62 genes, this approach adequately e
To determine the specific contribution of TrmB-promoter expeflfained the dynamics and predicted the level of gene r 20
1
binding to gene expression dynamics, we assayed TrmB expression in the (cid:2)ura3DtrmB mutant strain (Figure 6A 8
enrichment over time in response to glucose at the ppsA and B and Supplementary Figure S1). For other genes,
promoter using chromatin immunoprecipitation followed this method explained certain aspects of the dynamics,
byqPCR(ChIP-qPCR,‘MaterialsandMethods’section). but not all. For yet other genes, the model was unable
Binding was assayed from 60min before glucose addition to explain the changes in gene expression (Figure 6C
through180minafteraddition.Thesametimepointsused and D). The modeling and prediction for each of these
for gene expression measurement were sampled with five groups is described in turn later in the text.
additional time points for increased resolution. During To characterize both the role and the temporal
mid-logarithmic growth on amino acids as a sole source dynamics of additional regulators working with TrmB
of carbon and energy, TrmB showed significant binding over time in response to glucose, we calculated the differ-
enrichment(Figure5A).Ontheadditionofglucose,TrmB ence in gene synthesis rate between the model prediction
dissociated from the ppsA promoter within 2min, and and the observed data. We then clustered the synthesis
binding dynamics exhibited slight but reproducible rate residuals (actual minus predicted) in the (cid:2)ura3
8552 NucleicAcidsResearch,2013,Vol.41,No.18
C
A D
D
o
w
n
B E lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
ic
.o
u
p
.c
o
m
/n
a
r/a
F rticle
-a
b
s
tra
c
t/4
1
/1
8
/8
5
4
G H 6/1
0
3
5
5
4
2
b
y
g
u
e
s
t o
n
1
6
N
o
v
e
m
b
e
r 2
0
1
8
Figure 4. TrmB is required for metabolic modularity in H. salinarum metabolism. Figure depicts k-means clustering and using eight clusters and
Pearsoncorrelationasascoringmetricsonthe(cid:2)ura3(C)and(cid:2)ura3DtrmB(F)data.Onlyclustersenriched(P<0.01)forCOGbiologicalfunction
categoriesareshownforeachstrain.Fourclustersareenrichedingenesinvolvedinspecificpathwaysin(cid:2)ura3(A,B,DandE),whereasonlytwo
clusters are thus enriched in (cid:2)ura3DtrmB (G and H). Metabolic maps (C and F) depict enzyme-coding genes as colored squares, with colors
corresponding to cluster graphs. Each cluster graph depicts gene expression data for individual genes (thinner lines) and the mean expression
profile for the cluster (thicker lines). The number of genes in every cluster that shows significant (P<0.05) 1.5-fold over expression by the one-
sided Welch t-test is listed in the lower right-hand corner of each expression graph. Metabolite abbreviations used in (C and F) are listed in
Supplementary Table S7. In all, 30 of the 100 genes analyzed by NanoString are not represented on the metabolic map for the sake of clarity.
NucleicAcidsResearch,2013,Vol.41,No.18 8553
Table 1. P-values of enrichment of COG terms in k-means clusters of dynamic expression profiles of the (cid:2)ura3 and (cid:2)ura3DtrmB strains after
glucose stimulus
COG term Description Dura3 Dura3DtrmB
Carbohydrate transport and metabolism Glycolysis/gluconeogenesis 1.08E-03 NS
Amino acid transport and metabolism 2.95E-03 7.32E-04
Nucleotide transport and metabolism Purine biosynthesis 1.35E-04 1.23E-09
Coenzyme transport and metabolism Cobalamin biosynthesis 4.47E-03 NS
NS=not significant.
D
o
A functions (Table 2). In two of the five clusters (Cluster 4 w
n
and5,Figure 7D), containing 62of100 genes,no pattern lo
a
was found in the residuals. This suggests that these genes d
e
eexxhpirbeistseidoneipthreedrilcittitolenschfaronmgetohveerTrthmeBt-iemnericchomuresnetomrothdaetl d from
were accurate (Figure 7D). As expected, all three of the h
no-change control genes were members of these two ttps
clusters. Further, carbohydrate metabolism genes were ://a
c
enriched in one of these clusters: they are predominantly a
d
TrmB regulated and were therefore accurately predicted em
(Figures 5 and 6). ic
.o
The three remaining clusters (Clusters 1, 2, 3; Figure 7 u
p
A–C) showed significant dynamic patterns in their model .c
o
residuals between 0 and 180min after perturbation with m
/n
glucose. This indicates that the model prediction deviated a
spioginnitfisc.aFnutlrythferro,mthitshdeevoibasteiorvnesdugdgaetsatedovtehratceortthaeinr rteigmue- r/article
lators besides TrmB may be involved. -a
b
B s
tra
Inferring regulatory network topology from dynamic gene c
expression output t/41
/1
To identify how such regulators may be involved, we 8/8
reversed biological circuit design principles (2) to infer 54
6
GRN structure from dynamical output. Patterns of regu- /1
0
lation that could not be predicted from TrmB-promoter- 3
5
binding data were analyzed using logic approximations. 54
2
We applied this to Clusters 1 and 3, which were b
y
enrichedforfunctionsinpurineandvitaminB-12biosyn- g
u
thesis, respectively (Figure 7 and Table 2). e
s
The (cid:2)ura3 residuals in Cluster 3 (enriched for genes t o
n
encoding cobalamin biosynthesis proteins) exhibited a 1
6
peak centered at 20–45min (Figure 7C and Table 2). N
Consistent with the gene expression data, the spike was ov
e
absent in the (cid:2)ura3DtrmB strain. The global time delay m
b
from changes in transcript to changes in protein level fol- e
lowingperturbationinH.salinarumhasbeenestimatedat r 20
1
(cid:2)16min based on parallel mRNA and proteomics time 8
course data (13). The similarity between global transcrip-
Figure 5. (A) TrmB-binding enrichment at the ppsA promoter is tion-translation time lag and the lag between glucose
correlated with (B) predicted mRNA synthesis rate inferred from the additionandthepeakobservedincobalaminbiosynthesis
gene expression data (‘Materials and Methods’ section) at the TrmB-
gene residuals at 20–45min suggest that TrmB and a
dependentppsApromoteraftertheadditionof5%glucose.Errorbars
represent the standard error from the average of biological replicate second TrmB-dependent regulator may be involved in a
experiments. See also Supplementary Figure S3. FFL to regulate cobalamin biosynthesis. As we observed
an immediate decrease in gene expression after TrmB
dissociated from the DNA, we reasoned that either
strainforeachgeneovertheglucoseresponsetimecourse TrmB or the FFL regulator could activate cobalamin
using k-means. This approach classified the accuracy of gene expression. This represents an OR logic gate in the
modelfitstothedataintofivedifferentgroups(Figure7). FFL.Totestwhetherourdatawereconsistentwithsucha
All five clusters were highly enriched for specific COG feed-forward regulatory topology, we simulated the
8554 NucleicAcidsResearch,2013,Vol.41,No.18
A gpm B pykA
0 0
00 Observed Δura3 00 Observed Δura3
2 1
Observed Δura3ΔtrmB Observed Δura3ΔtrmB
PredictedΔura3 PredictedΔura3
00 PredictedΔura3ΔtrmB 800 PredictedΔura3ΔtrmB
5
1
0
0
6
0
0
0
1
0
0
4 D
o
w
nts) 500 00 nloa
u 2 de
o d
c fro
S m
0 0
N h
d ttp
ze −200 0 100 200 300 −200 0 100 200 300 s://a
aliC cbiG D purD ca
norm 15000 OObbsseerrvveedd ΔΔuurraa33ΔtrmB 25000 OObbsseerrvveedd ΔΔuurraa33ΔtrmB demic.o
NA ( PPrreeddiicctteeddΔΔuurraa33ΔtrmB 000 PPrreeddiicctteeddΔΔuurraa33ΔtrmB up.co
mR 00 20 m/na
100 5000 r/artic
1 le
-a
0 b
0 s
0 00 tra
0 1 c
50 t/4
1
5000 /18/8
5
4
6
/1
0 0 03
5
5
4
−200 0 100 200 300 −200 0 100 200 300 2
b
y
Time after 5% glucose addition (min) g
u
e
s
Figure 6. ODEmodelfittoNanoStringgeneexpressiondatainresponsetoglucoseispredictiveofgeneexpressionintheDura3(parentstrain)and t o
(cid:2)ura3DtrmB mutant background for some but not all genes. (A) Phosphoglycerate mutase (gpm). (B) Pyruvate kinase (pykA). (C) Cobalamin n
biosynthesisgene(cbiG).(D)5-phosphoribosylglycinamide(GAR)synthetase(purD,purinebiosynthesispathway).BlacklinesrepresentmRNAlevel 16
in the Dura3 strain. Red lines represent mRNA level in the DtrmB strain. The gray line and the dotted orange line show model fits for Dura3 and N
(cid:2)ura3DtrmB,respectively. Errorbarsdepictstandarderrorfromtheaverageoftwobiologicalreplicatesofthegeneexpressiondata. Modelfitsto ov
e
data for the remaining 96 genes are exhibited in Supplementary Figure S1. m
b
e
r 2
0
proposed network using logic approximations (30) dynamic profiles of these genes make it clear that TrmB 18
(Figure 8A). We then fitted the parameters of the simula- is acting as an activator of cob/cbi genes, and that the
tion to our cobalamin biosynthesis (cob/cbi) gene expres- second regulator is acting as an activator as well.
sion data over the first 60min following glucose addition In contrast to the lagged residuals of the cobalamin
(‘Materials and Methods’ section). The fitted degradation biosynthesis cluster, the residuals of Cluster 2, enriched
rateoftheFFLeffectorwas0.034min(cid:4)1.Thisindicatesa in arginine and serine metabolism genes, exhibited an im-
response time of (cid:2)20min, which is consistent with the mediately changing dynamic residual profile in the (cid:2)ura3
16min lag between transcription and translation (13), parent strain (Figure 7B and Table 2). Model residuals
andsupportsthetranscription-translationFFLhypothesis increase immediately on glucose addition in the (cid:2)ura3
(Figure 8B). From this model, the degradation constant parent strain. This suggests that TrmB is not a primary
forthecob/cbimRNAwasestimatedat(cid:2)5.5min,whichis regulator of these genes. In contrast to the previous
similar to the average empirically determined half-life of cluster, because of the diversity of amino acid metabolic
the cob/cbi mRNA [8.0min; (29)]. Furthermore, the pathways and the likelihood of multiple independent
NucleicAcidsResearch,2013,Vol.41,No.18 8555
A Cluster 1 B Cluster 2
Δura3 Δura3
4 Δura3ΔtrmB 4 Δura3ΔtrmB
0 0
0. 0.
2 2
0 0
0. 0.
)
s 0 0
si 0 0
e 0. 0.
h
t
n
sy 0.02 0.02 Dow
d − − n
e lo
ct 4 4 ad
redi −0.0 −0.0 ed fro
p m
- h
al −200 0 100 200 300 −200 0 100 200 300 ttp
Actu C Cluster 3 D Clusters 4&5 s://aca
( d
del fit 0.04 ΔΔuurraa33ΔtrmB 0.04 ΔΔuurraa33ΔtrmB emic.ou
o p
m .c
o
of 0.02 0.02 m/n
a
ss r/a
odne 0.00 0.00 rticle-a
o b
G −0.02 −0.02 stract/4
1
−0.04 −0.04 /18/854
6
/1
0
3
−200 0 100 200 300 −200 0 100 200 300 5
5
4
Time after 5% glucose addition (min) 2
b
y
Figure 7. Thedifferencebetweenmodel-predictedandactualgenesynthesisratesin(cid:2)ura3(black)and(cid:2)ura3DtrmB(gray)inclustersovertimeinresponse g
u
toglucose.(A)Cluster1,enrichedingenescodingforenzymesinvolvedinpurinebiosynthesis;(B)Cluster2,enrichedingenescodingforenzymesinvolved e
s
inaminoacidmetabolism;(C)Cluster3,enrichedingenescodingforenzymesinvolvedincobalaminbiosynthesis;(D)Clusters4and5,containinggenes t o
whoseexpression iswell-fitted byour model. The errorbarsrepresent thestandard errorfrom the averageof themodel fit residualsin each cluster. n 1
6
N
o
v
the residuals were similar in both (cid:2)ura3 and e
Table 2. P-values of enrichment of COG terms in k-means clusters m
(cid:2)ura3DtrmB. As the response time through a TF would b
of (cid:2)ura3 model residuals e
beontheorderofthehalf-lifeofthisfactor[(cid:2)10min;(2)], r 2
0
Cluster number COG term P-value thesedataimplicateasecondaryregulatoralreadypresent 1
8
in the cytoplasm. Although TrmB may ultimately control
Cluster 1 Nucleotide transport and metabolism 6.30E-03
the concentration of such a secondary regulator
Cluster 2 Amino acid transport and metabolism 3.19E-02
Cluster 3 Coenzyme transport and metabolism 1.39E-06 (Supplementary Figure S2), initial dynamics are likely
Cluster 4 Carbohydrate transport and metabolism 7.24E-03 TrmB-independent. As theinputsignals tothis secondary
Cluster 5 Energy production and conversion 4.73E-02 regulatorareunknown,thisnetworkisnotamenabletoa
logic simulation. Nevertheless, these data support the hy-
pothesisthatasecondary regulatormayworkthrough an
regulators, further understanding of network topology of independentpathwaytoimpingeonthepurinebiosynthe-
this cluster will require additional studies. sis genes coordinately with TrmB (Figure 6D).
The residuals of Cluster 1, enriched for functions in In summary, predictions from ODE and logic models
purine biosynthesis (Figure 7A and Table 2), increased suggestthatTrmBisresponsible,invariouscapacities,for
immediately following glucose addition. Furthermore, the temporal regulation of metabolism as well as the