extfieldsofstudy: Biology Medicine

text:
Transcription is regulated by NusA:NusG interaction

NusA and NusG are major regulators of bacterial transcription elongation, which act either in concert or antagonistically. Both bind to RNA polymerase (RNAP), regulating pausing as well as intrinsic and Rho-dependent termination. Here, we demonstrate by nuclear magnetic resonance spectroscopy that the Escherichia coli NusG amino-terminal domain forms a complex with the acidic repeat domain 2 (AR2) of NusA. The interaction surface of either transcription factor overlaps with the respective binding site for RNAP. We show that NusA-AR2 is able to remove NusG from RNAP. Our in vivo and in vitro results suggest that interaction between NusA and NusG could play various regulatory roles during transcription, including recruitment of NusG to RNAP, resynchronization of transcription:translation coupling, and modulation of termination efficiency.

INTRODUCTION
Transcription, the first step in gene expression, is highly regulated by a multitude of transcription factors. The core transcription machinery is RNA polymerase (RNAP), which consists of five subunits in bacteria (2 x ␣, ␤, ␤', and ) (1). RNAP initiates RNA synthesis at a promoter (initiation), extends the nascent RNA (elongation), and releases the RNA at a terminator (termination) (2). Among Escherichia coli (E. coli) transcription factors, N-Utilization Substances (Nus) A and G are the only ones known to affect both the speed of RNA chain elongation and termination (3). NusG (Spt5 in archaea and eukaryotes), the only universally conserved transcription factor, is composed of an N-terminal domain (NTD) flexibly connected to a Cterminal domain (CTD) (Supplementary Figure S1A) (4,5). NusG-NTD interacts with the ␤' clamp helices (␤'CH) and the ␤ gate loop (␤GL) of RNAP to increase RNAP processivity (4,6,7). NusG-CTD is target of at least two cellular partners, termination factor Rho and antitermination factor NusE, which is identical to ribosomal protein S10 (8,9). NusA comprises an NTD that binds to the ␤ flap tip helix of RNAP at the RNA exit channel, three RNA binding motifs (S1, KH1, KH2) that together form the SKK domain, and, in E. coli and other ␥ -proteobacteria, two additional C-terminal acidic repeat domains (AR1, AR2; Supplementary Figure S1B) (10)(11)(12)(13)(14). NusA-AR1 interacts with N protein of phage , but is not essential for N-mediated suppression of transcription termination (antitermination) (15)(16)(17). NusA-AR2 can bind either to the CTD of the RNAP ␣-subunit (␣CTD) or to the NusA-SKK. NusA-AR2 attached to NusA-SKK autoinhibits NusA activity by preventing RNA binding (15,18).
NusA and NusG differentially alter the properties of the transcription elongation complex (TEC) via direct and independent interactions (3). NusG increases TEC processivity whereas NusA slows RNAP by either increasing pause times or by introducing new pause sites (19). The two factors have context-dependent effects on termination and act either in concert or as antagonists. On the one hand, NusG and NusA are proposed to support Rho cooperatively to suppress the toxic functions of foreign genes. On the other hand, both are part of antitermination complexes on ribosomal RNA and phage nascent transcripts (20)(21)(22)(23). Furthermore, NusA can enhance or decrease both Rho-dependent and intrinsic termination efficiency, depending on the specific terminator (reviewed in (24,25)). NusG, in contrast, enhances termination exclusively at Rho-dependent sites (26,27). Importantly, NusG serves as the physical linker between the RNAP and the ribosome by binding RNAP via NusG-NTD and S10 via NusG-CTD, thus coupling transcription and translation (8).
NusA and NusG bind to different sites on RNAP (6,7,10). Although these sites are in close proximity, a direct connection between the two factors has not been reported thus far. With nuclear magnetic resonance (NMR) spectroscopy we here demonstrate that NusA and NusG do specifically mutually interact and, supported by in vivo and in vitro data, we propose that this interaction may have key regulatory roles in diverse steps of transcription.

Cloning
nusA was cloned into the pTKK19 expression vector (28) via NdeI and BamHI restriction sites resulting in the recombinant plasmid pTKK19 nusA . The recombinant protein carried a deca-histidine tag followed by a PreScission cleavage site at its N-terminus.
nusA was expressed in E. coli BL21 (DE3) harboring pTKK19 nusA . Cells were grown in lysogeny broth (LB) medium containing 30 g/ml kanamycin at 37 • C. At an optical density at 600 nm (OD 600 ) of ∼0.7 expression was induced by 1 mM isopropylthiogalactoside (IPTG). Cells were harvested after 4 h (9,000 x g, 15 min, 4 • C), resuspended in buffer A (20 mM Tris(hydroxymethyl)aminomethane (Tris)/HCl, pH 7.9, 500 mM NaCl, 5 mM imidazole, 1 mM ␤-mercaptoethanol) and disrupted by a microfluidizer (Microfluidics, Newton, MA, USA). After centrifugation (12,000 x g, 30 min, 4 • C) the crude extract was applied to a 5 ml HisTrap HP column (GE Healthcare, Munich, Germany) and eluted using a step gradient from 5 mM to 1 M imidazole in buffer A. Fractions containing the His 10 -NusA fusion protein were combined and the protein was cleaved by PreScission protease during dialysis against buffer B (20 mM Tris/HCl, pH 8, 1 mM ␤-mercaptoethanol) at 4 • C overnight. The protein solution was applied to a 5 ml GSTrap FF column (GE Healthcare, Munich, Germany) and the flow-through subsequently to a 5 ml QXL column (GE Healthcare, Munich, Germany). NusA was eluted using a step gradient from 0 to 1 M NaCl in buffer B. Fractions containing pure NusA were combined, dialyzed against 5 l 20 mM Tris/HCl, pH 7.5, 50 mM NaCl, 1 mM dithiothreitol (DTT), concentrated using ultrafiltration units (Viva Science, molecular weight cutoff (MWCO): 10 kDa), shock frozen in liquid nitrogen, and stored at −80 • C.
Production and purification of RNAP for NMR experiments was based on Ref. (34). In brief, E. coli BL21(DE3) (Novagen, Madison, WI, USA) harboring a plasmid containing rpoA, rpoB, rpoC and rpoZ in one operon under control of T7 promoter was used for gene expression. The ␤' subunit was produced as fusion protein carrying a His 6 tag at its C-terminus. 2 l of LB medium in a 5 l flask supplemented with 100 g/ml ampicillin were inoculated with an overnight culture to an OD 600 of 0.02 and incubated at 37 • C and 150 rpm. Having reached an OD 600 of 0.2 the temperature was decreased to 20 • C. After 2 h, IPTG was added to a final concentration of 0.5 mM for induction, and the culture was incubated overnight. Cells were then harvested by centrifugation for 15 min at 4 • C and 9,000 x g, resuspended in buffer C (50 mM Tris/HCl, pH 6.9, 0.5 M NaCl, 5 % (v/v) glycerol, 1 mM DTT) and lysed with a microfluidizer (Microfluidics, Newton, MA, USA). RNAP was purified by nickel affinity chromatography with 2×5 ml Ni 2+ -nitrilotriacetic acid (NTA) Superflow columns (QIAGEN) and eluted by a constant gradient from 0 to 1 M imidazole in buffer C. Peak fractions containing RNAP were dialyzed against buffer D (50 mM Tris/HCl, pH 6.9, 0.5 mM ethylenediaminetetraacetic acid (EDTA), 5 % (v/v) glycerol, 1 mM DTT) overnight at 4 • C and then applied to a 5 ml Heparin HP column (GE Healthcare), followed by elution with a constant NaCl gradient from 0 to 1 M NaCl in buffer D. Fractions containing RNAP were pooled, dialyzed against buffer B overnight at 4 • C, applied to a 25 ml Q-Sepharose FF column, and eluted with a constant gradient from 0 to 1 M NaCl in buffer B. RNAP containing fractions were concentrated using ultrafiltration units (Viva Science, MWCO: 10 kDa), shock frozen in liquid nitrogen, and stored at −80 • C.
Production and purification of RNAP for in vitro transcription assays was carried out as described (35), with minor modifications.

Isotopic labeling of proteins
Proteins were uniformly labeled with 15 N or 15 N, 13 C by growing E. coli in M9 minimal medium (36,37) supplemented with ( 15 NH 4 ) 2 SO 4 (Campro Scientific, Berlin, Germany) or ( 15 NH 4 ) 2 SO 4 and 13 C-D-glucose (Cambridge Isotope laboratories, Inc., Andover, MA, USA) as the only nitrogen and carbon source, respectively. Expression and purification procedures were identical to those used for proteins produced in LB medium.

Pull-down assay
The pull-down assay was performed with a 1 ml HisTrap HP column (GE Healthcare) equilibrated with buffer E (10 mM potassium phosphate, pH 6.9, 50 mM NaCl). The proteins were dialyzed against buffer E overnight at 4 • C before application to the column. The application volume was always 1 ml. After extensive washing with buffer E, elution was carried out with 100 or 400 mM imidazole in buffer E, respectively, and resulted in a mixture of His 10 -NusA-AR2 and NusG-NTD (400 M and 200 M, molar ratio: 2:1) as assayed, with His 10 -NusA-AR2 (200 M) alone and NusG-NTD (400 M) alone as controls.

In vitro transcription assay
RNA and DNA oligonucleotides were commercially synthesized by Integrated DNA Technologies with sequences derived from the T7A1 promoter sequence. Assembly of the TEC and the in vitro transcription were carried out as described (38). Briefly, a 65mer template DNA strand was hybridized to an 11mer RNA labeled with 32 P at the 5 end. 50 pmol RNAP in transcription buffer (TB, 20 mM Tris/HCl, pH 7.9, 5 mM MgCl 2 , 40 mM KCl, 2 mM ␤mercaptoethanol) were mixed with equimolar concentrations of the DNA:RNA hybrid, followed by addition of the nontemplate DNA strand. The assembled TECs were purified by affinity chromatography using Ni 2+ -nitrilotriacetic acid agarose (QIAGEN) and subsequent membrane filtration with Ultrafree C 0.65 m PVDF centrifugal filters (Millipore). 2.5 M TEC were incubated with 50 M of one or both transcription factors or TB for 10 min at 25 • C. When two proteins were tested they were added simultaneously, as previous tests indicated that the order of addition had no effect. Transcription was initiated by addition of 1 mM NTPs and stopped after 60 seconds by addition of an equal volume of 2x loading buffer (10M urea, 50 mM EDTA, pH 7.9, 0.05 % (w/v) bromophenol blue, and 0.05 % (w/v) xylene cyanol). RNA products were resolved on a 23 % denaturing polyacrylamide gel containing 7 M urea. Gels were exposed to phosphor screens and scanned by Typhoon Phosphorimager (GE Healthcare Life Sciences).

NMR experiments
All NMR spectra were recorded at 298 K on Bruker Avance 700 MHz and Avance 800 MHz spectrometers with cryogenically cooled triple-resonance probes equipped with pulsed field-gradient capabilities. Processing of NMR data was carried out using in-house routines and visualized by NMRView (39). For all NMR experiments the proteins were in 10 mM potassium phosphate buffer, pH 6.4, 50 mM NaCl, at 298 K. The initial sample volume was 550 l if not stated otherwise. Backbone assignments of NusG-NTD, NusG-CTD, NusA-AR2, and ␣CTD were taken from previous studies (4,18,32).
Dissociation constants (K D ) were calculated from [ 1 H, 15 N]-HSQC titrations by analyzing the chemical shift changes and fitting a two-state model as in Equation (2) to the chemical shift change of amide protons showing fast exchange on the chemical shift timescale.
where ν is the normalized resonance frequency difference in Hz, ν End the normalized resonance frequency difference between free and fully bound protein in Hz, K D the dissociation constant, r the protein:labeled protein ratio and [P] 0 the total concentration of 15 N-labeled protein. K D and ν End were used as fitting parameters. The reduction of [P] 0 due to dilution was accounted for during fitting.
For the displacement experiment of 15 N-NusA-AR2 from ␣CTD by NusG-NTD separate samples were prepared for 15 N-NusA-AR2 (100 M, 500 l) and 15 N-NusA-AR2: ␣CTD (100 M each, 500 l). NusG-NTD was then added to the latter sample from a 287 M stock solution.
For the quantitative analysis of signal intensities in the displacement experiments signal intensities were normalized by the number of scans, the concentration, and the length of the 90 • proton pulse.

Strains and ␤-galactosidase assays
Strains were derivatives of MDS42, which lacks prophages and insertion elements (41) containing fusions (42). lacZ is expressed from the fusions λcI857 -pR -cro(ΔRBS) -nutR -tR1 -cII::lacZ or λcI857 -pR -cro 27 -nutR -tR1 -cII::lacZ. The TAAGGAGGTTGT to TaccctccTTGT substitution in the cro ribosome binding site (RBS), blocks cro translation in the former strain (cro(ΔRBS)). cro 27 has an RBS, but Cro 27 is non functional due to an amino acid exchange (R27Q). The creation of strains 10323 and 10881 was described previously (18). The rhoE134K mutation was introduced by phage P1 transduction, resulting in a Rho variant which is non functional at tR1. The strains carrying nusA variants were constructed by recombineering.
Cells were assayed for ␤-galactosidase activity (Miller units) after overnight growth at 37 • C. ␤-galactosidase activity of cells with defective Rho (rho15) was set to 100 % (strains 11633 and 11634), since termination was completely abolished at λtR1. Assays were performed independently four times and resulting activities were averaged.

Programmes
Graphical representations of protein structures were created with PyMOL (43). Sequence alignments were done with Clustal omega (44).

NusG interacts specifically with NusA
First, we probed a possible NusA:NusG interaction by NMR spectroscopy with full-length proteins. Addition of NusA to 15 N-labeled NusG to equimolar concentration resulted in a strong decrease of NusG-NTD signals in the [ 1 H, 15 N]-HSQC spectrum, whereas NusG-CTD signals were weakened only marginally ( Figure 1A). The high transversal relaxation rate of the 54.9 kDa NusA strongly affects magnetization transfer efficiency upon binding, which leads to line broadening and ultimately to a decrease of signal intensity. Thus, the observed loss of NusG-NTD signals suggests direct NusA:NusG-NTD interaction. Specific NusA:NusG-NTD complex formation was confirmed by addition of NusA to either isolated 15 N-NusG-NTD or 15 N-NusG-CTD as signal loss was only observed for 15 N-NusG-NTD ( Figure 1B and C). [ 1 H, 15 N]-HSQC spectra of a 2-fold molar excess of 15 N-NusG-NTD in the presence of isolated NusA-NTD, NusA-SKK, NusA-AR1, or NusA-AR2 clearly showed that, of these, only NusA-AR2 interacted directly with NusG-NTD (Figure 2).
To corroborate the results of the NMR experiments, a pull-down assay was performed with decahistidine-tagged NusA-AR2 (His 10 -NusA-AR2) and untagged NusG-NTD (Supplementary Figure S2). NusG-NTD eluted from the nickel column together with His 10 -NusA-AR2, thus confirming the direct interaction between the two domains.  Figure S3). By mapping the normalized chemical shift changes on the three-dimensional structures of NusG-NTD and NusA-AR2 the interaction interface can be located in the C-terminal part of helix ␣5 of NusA-AR2, markedly involving W490 and F491 ( Figure  3B and C). Although NusA-AR1 and NusA-AR2 share high sequence identity (31.5%) and have virtually identical three-dimensional structures with a root mean square deviation of main chain atoms of 1.2Å (Supplementary Figure  S4), NusG-NTD recognizes NusA-AR2 exclusively ( Figure  2). This selectivity can probably be attributed to the presence of W490 and F491 in NusA-AR2, since Leu and Ala are found at corresponding positions in NusA-AR1 (Supplementary Figure S4). Different binding specificities for NusA-AR1 and NusA-AR2 to protein N and the ␣CTD of RNAP have been noted earlier (13,16,18). NusG-NTD signals from residues in the C-terminal helix ␣3' (aa 104-117) and in the elongated loop region between ␣1' and ␤1' (aa 48-68) are mainly affected by the NusG-NTD:NusA-AR2 interaction ( Figure 3B and C). Based on these binding surfaces a docking model without conformational rearrangements of the complex was generated ( Figure 3C). Remarkably, the NusA-AR2 binding site on NusG-NTD is also involved in the NusG-NTD:RNAP ␤'CH interaction (Supplementary Figure S5) (7,45). Furthermore, NusA-AR2 residues responsible for NusG-NTD binding are necessary for the interaction with ␣CTD ( Supplementary Figure S6) (18). Thus, formation of the NusG-NTD:RNAP and NusG-NTD:NusA-AR2 complex is mutually exclusive, as is formation of the NusG-NTD:NusA-AR2 and NusA-AR2:RNAP complex.

NusG-NTD:RNAP versus NusG-NTD:NusA-AR2 versus NusA-AR2:RNAP
We asked if NusG-NTD:NusA-AR2 interaction can take place in the presence of RNAP. Therefore a [ 1 H, 15 N]-HSQC displacement experiment with isolated ␣CTD was performed ( Figure 4A,B). Addition of NusA-AR2 to 15 N-␣CTD to equimolar concentration induced the chemical shift perturbations of 15 N-␣CTD signals observed earlier for this interaction (18). Stepwise addition of NusG-NTD to a final three-fold molar excess of NusG-NTD partially reversed these shifts, indicating that NusG-NTD can displace NusA-AR2 from the ␣CTD. Although the displacement was incomplete due to the lower K D of NusA-AR2:␣CTD (K D < 5 M) (18) versus NusA-AR2:NusG-NTD (K D : 24 M), it confirms that the RNAP/NusA-AR2 binding sites on NusG-NTD overlap as do the ␣CTD/NusG-NTD interaction interfaces on NusA-AR2. This finding was counter-checked by the displacement of 15 N-NusA-AR2 from ␣CTD by NusG-NTD (Supplementary Figure  S7). Addition of ␣CTD to 15 N-NusA-AR2 to equimolar concentration resulted in chemical shift changes of signals from 15 N-NusA-AR2 typical for 15 N-NusA-AR2:␣CTD complex formation (18). Subsequent addition of NusG-NTD caused the 15 N-NusA-AR2 signals to shift towards the resonance positions of the 15 N-NusA-AR2:NusG-NTD complex. Again, the displacement was incomplete owing to the difference in the affinities of NusA-AR2 to ␣CTD and NusG-NTD.
We extended these studies to full-length proteins stepby-step. In an initial experiment we added RNAP to 15 N-NusG, which led to a loss of almost all 15 N-NusG signals in the [ 1 H, 15 N]-HSQC spectrum owing to the dramatic increase of the NusG rotational correlation time upon formation of the NusG:RNAP complex (Supplementary Figure S8A). Although NusG-NTD is only flexibly linked to NusG-CTD and NusG is supposed to interact with RNAP via NusG-NTD (4), NusG-CTD signals were not observable in the 15 N-NusG:RNAP complex. Thus, either NusG-CTD is sterically hindered in the complex so that it cannot move freely, or NusG-CTD interacts directly with RNAP. To exclude such direct NusG-CTD:RNAP interac-  tion, we monitored isolated 15 N-NusG-CTD in the presence of RNAP in a separate experiment and found no changes in the [ 1 H, 15 N]-HSQC spectra as compared to 15 N-NusG-CTD alone (Supplementary Figure S8B). Consequently, the loss of all 15 N-NusG signals upon RNAP addition solely originates from formation of the NusG-NTD:RNAP complex.
To probe the NusA:NusG interaction in the presence of full-length RNAP we added NusA to 15 N-NusG, leaving only NusG-CTD signals visible ( Figure 4C and Supplementary Figure S8C). On addition of RNAP all signals decreased ( Figure 4C and Supplementary Figure S8C). Thus, either (i) NusA:NusG interaction is disrupted by RNAP and both NusA and NusG bind individually to RNAP, or (ii) NusA:NusG remains intact and interacts with RNAP via NusA-NTD, or (iii) both. To eliminate interference by NusA-NTD:RNAP interactions, we repeated the experiment using isolated NusA-AR2 instead of fulllength NusA ( Figure 4D and Supplementary Figure S8D). When NusA-AR2 was present, the [ 1 H, 15 N]-HSQC spectrum of 15 N-NusG showed chemical shift changes corresponding to NusG-NTD:NusA-AR2 complex formation (see also Figure 2D). On addition of RNAP the intensity of the NusG signals decreased, however, in contrast to the experiment with full-length NusA, both NusG-NTD and NusG-CTD signals remained visible with the chemical shifts of the NusG-NTD:NusA-AR2 complex. We conclude (i) that at least some of the NusG-NTD:NusA-AR2 complexes remain intact in the presence of RNAP and (ii) that these complexes cannot bind to RNAP in the absence of NusA-NTD. This confirms that NusA-AR2:NusG and NusG:RNAP formation are mutually exclusive. The decrease in signal intensity may be explained by dissociation of a certain portion of the NusG-NTD:NusA-AR2 complex so that NusG binds to RNAP, and NusA-AR2 either interacts with the ␣CTD of RNAP or remains free. Thus, with full length NusA, the NusG:NusA complex is stable and is connected to RNAP via NusA-NTD, although a fraction of NusG and NusA might interact with RNAP individually.

NusA-AR2 induces pausing and blocks NusG suppression of pausing in vitro
We then asked if the NusA-AR2:NusG-NTD interaction affects transcription in vitro. For this we utilized a nucleic acid scaffold to generate a transcription elongation complex (TEC) which carries an 11 nt 32 P-labeled RNA primer basepaired to template DNA and flanked by non-template DNA (TEC11A, for details see Materials and Methods). Transcription was initiated by the addition of the four NTPs and stopped after 60 seconds. The TEC paused at several intrinsic pause sites in the template in the absence of additional transcription factors, with pause 1 being the most prominent ( Figure 5; lane 2). Full-length NusG suppressed pause Full-length NusA enhanced pause 1, but did not, however, induce pausing at pause 3 (lane 8). This suggests that the NusA-AR2 domain in full-length NusA was still bound to the NusA-SKK domain, and was not free to interact with the initiating TEC. Unlike NusA-AR2, suppression of pause 1 by full-length NusG or NusG-NTD was not completely abrogated by full-length NusA, possibly because NusA-AR2 remains bound to the NusA-SKK domain during elongation, and is not available to interact with NusG-NTD (lanes [8][9][10]. The ability of isolated NusA-AR2 to pause the TEC at 12C (pause 3) was unexpected. It suggests that at least early in elongation, when RNA has not yet extruded from the exit channel, NusA-AR2 may still be bound to the SKK domain and may thus be unavailable to interact with ␣CTD. In addition, the ability of NusA-AR2 to induce a strong pause implies that NusA-AR2 might act as a regulatory element during elongation if dissociated from ␣CTD.

NusA-AR2 blocks Rho-dependent termination at tR1 in vivo
NusA suppresses termination at certain Rho-dependent sites, e.g. within tR1 (46). We propose that the NusA:NusG interaction contributes to this effect. To support this hypothesis we asked if a deletion of NusA-AR2 (NusA-AR2) affected termination at tR1 in vivo ( Figure 6). We performed ␤-galactosidase assays using two fusions to measure termination: λcI857 -pR -cro(ΔRBS) -nutR -tR1 -cII::lacZ and λcI857 -pR -cro 27 -nutR -tR1 -cII::lacZ. Termination at tR1 is indicated by low ␤-galactosidase activity. Strains carrying a mutation in the rho gene show no termination at tR1; ␤-galactosidase activity of these strains was thus set to 100%. The efficiency of termination at tR1 was 93% when cro translation was prevented by an RBS mutation (strain 10323), and reduced to 74 % when cro was translated (strain 11149). Ribosomes reduce the amount of free RNA upstream to tR1 that is available to Rho, and thus block a Rho-binding site (rut) in cro (47).
We found that termination efficiency was significantly impaired (50%) in the nusA-ΔAR2 mutant only when cro mRNA was translated (compare strain 10881 to strain 10323 and strain 11699 to strain 11149). Our results suggest that NusA-AR2 may compete with Rho binding near nutR, the only rut site available when the cro transcript is occluded by ribosomes. When cro is not translated, Rho can attach to the free cro rut site (48). Competition with Rho by NusA-AR2 can be explained by constitutive binding of NusA-AR2 to RNA via NusA-SKK. Consistent with the in vitro studies described above, this implies that full length NusA may still be, at least partially, in the autoinhibited state and unable to bind rut RNA at tR1. It also raises the possibility that NusG may activate RNA binding of full-length NusA by displacing NusA-AR2 from the NusA-SKK domain. Further experimentation will be needed to address these questions.

Possible regulatory functions of the NusG-NTD:NusA-AR2 interaction
Transcription factors NusA and NusG act independently to slow or accelerate transcription elongation, respectively. The two factors can also function synergistically in restoring termination by the partially defective RhoE134K mutant (49). Here, however, we demonstrate a direct interaction between NusA and NusG. Formation of a complex between NusG-NTD and NusA-AR2 may explain their various in vivo and in vitro combinatorial regulatory effects (Figure 7).
NusG:NusA interaction may be involved in recruiting NusG to the TEC. ChIP-chip data suggest that NusA and NusG associate with RNAP when the TEC has exited the promoter region, with NusG attaching after NusA does (50). The delay in NusG association may be due to a competition between 70 region 2 and NusG-NTD, since both bind to the RNAP ␤'CH (7,(51)(52)(53).
Recall that the 70 region 2 can remain bound to the TEC even when 70 region 4 has dissociated from the ␤ flap, allowing NusA-NTD binding (54,55). After promoter escape, NusA attaches to the ␣CTD via NusA-AR2, to the ␤ flap via NusA-NTD, and to nascent RNA via NusA-SKK (10,15,18,56). Thus, stable tethering of NusA to the TEC may not require continuous binding of NusA-AR2 to ␣CTD. NusA-AR2 could, therefore, bind NusG-NTD without disrupting the NusA:TEC complex. In this model, NusA serves as a long linker to increase the local concentration of NusG, facilitating NusG recruitment to the TEC and displacement of 70 region 2 from the ␤'CH at the transition from initiation to elongation. 2. NusG-NTD:NusA-AR2 interaction assists transcription:translation coupling. As a second function, we suggest that NusG:NusA interaction plays a role in coupling transcription and translation. NusG connects these two processes by physically linking RNAP and the leading ribosome via NusG-NTD:RNAP and NusG-CTD:S10 interactions (4,8).
The NusA:NusG interaction could serve to resynchronize transcription and translation by coordinating the movements of RNAP and the ribosome. If translation is slowed, transcription could likewise be slowed by the temporary removal of NusG-NTD from the TEC by NusA-AR2. Since the NusG:ribosome connection remains intact, transcription and translation can be kinetically resynchronized. Also, the initial coupling of transcription and translation may occur via the NusA:NusG linker. This would explain the apparent late association of NusG with the TEC, as detected by ChIP-chip experiments (50). 3. NusG-NTD:NusA-AR2 interaction regulates transcription termination. The NusA:NusG interaction could influence contextdependent intrinsic or Rho-dependent transcription termination. In the former, the TEC pauses at an intrinsic termination signal, enters an elemental pause state, and then isomerizes into the termination state where the termination hairpin is formed (57). NusA-AR2 can remove NusG-NTD from RNAP, resulting in loss of both NusA-AR2 and NusG-NTD contacts to the TEC ( Figure 4E). Loss of these interactions might partially destabilize the TEC, open the clamp around the nucleic acids, and facilitate intrinsic termination. In contrast, the NusA-AR2:NusG-NTD interaction might decrease Rho-dependent termination. Binding of NusA-AR2 to NusG-NTD would release NusA autoinhibition, enhancing binding of NusA-SKK to nascent RNA to block Rho recruitment.
The K D values for the various interactions suggest that scenario 1 is the most probable. This scenario is also con-sistent with a global survey of distribution of transcription factors (50). The relevance of NusG-NTD:NusA-AR2 interaction in detail will need further experimental clarification, but the finding that NusG interacts directly with NusA may explain the various effects of these transcription factors on elongation and termination reported here and earlier.


s2fieldsofstudy:
Biology
Agricultural And Food Sciences

text:
Sequence diversification in recessive alleles of two host factor genes suggests adaptive selection for bymovirus resistance in cultivated barley from East Asia

Key message Two distinct patterns of sequence diversity for the recessive alleles of two host factors HvPDIL5 - 1 and HvEIF4E indicated the adaptive selection for bymovirus resistance in cultivated barley from East Asia. Abstract Plant pathogens are constantly challenging plant fitness and driving resistance gene evolution in host species. Little is known about the evolution of sequence diversity in host recessive resistance genes that interact with plant viruses. Here, by combining previously published and newly generated targeted re-sequencing information, we systematically analyzed natural variation in a broad collection of wild (Hordeum spontaneum; Hs) and domesticated barleys (Hordeum vulgare; Hv) using the full-length coding sequence of the two host factor genes, HvPDIL5-1 and HvEIF4E, conferring recessive resistance to the agriculturally important Barley yellow mosaic virus (BaYMV) and Barley mild mosaic virus (BaMMV). Interestingly, two types of gene evolution conferred by sequence variation in domesticated barley, but not in wild barley were observed. Whereas resistance-conferring alleles of HvEIF4E exclusively contained non-synonymous amino acid substitutions (including in-frame sequence deletions and insertions), loss-of-function alleles were predominantly responsible for the HvPDIL5-1 conferred bymovirus resistance. A strong correlation between the geographic origin and the frequency of barley accessions carrying resistance-conferring alleles was evident for each of the two host factor genes, indicating adaptive selection for bymovirus resistance in cultivated barley from East Asia. Electronic supplementary material The online version of this article (doi:10.1007/s00122-016-2814-z) contains supplementary material, which is available to authorized users.

Introduction
Plant pathogens such as viruses, bacteria and fungi are constantly challenging crop plant fitness since infections may cause yield losses (Oerke 2006). Resistant genotypes have the selective advantage of producing (more) seeds, thus their offspring may predominate in populations after some generations under strong selection. Some pathogen isolates may counteract this via own rapid evolution through mutation-this triggers the 'arms race' or the antagonistic co-evolution between plants and their adapted pathogens (Jones and Dangl 2006;Obbard and Dudas 2014;Rausher 2001). In natural populations or wild species a higher level of genetic diversity could maintain natural resistances that prevent the epidemic outbreak of diseases and their severe consequences for the populations or species (Tanksley and McCouch 1997). By contrast, the monoculture in agricultural systems destabilizes the balance, and the densely grown and genetically uniform crop populations facilitate the success of individual newly evolved and virulent strains of a pathogen supporting the chance of pathogen transmission and epidemic outbreak of respective diseases (Stukenbrock and McDonald 2008). Therefore, understanding the patterns of sequence variation underlying the evolution of dominant resistance genes and recessive host factors in cultivated crops and their wild progenitors should contribute to reveal the evolutionary processes during domestication, improvement and local adaptation of crops and thus may provide insight into a more sustainable use of resistance resources in agriculture.
Unlike bacterial and fungal pathogens that maintain relatively large genomes, plant RNA viruses have small genomes and commonly encode only for about four to ten viral proteins. Such viruses cannot complete their entire proliferation cycle independently of their plant hosts. For instance, they are dependent on the host's cellular protein synthesis machinery and also often employ other hostencoded proteins to complete their lifecycle (Ahlquist et al. 2003;Robaglia and Caranta 2006). These host components are known as 'susceptibility' or 'host' factors. Mutations in the genes coding for host factors may result in virus resistance (Soosaar et al. 2005;Whitham and Wang 2004). These recessive resistance genes, therefore, function very differently compared to the dominant resistance genes such as the nucleotide binding leucine rich repeat (NB-LRR) protein encoded genes (Sacristan and Garcia-Arenal 2008), which recognize directly or indirectly the so-called effector proteins produced by the pathogens and then trigger an active defense response often culminating in plant cell death and halting the pathogen growth (Tiffin and Moeller 2006). Depending on the selection pressure imposed by the frequency and/or severity of the respective virus disease, sequence variation in the host genes conferring disease resistance may be quickly selected and maintained at higher frequency in populations under conditions of natural selection. Artificial selection during domestication, adaptive improvement and breeding could even accelerate this selection process. Rapid evolution of plant virus genomes, which occurs constantly at least in RNA viruses due to the lack of an efficient RNA sequence repair system (Domingo and Holland 1997), may result in emergence of new virus pathotypes being able to overcome the newly evolved plant defense mechanism, e.g. by cooperating with alternate host factors or by adapting to the modified host proteins (Kühne et al. 2003;Okada et al. 2003). In agricultural systems, infectious diseases reduce yields and can threaten the sustainability of crop production (Oerke 2006). Selection for resistance forms of host factors was probably unintentional during the early domestication and an adaptation to local environments (Kang et al. 2005). Understanding the natural variation of such host factors could reveal the mechanisms underlying the 'arms race' between cultivated plants and their viruses.
The barley yellow mosaic virus disease, wherever it occurs, is a severe, long-term and persistent biotic stress, threatening winter barley production in a majority of European and Eastern Asian countries (Kühne 2009). It was first detected in Japan in the 1940s (Ikata and Kawai 1940), and later in the 1950s observed in China and Korea (Lee et al. 1996;Zhou and Cao 1985), and since 1978 reported to occur in European countries (Kühne 2009 (Adams et al. 1988). Persistent resting spores of this microorganism provide protection for virus particles, which can remain infectious for many years (Huth 1991). There are currently no chemicals available for control of plant viruses, whereas use of chemical agents (i.e. soil fumigants) effective against the virus vector, P. graminis, is undesirable due to ecological implications Kühne 2009). After infection, yield losses of up to 50% have been observed in Europe (Plumb et al. 1986). In Asia, the disease can be even more severe where complete yield losses were recorded in China in the mid-1970s (Chen 2005;Kühne 2009). Breeding for genetic resistance to bymoviruses, therefore, remains the most viable and environmentally safe option for crop protection, and over time eighteen bymovirus resistance loci were genetically characterized in barley (Kai et al. 2012;Ordon et al. 2005). We previously isolated by positional cloning two recessive resistance genes rym1/11 and rym4/5, corresponding to two different host factor genes Protein Disulfide Isomerase Like 5-1 (HvPDIL5-1) and Eukaryotic Translation Initiation Factor 4E (HvEIF4E), respectively (Kanyuka et al. 2005;Stein et al. 2005;Yang et al. 2014a, b). The resistance genes/alleles rym1 and rym5 were originally found in a Chinese landrace "Mokusseko 3" and shown to be inherited independently (Konishi et al. 1997). The HvEIF4E protein putatively functions in assisting the translation initiation of bymovirus precursor protein (Kanyuka et al. 2005;Stein et al. 2005), and HvPDIL5-1 is speculated to function as chaperone in correct folding of virus proteins (Yang et al. 2014b). Our previous work revealed sequence diversity of HvEIF4E in domesticated barley accessions and HvP-DIL5-1 in wild/domesticated accessions (Hofinger et al. 2011;Yang et al. 2014b). However, the sequence diversity of HvEIF4E in wild barley was unknown, and importantly, the responses of most of these HvEIF4E and HvPDIL5-1 haplotypes to the bymovirus disease of barley were not yet investigated. This limited our understanding of the relevance of naturally occurring sequence variation of both host factor genes in eco-system vs agricultural system, especially in the regions severely threatened by the bymovirus disease. By combining previously published and newly generated sequences and by testing with a viral assay the identified haplotypes, here we systematically investigated the diversity and frequency of sequence variation patterns in resistanceconferring haplotypes of these two host factor genes in barley genotypes from different geographic regions.

Plant materials
365 wild and 2557 domesticated barleys (mostly landraces or historic cultivars, but also including six-rowed barley with brittle rachis Hordeum vulgare subsp. agriocrithon) with the GPS-coordinated geographical information, which were mainly collected by IPK-Genebank or obtained from the National Small Grains Collection, USDA-ARS (Aberdeen, Idaho, USA), were used in the present study (Supplementary Fig. S1; Supplementary Tables S1, S2). The wild barley accessions were mainly collected from the Near East, and the cultivated barley accessions originated from 89 countries representing most of the barley cultivation areas worldwide. For a small number of accessions the precise collection site was unclear, and for these we specified the latitude and longitude of the capital of the source state/province. Details of the barley accessions and the collection sites (Country, latitude and longitude) are given in Supplementary Table S2.

Amplicon sequencing
Plant DNA and RNA extraction, synthesis of first-strand cDNA, PCR amplification, purification of PCR products and Sanger sequencing were performed as described previously (Yang et al. 2014b). Purification of PCR amplicons was carried out using the NucleoFast 96 PCR Kit (Macherey-Nagel, Germany). Two µl of each purified PCR product was analyzed and quantified by agarose gel electrophoresis. The optimal amount of each purified amplicon (10 ng per each 100 bp) was subjected to cycle-sequencing and the products analyzed on an ABI-3730xl DNA Analyzer (Applied Biosystems, Darmstadt, Germany). The trimmed sequence reads were assembled using Sequencher 4.7 (Gene Codes, USA). The primers used in this study are given in Supplementary Table S3.

Resistance tests
Plant material required for performing virus resistance tests was cultivated in climatic chamber as previously described (Habekuß et al. 2008;Yang et al. 2014b). A subset of wild and domesticated barley accessions carrying non-synonymous mutations in either HvPDIL5-1 or HvEIF4E was analyzed after mechanical inoculation with the Aschersleben isolate of BaMMV (BaMMV-ASL). At least 10 plants per accession were tested. At three-leaf stage (2 week after sowing) plants were mechanically inoculated twice at an interval of 5-7 days. Four to six weeks later, the phenotypes were collected by detecting viral RNAs by RT-PCR (via amplification of a 525-bp viral VPg gene fragment), or DAS-ELISA with a viral coat protein-specific antibody (Clark and Adams 1977). For instance, if absorbance at 405 nm was below 0.1, the respective plant was scored as being resistant (no virus particles detected). Otherwise, plants were scored susceptible (absorbance >0.1).

Data analysis
Sequence assembly and alignments were performed using Sequencher v4.7 (Gene Codes, USA). Sequence reads that originated from introns, and those of poor quality or indication of heterozygosity were removed manually from the analysis. The genotypes for which only the partial gene coding sequences were obtained were also excluded from the analysis. The polymorphic loci including loss-of-function (LoF), synonymous (S) and non-synonymous (Ns) were counted manually. The allelic haplotypes, the haplotype diversity (H), the nucleotide diversity (π) as well as statistics for selection (Fu and Li's D* and F*, Tajima's D) were calculated using DNASP v5.10.01 (Librado and Rozas 2009). The generated files, containing polymorphism and haplotype information, were subjected to the Median-Joining network analysis using DNA Alignment v1.3.1.1, Network v4.6.1.1 and Network Publisher v1.3.0.0 software (Fluxus Technology, UK). The significance of type of mutations between wild and domesticated barleys was analysed using Fisher's Exact Test.
The topographic maps were drawn using ArcGIS 10 software (ESRI, Redlands, CA, USA). The Base Map at 10-km resolution developed by the US Geological Survey (USGS) according to the GTOPO30 global digital elevation dataset (http://eros.usgs.gov) was used. These accessions used in this study were assigned to the Base Map according to the specified collection sites (Supplementary Table S2).

HvPDIL5-1 and HvEIF4E exhibit distinct patterns of sequence diversity in domesticated but not in wild barley
To assess the extent and pattern of global genetic diversity of the two host factor genes HvPDIL5-1 and HvEIF4E, sequence variation in the full-length open reading frames (fl-ORFs) of both genes was surveyed in 365 wild and 2557 domesticated geographically referenced barley accessions (Supplementary Tables S1, S2). The dataset included (i) previously published sequences of HvEIF4E obtained from 1090 cultivated barleys (Hofinger et al. 2011) and HvPDIL5-1 obtained from 350 and 1382 wild and domesticated barleys, respectively (Yang et al. 2014b), as well as (ii) newly generated sequences of HvPDIL5-1 and HvEIF4E obtained in this study from 856 cultivated and 320 wild barley accessions, respectively (Table 1). Furthermore, we sequenced in a subset of the total of 2922 barley accessions the fl-ORFs of HvGT43 and HvMCT-1 (Supplementary Table S1). These genes co-segregated and are located in close physical proximity to HvPDIL5-1 and HvEIF4E, respectively (Fig. 1a, Stein et al. 2005;Yang et al. 2014a, b), and were included in this study for the purpose of determining the demography of population and testing the effect of genomic location on patterns of sequence diversity. The relative genomic positions of HvGT43/HvPDIL5-1 and HvMCT-1/HvEIF4E, besides allocating to different chromosomes, are significantly different. The first gene pair is located close to the centromere (low recombination frequency region) on the long arm of chromosome 4H (4HL) whereas the second gene pair is located close to the telomere (high recombination frequency region) on the long arm of chromosome 3H (3HL) (IBRC 2012). Importantly, fl-ORF sequences of HvPDIL5-1 and HvEIF4E genes were obtained from a shared set of 1152 genotypes, and fl-ORF sequences of all four genes HvPDIL5-1, HvEIF4E, HvGT43 and HvMCT-1 were obtained from a shared set of 392 barley accessions (Supplementary Table S2).
We further analyzed the coding sequence diversity of these four genes in wild barley in comparison to domesticated barley.  (1) and (2) This study

HvPDIL5-1
Nineteen wild barley haplotypes (350 analyzed sequences) were observed based on ten non-synonymous (refers to non-synonymous amino acid change or in-frame deletions/ insertions) and eight synonymous mutations, while only 15 haplotypes in domesticated barley (from 2238 sequences) were revealed with five containing loss-of-function mutations/deletions, another five containing non-synonymous mutations and remaining four containing synonymous mutations (Table 1). Synonymous and non-synonymous mutations in wild barley occurred evenly along the coding gene sequence, whereas loss-of-function sequence variation in domesticated barley were enriched in the central part of the gene (Fig. 1b), coding putatively for the functional thioredoxin domain (Yang et al. 2014b). Four haplotypes (hap-I, III, IV and XXIV) were shared between wild and domesticated barleys (Fig. 2). Haplotype I predominated in the collection and other haplotypes were derived from this ancestral haplotype by one or two mutations (Yang et al. 2014b). Interestingly, loss-of-function haplotypes were only found to be significantly accumulated in domesticated, but not in wild barley (Fisher's Exact Test, P = 0.0099). The number of polymorphisms, haplotypes, as well as the values of the parameters H and π were lower in domesticated barley compared to those in wild barley (Table 1), indicating the reduction of gene sequence diversity. Statistical tests employing Fu and Li's D* and F* and also Tajima`s D indicated low frequencies of rare alleles of HvPDIL5-1 in wild and domesticated barley, possibly being related to selection or population size expansion.

HvGT43
The close neighbor gene of HvPDIL5-1 exhibited thirteen and five polymorphic nucleotide positions (including nonsynonymous and synonymous mutations) in 192 wild and 221 domesticated barley accessions, respectively, defining 14 and seven haplotypes (Table 1). Since the same three major haplotypes (hap-I, II and III) were also present in domesticated barley accessions (95.5%, 211 out of 221 accessions), no strong differences were found for H and π values between the domesticated and wild barleys ( Table 1). The same set of cultivated barley vs wild barley revealed a strongly decreased sequence diversity at the HvPDIL5-1 locus (Supplementary Table S9). This indicated that the observed sequence diversity of HvPDIL5-1 and HvGT43 was gene specific and irrespective of chromosomal position. Importantly, no defective haplotype was found at the HvGT43 locus.

HvEIF4E
In 320 wild and 1090 domesticated barley accessions, 26 and 47 haplotypes were defined, respectively. Fourteen nonsynonymous and eleven synonymous nucleotide polymorphisms were detected in wild barley (Table 1), and these were evenly distributed along the coding gene sequence (Fig. 1b). This pattern was similar to the observations at the HvPDIL5-1 locus, which implied a lack of specific selection at a particular region of both genes in wild barley. As to domesticated barley, an extreme situation was observed: sequence diversity was caused by a total of 30 non-synonymous mutations, whereas not even a single silent (synonymous) mutation was detected (Fisher's Exact Test, P = 0.0001, Table 1). Non-synonymous mutations were enriched at three regions (nucleotides 336-384, 480-528, and 600-648) (Fig. 1b), coding for amino acid sequences located in proximity to the putative cap-binding domain of the protein (Kanyuka et al. 2005;Stein et al. 2005). In contrast to HvPDIL5-1 that showed a simple evolutionary network, the Median-Joining (MJ) network revealed a complex status of sequence variation for HvEIF4E (Fig. 2). For instance, it seems likely that haplotype XII originated by one additional mutation from one of haplotypes wt0A, XIV, XXI, XXXI or XXXII. The number of polymorphic loci and haplotypes, and H and π values in cultivated barley were not lower than in wild barley.

HvMCT-1
Three synonymous and three non-synonymous mutations of this gene were detected in wild barley accessions, while only three non-synonymous mutations were found in domesticated barley (Table 1). In contrast to HvEIF4E that represented a higher sequence diversity in domesticated barley, the same set of barley accessions revealed  Table  S9). This indicated that the higher sequence diversity of HvEIF4E was not a function of recombination frequency at the 3HL telomeric region. Collectively, the sequence variation analyses revealed two qualitatively distinct patterns of gene sequence diversity with a bias for loss-of-function mutations in HvP-DIL5-1 and non-synonymous mutations in HvEIF4E-both occurring only in domesticated but not in wild barley accessions. Importantly, analyzing the sequence diversity of the next neighboring genes HvGT43 and HvMCT-1 indicated that the two distinct patterns observed for HvPDIL5-1 and HvEIF4E are characteristic for the respective positions on the chromosome.
Bymovirus resistance correlated with the increased diversity of the host factor genes HvEIF4E and HvPDIL5-1.
The barley yellow mosaic virus disease is a highly virulent and persistent biotic stress for cultivated barley in East Asia and Europe (Kühne 2009), and selection for resistance was an important target in recent barley breeding and cultivation in these regions (Chen 2005;Ordon et al. 2005). To determine whether appearance of haplotypes conferring virus resistance was associated with the occurrence of loss-of-function and non-synonymous exchange mutations in HvPDIL5-1 and HvEIF4E (the corresponding recessive alleles of these genes conferring resistance to bymoviruses are known as rym1/11 and rym4/5, respectively), we tested accessions carrying different haplotypes at the HvPDIL5-1 or HvEIF4E locus for their reaction to a common isolate of Barley mild mosaic virus (BaMMV-ASL). A total of 30 and 65 different haplotypes of HvPDIL5-1 and HvEIF4E, respectively, were defined in wild and domesticated barley accessions (Supplementary Tables S4, S5). Of these, a subset of accessions representative for haplotypes carrying either lossof-function and/or non-synonymous mutations were subjected to mechanical virus inoculation (Supplementary Tables S10, S11). In these assays, primarily the accessions carrying known, previously characterized susceptibility haplotypes at one locus and novel haplotypes with the unknown susceptibility/resistance status at another locus were used. For instance, when one accession contained a known susceptibility haplotype of HvPDIL5-1, testing the susceptibility/resistance status of a novel haplotype of HvEIF4E was allowed. However, most of accessions carrying haplotypes of unknown susceptibility/ resistance status at both loci (HvPDIL5-1 and HvEIF4E) were excluded (or the phenotype of respective haplotype was marked as 'unknown'). It needs to be noted, that this approach cannot rule out completely the remote possibility that an accession may contain additional independent bymovirus resistance loci. Less than 1% of nearly 10,000 domesticated barley accessions were reported to be completely resistant against isolates of BaYMV and BaMMV in the field in earlier studies (Ruan et al. 1984;Zhou and Cao 1985).

HvEIF4E
Accessions representing either one of the identified nineteen HvEIF4E haplotypes were tested as conferring resistance to BaMMV-ASL (Supplementary Tables S5, S11). This included three previously reported resistance alleles (rym4, rym5, and rym4 HOR4224 ) (Perovic et al. 2014;Stein et al. 2005), and 15 new haplotypes for which bymovirus resistance had not been reported before (hap-IX, X, XIII, XIV, XV, XVII, XIX, XXIII, XXIII, XXIV, XXVII, XXVIII, XXXI, XXXVII, XLI, and XLV) (Fig. 3). Wild barley accessions carrying HvEIF4E hap-XLV were also resistant to BaMMV-ASL (Supplementary Table S11), confirming that H. spontaneum is a useful source of resistance to bymoviruses. In domesticated barley, sixty-eight (5.31%) of 1090 accessions carried any of the remaining 18 resistance-conferring HvEIF4E haplotypes (Fig. 3). These haplotypes were found in many barley cultivation areas except the Americas, Near East and Oceania (Fig. 4b), where no reports about the occurrence of the bymovirus disease are available. Interestingly, all resistance-conferring haplotypes contained non-synonymous mutations in the coding sequence of HvEIF4E.
Collectively, the data presented above revealed a strong correlation between bymovirus resistance and either the loss-of-function mutations in HvPDIL5-1 or the non-synonymous mutations in HvEIF4E. Importantly, resistanceconferring haplotypes were found preferentially in cultivated barley rather than wild barley, suggesting the more frequent rise of bymovirus resistance after domestication.
Type and frequency of mutations in the two host factor genes in context of geographical origins of the domesticated barley accessions.
The large collection of geographically referenced domesticated barley accessions allowed us to investigate if any geographic region contributed disproportionally to gene sequence diversity (loss-of-function mutations of HvPDIL5-1 and non-synonymous mutations of HvEIF4E), and if there were any preferred geographic origins of resistance-conferring haplotypes. Seven sub-populations (Africa, Americas, Central Asia, East Asia, Europe, Near East and Oceania) of the domesticated barley accessions were defined (Supplementary Table S1), and analyzed for diversity in the host factor genes HvPDIL5-1 and HvEIF4E (Tables 2, 3).

HvPDIL5-1
The polymorphic loci, the number of haplotypes, H and π of HvPDIL5-1 decreased in all seven sub-populations of domesticated barley compared to those in wild barley (Tables 1, 2). Six HvPDIL5-1 haplotypes (hap-II, VII, VIII, IX, XXVIII and XXIX) were found to be associated with bymovirus resistance, with five of those containing a lossof-function mutation/deletion and one other carrying a nonsynonymous mutation. One of these resistance-conferring haplotypes of HvPDIL5-1, rym11-a (hap-XXVIII), carrying a 1375-bp deletion was found only in one accession from Near East (Yang et al. 2014b). The other five haplotypes were found in 36 accessions originated from East Asia and one accession (Russia 57, hap-II, rym11-b) for which the exact location of a collection site was uncertain (Fig. 3). Of these, 35 contained haplotypes with loss-offunction mutations, indicating a disproportional frequency of resistance-conferring haplotypes in East Asia compared to other barley cultivation regions (Fig. 4a). The frequency of these haplotypes varied between 0.23% (e.g. hap-VIII) and 6.56% (hap-II) in the set of 442 accessions from East Asia used in this study. In comparison, no accessions carrying resistance alleles were found in Europe, where the bymovirus disease is also known to occur (Kühne 2009).

HvEIF4E
The genetic diversity of HvEIF4E varied among cultivated barley from the different regions of origin (Table 3). The highest diversity was observed in barley from East Asia. Most of the bymovirus resistance-conferring haplotypes were found exclusively in a single barley cultivation area, except hap-XXII and XXVIII, which were present in more than two major geographic regions (Fig. 3). Two haplotypes, rym4 and rym5, were extensively found in Europe and East Asia, respectively (Fig. 3). This is consistent with the known origin of barley landraces Ragusa (from Croatia) and Mokusekko 3 (from China) that were originally identified as carriers of rym4 and rym5, respectively, and used for introgression of bymovirus resistance into elite barley germplasm by breeding (Ordon and Friedt 1993). In East Asia, fourteen resistance-conferring haplotypes were found (Fig. 3). They were present in 46 (24.2%) of 190 cultivated barley accessions from this region in our collection (Table 3). Each haplotype was present at relatively low frequency in the set of accessions tested in this study, ranging between 0.53% (1 of 190, e.g. hap-rym5) and 8.95% (14 of 190, hap-XIV). By contrast, in Europe only three haplotypes (hap-rym4, XL and XXII) conferring resistance to bymoviruses were found. These haplotypes were present in 15 out of 304 (4.93%) barley accessions from this region (Table 3).
Taken together, the bymovirus resistance-conferring haplotypes caused by loss-of-function mutations in HvP-DIL5-1 and non-synonymous mutations in HvEIF4E were greatly overrepresented in accessions from East Asia (Tables 2, 3), and the predominant occurrence as minor haplotypes together with statistics indicating selection suggest an evolution of bymovirus resistance alleles in barley germplasm from this geographic region.

Discussion
This work, by combining previously published data and newly generated sequence information, revealed sequence variation in a broad collection of wild and domesticated barley accessions of the two host factor genes, HvEIF4E and HvPDIL5-1, involved in bymovirus infection. Importantly, by testing haplotypes for resistance, this work demonstrated the two distinct types of sequence diversity patterns occurring in resistance forms of the two host factor genes: non-synonymous mutations in HvEIF4E and lossof-function mutations in HvPDIL5-1. Re-sequencing of Fig. 4 Collection sites of the barley accessions carrying resistance-conferring haplotypes of HvPDIL5-1 (a) and HvEIF4E (b). Each accession is allocated according to latitude and longitude coordinates of its geographic collection site (represented as red or blue circles). Red and blue circles indicate collection sites for wild and domesticated barleys, respectively. Size of each circle represents the number of accessions collected at each particular collection site (color figure online) the neighbor genes HvGT43 and HvMCT-1 in a subset of barley accessions indicated that the two distinct patterns of preferential resistance-conferring sequence polymorphism of the two host factor genes were not caused by population demographic history or their chromosomal positions on the respective chromosomes. Haplotypes associated with bymovirus resistance occurred predominantly in domesticated barley accessions but were rare in wild barley, indicating a rapid rise of resistance in cultivated barley post domestication. Interestingly, sources of bymovirus resistance were mainly found in barley germplasm from East Asia where the barley yellow mosaic virus disease is omnipresent (Kühne 2009), indicating the existence of interaction between bymoviruses and host factors. In addition, although the overall resistance spectrum of the identified HvEIF4E and HvPDIL5-1 haplotypes conferring resistance to BaMMV-ASL warrants further investigation (e.g. by testing using a range of BaMMV and BaYMV isolates), these haplotypes represent useful new sources of bymovirus resistance for exploitation in barley breeding programs globally.
The differences in gene diversity patterns observed for HvEIF4E vs HvPDIL5-1 may be explained by the different functions of these two host factor genes.
The eukaryotic translation initiation factor 4E (EIF4E) was repeatedly reported as a host factor involved mainly in potyvirus infections in different plant species (Charron et al. 2008;Wang and Krishnaswamy 2012). The EIF4E protein is of crucial importance in organisms, as it is a component of the conserved eukaryotic translation initiation complex that interacts with the cap-structure of the host mRNAs and enables initiation of the translation of cellular proteins (Browning 1996;Kawaguchi and Bailey-Serres 2002). In previous EcoTILLING (a method for the discovery of natural nucleotide diversity) studies, no haplotypes carrying loss-of-function mutations in EIF4E were found in broad collections of accessions of a number of different crops and related wild plant species (Charron et al. 2008;Hofinger et al.2009Hofinger et al. , 2011Ibiza et al. 2010;Konecna et al. 2014;Nieto et al. 2007;Rubio et al. 2009;Stein et al. 2005;Stracke et al. 2007). In the coding sequence of HvEIF4E, only non-synonymous mutations were previously identified in an Ethyl Methane Sulfonate (EMS) induced barley population (Gottwald et al. 2009). Consistent with this, a previous EcoTILLING study using cultivated barley also revealed only non-synonymous and small in-frame insertion/deletion mutations in coding sequence of the HvEIF4E gene (Hofinger et al. 2011). Similarly, in this study surveying the variation in HvEIF4E also in a wild barley collection detected only synonymous and non-synonymous mutations and no loss-of-function mutations. This cumulative evidence strongly suggested that HvEIF4E is indispensable in barley, although a highly conserved paralogous gene,  HvEIF(iso)4E, exists in barley and may potentially provide some functional compensation upon loss-of-function of HvEIF4E (Yang et al. 2013). Since the VPg protein of plant viruses in the Potyviridae family plays a role similar to that of the eukaryotic mRNA cap-structure and competes with the plant cellular mRNAs for EIF4E binding (Leonard et al. 2000), the changed compatibility by amino acid changes may block the physical interaction with VPg resulting in disease resistance (Charron et al. 2008). This mechanism reported for potyviruses in an interaction with Arabidopsis may explain the situation in barley, since the amino acid substitutions in either HvEIF4E or the bymovirus VPg protein may influence the efficiency of interaction between these proteins (Habekuß et al. 2008;Kanyuka et al. 2005;Kühne et al. 2003;Stein et al. 2005). Thus, we speculate that amino acid substitutions in HvEIF4E represent the optimal way of evolution of this gene, to maintain the function for protein translation initiation in barley while simultaneously conferring resistance to bymoviruses. In contrast to the complex evolutionary haplotype network of HvEIF4E in domesticated barley, the HvPDIL5-1 haplotype variation was rather simple. Only 15 haplotypes were found in 2238 cultivated barley accessions. Five haplotypes encoded either no transcript (hap-XXVIII, 1375-bp deletion including first three protein coding exons leads to complete lack of transcript (Yang et al. 2014b), truncated proteins (hap-VII, VIII and XXIX) or shifted coding frames (hap-II). Thus all five identified bymovirus resistance-conferring haplotypes of HvPDIL5-1 lead to non-functional proteins or no protein at all. HvPDIL5-1 belongs to a protein disulfide isomerase (PDI) gene family (Yang et al. 2014b). The number of the family members varied in eukaryotic species. Arabidopsis, rice and maize contain 22, 19 and 22 PDIs, respectively (Houston et al. 2005). Thirty-five barley cultivars from East Asia contained one of these loss-of-function versions of HvPDIL5-1, suggesting that the non-functional alleles unlikely resulted in yield penalties. In addition, the EMS-induced non-functional alleles (rym11-9699 and -10253) of HvPDIL5-1 were not associated with any visible phenotypic changes in plant growth or development (Yang et al. 2014b), suggesting a certain level of dispensability for HvPDIL5-1 in barley. Although an amino acid substitution (rym11-e) of HvPDIL5-1 can confer bymovirus resistance in domesticated barley, the resistance spectrum of this allele to multiple isolates of bymoviruses still remains to be determined (Yang et al. 2014a). Overcoming resistance conferred by this allele by new bymovirus isolates may be rather likely in contrast to the other loss-of-function alleles (rym11-a, b, c, d, and f), but this is similar to the situation observed for several HvEIF4E alleles (Kühne et al. 2003;Stein et al. 2005). Complete knock-out of the HvPDIL5-1 gene likely terminates the co-evolution between bymoviruses and this host factor gene.  Therefore, we postulate that the specific biological functions of both host factors HvPDIL5-1 and HvEIF4E has driven the two observed distinct patterns of gene diversity. Importantly, the highest diversity of the resistance alleles for both host factor genes was found in East Asia, a region with high incidence of barley yellow mosaic virus disease (Kühne 2009).

Did natural selection for bymovirus resistance also occur in European cultivated barley?
The barley yellow mosaic virus disease occurs in a majority of European countries (Kühne 2009). Given the predominance of resistance-conferring alleles in East-Asian accessions, it is unclear whether natural selection for bymovirus resistances also independently happened in Europe, perhaps even before the active breeding for bymovirus resistance in barley had been initiated (since 1980s). No HvPDIL5-1 haplotypes conferring bymovirus resistance (rym1/11) were found in the current collection of European barley accessions, supporting the hypothesis that rym1/11 evolved in East Asia. Three HvEIF4E haplotypes (hap-rym4, XXII, and XL) that also confer bymovirus resistance were found in 15 of 304 European barley accessions. The eleven historical accessions (landraces) containing hap-rym4 or hap-XL were particularly collected before 1970 s from the Balkans (not shown), a region where the bymovirus disease was not reported until the early 1990s (Katis et al. 1997). The remaining four accessions carrying hap-XXII were from countries of Western and Northern Europe. The same haplotype (hap-XXII) was found also in Central and East Asia. In some geographic regions, where the barley yellow mosaic virus disease was not reported, twelve accessions (including wild and domesticated) showing bymovirus resistance were found. While there was an enrichment of HvEIF4E haplotypes conferring bymovirus resistance in European barley germplasm, from the results of the present study we cannot confidently conclude whether the resistance-conferring mutations in hap-XXII, and hap-XL originated in Europe or were introduced to Europe from elsewhere with the seed trade.


extfieldsofstudy: Biology Mathematics Physics Medicine

text:
Law of localization in chemical reaction networks

In living cells, chemical reactions are connected by sharing their products and substrates, and form complex networks, e.g. metabolic pathways. Here we developed a theory to predict the sensitivity, i.e. the responses of concentrations and fluxes to perturbations of enzymes, from network structure alone. Responses turn out to exhibit two characteristic patterns, $localization$ and $hierarchy$. We present a general theorem connecting sensitivity with network topology that explains these characteristic patterns. Our results imply that network topology is an origin of biological robustness. Finally, we suggest a strategy to determine real networks from experimental measurements.

I. INTRODUCTION
Cells have many chemical reactions, each of which is mediated by organic catalysts, enzymes. Reactions are not independent but connected and form complex networks.
The dynamics of chemical concentrations are considered as the origin of physiological functions. However, dynamical behavior based on the network is not understood well.
One experimental approach to study such network systems is sensitivity analysis where the amount or activity of enzyme is perturbed and responses (concentrations of chemicals in the system) are measured [1] (see Fig. 1). However, the results of such experiments are very difficult to interpret, because theoretical criteria to evaluate the results of perturbations from network structures are not established. There are other difficulties in understanding dynamical behaviors of reaction systems in biology. First, although a huge amount of information of reaction networks is available on databases [2][3][4], they provide no more than knowledge of identified reactions in biochemistry. It is possible that the information is incomplete, including many unidentified reactions or regulations. Second, in spite of the recent progress in biosciences, it is still difficult or almost impossible to determine quantitative details of the dynamics, such as functions for reaction rates, parameter values, or initial states.
In order to circumvent these difficulties, we introduce a mathematical method, named structural sensitivity analysis [13,14], to determine responses of chemical reaction systems to the perturbation of the enzyme amount or activity based only on network structure.
From analyses we found that qualitative responses at a steady state are determined from information of network structure only. We also found that response patterns, e.g. distribution of nonzero responses of chemical concentrations in the network, exhibit two characteristic features, localization and hierarchy depending on the structure of networks and position of perturbed reactions. Finally we found a general theorem connecting the network topology and the response patterns directly, and governing the characteristic patterns of responses. This theorem, which we call the law of localization, is not only theoretically important, but also practically useful for examining real biological systems.
In the context of adaptation, there were some previous studies, which reported confined nonzero responses in specific systems [9][10][11][12]. However, they did not find general laws of such response patterns, nor any topological conditions.

II. STRUCTURAL SENSITIVITY ANALYSIS
We study concentration changes in a reaction system under perturbation of reaction rate parameters, assuming that the system is in a steady state [13,14]. We label chemicals by m (m = 1, . . . , M ) and reactions by j (j = 1, . . . , R). A state of the system is specified by concentrations x m (t) and obeys the following differential equations [15,16] Here, ν is called a stoichiometric matrix. W i is called a flux, which depends on metabolite concentrations and also on a reaction rate k i . We do not assume specific forms for W i , but assume that each W i is an increasing function of its substrate concentration; Below, we abbreviate and emphasize nonzero ∂W i ∂xm as r im . In this framework, enzyme knockdown of the j-th reaction corresponds to changing the rate as k j → k j + δk j (triangles in FIG. 2). By assuming steady state [5][6][7][8], the flux is expressed, in terms of a basis { c n } of ker ν, as W = N k n=1 µ n c n , where N k is the dimension of the kernel and µ n are N k coefficients depending on reaction rates. Under the j-th knockdown, we have The i-th component of δ j W is also expanded as From (3) (4), the response of steady state concentration δ j x ≡ d x dk j δk j (circles in FIG. 2) and flux δ j W (arrows in FIG. 2) to each perturbation k j → k j + δk j is determined from network structure only [13,14]. The result for each perturbation is simultaneously obtained through the following matrix computation: where the matrix A is given as In (5) and (6), the horizontal and vertical lines are the partitions of the matrices. We then obtain the flux change δ j W from (3), or in a matrix notation. We call the inverse of A as the sensitivity matrix S. Note that δ j x, δ j µ, δ j W , c n are column vectors with M, N k , R, and R components respectively, and ∂W i ∂xm is an R-by-M matrix. We assume networks with dim ker ν T = 0 throughout this paper, which guarantees the matrix A is square, i.e. R = M + N k .
Comments are in order. First, our theory depends only on the structure of reaction networks. The network structure is reflected in the distribution of nonzero entries in the A-matrix, which determines the qualitative responses. Second, as a generalization of our method, we account for regulations such as allosteric effects by relaxing (2) as Then, regulations add additional r im in the A-matrix, and the response is still determined through (5).

III. LOCALIZATION AND HIERARCHY
Let us see some results of structural sensitivity analysis.
The flux changes only when we perturb the top reaction 1 (the 1st column of S). The perturbation to reactions 2 or 3 changes only its substrate concentration (the 2nd, 3rd column of S).
Again, only the perturbation to the input rate, corresponding to the 1st column in   Proof.-The theorem is proved from the distribution of nonzero entries of the A-matrix.
(i) Suppose a subnetwork Γ is a buffering structure. Then by appropriately choosing a basis of the kernel of ν and the orderings of the indices of the A-matrix, we can always rewrite the A-matrix as The lower left block vanishes because A is output-complete. (ii) As explained already, the concentration change δ j x m is proportional to A −1 mj ∝ DetÂ (j;m) , whereÂ (j;m) is the minor matrix associated with the row of the j-th reaction and the column of the m-th metabolite. Then, DetÂ (j;m) = 0 for i ∈ r, m / ∈ m follows because the upper left block in the minorÂ (j;m) , which was originally square in (12), is horizontally long.
We illustrate the law of localization in the example networks in FIG. 2 and FIG. 3 (Left). should be limited within Γ 4 , which is observed in the 4th column in (10). Similarly, the response to perturbation of reaction 6 is explained by Γ 5 . from network topology, surprisingly overlap biologically identified sub-circuits, the pentose phosphate pathway (yellow in FIG. 4, 5), the tricarboxylic acid cycle (blue) and the glycolysis (green). This correspondence may be understood from an evolutional point of view by considering the advantage of buffering structures.

V. DISCUSSIONS AND CONCLUSIONS
Here we discuss the biological significances of buffering structures (and nest of them) in two different levels. The first discussion is on the physiological importance. A buffering structure prohibits influence of given perturbation from expanding to the outside, like a "firewall." In other words, it is a substructure with robustness emerging from the network topology. The carbon metabolism network of E. coli possesses multiple nested firewalls (FIG. 5), and are expected to be robust to fluctuations of enzymes in it. We expect that such a topological characteristic of reaction networks could be the evolutionary origin of homeostasis in biological systems. A set of chemical reactions satisfying the condition of buffering structure by chance in evolutionarily early time would be positively selected as an advantageous circuit. We then expect that buffering structures in existing biological networks today might be generated and selected in such ways.