extfieldsofstudy: Biology Medicine

text:
Mathematical models of lignin biosynthesis

Background Lignin is a natural polymer that is interwoven with cellulose and hemicellulose within plant cell walls. Due to this molecular arrangement, lignin is a major contributor to the recalcitrance of plant materials with respect to the extraction of sugars and their fermentation into ethanol, butanol, and other potential bioenergy crops. The lignin biosynthetic pathway is similar, but not identical in different plant species. It is in each case comprised of a moderate number of enzymatic steps, but its responses to manipulations, such as gene knock-downs, are complicated by the fact that several of the key enzymes are involved in several reaction steps. This feature poses a challenge to bioenergy production, as it renders it difficult to select the most promising combinations of genetic manipulations for the optimization of lignin composition and amount. Results Here, we present several computational models than can aid in the analysis of data characterizing lignin biosynthesis. While minimizing technical details, we focus on the questions of what types of data are particularly useful for modeling and what genuine benefits the biofuel researcher may gain from the resulting models. We demonstrate our analysis with mathematical models for black cottonwood (Populus trichocarpa), alfalfa (Medicago truncatula), switchgrass (Panicum virgatum) and the grass Brachypodium distachyon. Conclusions Despite commonality in pathway structure, different plant species show different regulatory features and distinct spatial and topological characteristics. The putative lignin biosynthes pathway is not able to explain the plant specific laboratory data, and the necessity of plant specific modeling should be heeded.

Background
The recalcitrance of woody plant materials to enzymatic fermentation is the result of numerous molecular processes and features. At its core is the phenolic polymer lignin, which is interwoven with cellulose and hemicellulose, and thereby impedes access of cellulolytic enzymes, necessitating costly physico-chemical pretreatments before effective microbial fermentation can take place. With the exception of cellulose, lignin is the most abundant terrestrial biopolymer and accounts for roughly 30% of all organic carbon in the biosphere [1]. It gives a plant its structural stability, waterproofs the cell wall, thereby enabling water transport through the xylem, and protects the plant against pathogen threats. Lignin is an aromatic heteropolymer composed mainly of three types of hydroxycinnamyl alcohol monomers, namely the monolignols p-coumaryl alcohol, coniferyl alcohol, and sinapyl alcohol, which are commonly called H-, G-, and S-lignin, respectively.
Both the amount and composition of lignin are thought to be correlated with the hardness as well as the recalcitrance of structural plant materials. It is therefore important to the production and manipulation of bioenergy crops to understand the details of lignin synthesis and the deposition and polymerization of monolignols in the plant cell wall. In particular, the question arises whether

Open Access
Biotechnology for Biofuels *Correspondence: Eberhard.Voit@bme.gatech.edu 1 The Wallace H. Coulter, Department of Biomedical Engineering, Georgia Institute of Technology and Emory University, 313, Ferst Drive, Atlanta, GA 30332, USA Full list of author information is available at the end of the article it is possible and feasible to intervene in the phenylpropanoid pathway of lignin biosynthesis in a targeted and effective manner, for instance through gene knockdowns. The answer to this question is evidently preconditioned on detailed knowledge of this pathway and its control in situ. This knowledge in turn requires different types of biological data and, in cases where these are difficult to understand, the use of computational models that are capable of integrating small or large datasets of the same or different types and explaining observations that are sometimes unexpected.
In this article, we discuss computational models that are beneficial for explaining counterintuitive aspects of lignin biosynthesis and for making predictions regarding rational alterations in the molecular make-up of the pathway. We decided to present this material not in the form of a typical modeling paper, which would inform fellow modelers regarding all steps and technical details of model design, parameter estimation, methods of diagnostics and analysis, and interpretation of results. Instead, this paper is intended to address the practicing bioenergy scientist or engineer. It focuses on two overarching questions. First, what kinds and quantities of biological information are needed, or particularly beneficial, for setting up models of lignin synthesis and recalcitrance that have explanatory or predictive power? And second, if we succeed in constructing and implementing an effective model, what genuinely new insights might this model be able to offer? Guided by these questions, we will brush over most of the typical mathematical modeling steps and refer the reader to details in pertinent articles and reviews in the published literature.
It may surprise newcomers to the field of computational modeling that even within the limited scope of metabolic pathway modeling, the choices of mathematical formats and methods are all but unlimited. There is not "one" model that is somehow optimal, but there are many distinct options and numerous nuances. Even the representation of an enzyme catalyzed reaction can take a variety of mathematical formats, which are the result of different assumptions and focus either on molecular mechanisms or on the systemic behavior of a pathway system [2]. Studying these questions in detail, one comes to the conclusion that the selection of a model should ultimately be driven by the available data and by the scientific questions that the model is supposed to answer [3].
Although the structure of the lignin polymer is rather similar among different plant species, targeted experiments have revealed that the pathway of lignin biosynthesis exhibits variations among these same plant species. These variations are primarily manifest in the presence or absence of some of the involved enzymes, secondarily in different enzyme activities and substrate affinities, and third in possibly different regulatory control structures. As an illustration, Fig. 1 overlays the pathways of lignin biosynthesis in Populus trichocarpa (black cottonwood poplar), Medicago sativa (alfalfa), Panicum virgatum (switchgrass), and the model grass Brachypodium distachyon, as far as they are known or suspected today. A commonality among these species is that the pathway of lignin biosynthesis uses phenylalanine as its starting substrate; however, monocot grasses, including B. distachyon, and possibly P. virgatum as well, also use tyrosine, in addition to phenylalanine. It is presumably a biochemical necessity that most intermediates between these initial substrates and the final monolignols are by and large preserved, but the pathway systems in the species are connected in a slightly different manner through enzymatic reactions. These differences are not only of academic interest to the evolutionary biologist, but also of great significance to the biofuel researcher, because targeted interventions are almost always based on specific changes in gene expression, with concomitant alterations in fluxes through enzymatic reaction steps, such that a precise understanding of the details of the metabolic system is a prerequisite for targeted manipulations.
Predicting global effects of such manipulations on the ultimate lignin output and composition is not trivial, because the pathway utilizes the same enzymes for different reaction steps, but presumably with different substrate affinities (Fig. 1). Furthermore, the pathway is regulated, and some reactions occur in different locations of the cell and some may form functional metabolic channels. Details of the latter insights were actually derived from computational models that demonstrated that the absence of these features was inconsistent with experimental findings, as we will discuss later in this article.

An ideal dataset
In an ideal modeling world, experimental teams would be able to measure every piece of information needed to create a comprehensive model. The data would be of high quality, obtained in situ, from the same species and from multiple organisms. Obviously, this high bar cannot often be reached, and one must ask instead what compromises are still sufficient for modeling. We discuss this issue in the following.
To design and explore a model with computational methods, one needs to choose proper functional forms for the fluxes and determine their parameters. In a true mechanistic model, the mathematical format of a flux corresponds directly to the alleged biophysical or chemical mechanism, and typical parameters may be pH and temperature, and more specifically for metabolic models, may include quantities such as V max , K M , K cat , or K i , which correspond to rates and affinities in conceptual frameworks like the Michaelis-Menten mechanism.
In an idealized modeling situation, two scenarios can lead to a full model. First, knowledge of all metabolite concentrations and of all mechanisms, including input to the system, along with a complete set of physical and kinetic parameters, measured in vivo, can quite easily be converted into a comprehensive model. However, even in this quite unrealistic case, the model would ignore the spatial distribution of processes and stochastic events, which could, for instance, be due to environmental randomness or to very low numbers of enzyme or substrate molecules. Second, knowledge of all fluxes of the system and a complete set of measured physical parameters would allow the design of the model, again with the same limitation as before. At present, neither scenario is realistic, and missing information must be obtained from other sources, such as in vitro measurements, or inferred through computational means.
At this point, many modeling approaches and methods are readily available that could create functioning models out of such data, if they were available. However, they are not, and the more important point therefore is to realign the existing modeling techniques with the realities of data acquisition in a field where some of the key metabolic intermediates are below the level of solid quantification.
As a premier example, flux balance analysis (FBA) [4] and its extensions are based on a mathematical framework that allows assessments of the distribution of fluxes within a metabolic pathway at a steady state under the assumption of an alleged objective of the cell or organism, such as maximal growth, the maximal efflux of some metabolite, or the production of a compound like lignin. FBA formulates the operation of the pathway system as a so-called "linear programing problem" that optimizes the chosen objective, while satisfying biological constraints, such as non-negativity and maximal magnitudes of fluxes.
FBA is a computationally simple, yet powerful tool that has been widely used in many contexts, including plant systems. For instance, in a plant context, Paez et al. [5] analyzed biomass synthesis in Chlamydomonas reinhardtii under different CO 2 conditions, and Chang et al. [6] presented a genome-scale metabolic network model of the same organism. An interesting variation of FBA is the method of minimization of metabolic adjustment (MOMA) [7], which in a mutated organism tries to emulate a flux distribution that most closely mimics the wild type. Lee et al. [8] used MOMA to analyze data from knock-down experiments with genes associated with lignin biosynthesis in alfalfa.
While FBA and MOMA focus on the important distribution of fluxes at a steady state, dynamic modeling attempts to capture time-dependent changes in metabolites following any sort of perturbation. The hope is not only to understand short-term responses better, but also to capture regulatory features of the pathway system that are likely to become critical when the system is mutated. Expressed differently, FBA by and large assumes that everything in the organism remains the same, except for the mutated process and its direct derivatives, although it is to be expected that the organism will attempt to regain normalcy upon such a perturbation by evoking compensatory mechanisms. Thus, dynamic modeling is in principle more powerful but requires much more data support.
In the following, we describe case studies addressing lignin biosynthesis in different plants and with different methods. As stated before, we will focus primarily on data needs and different model uses.

Use of in vitro data
At present, metabolic modeling is far from having access to ideal comprehensive data obtained in vivo. To overcome this challenge, a common approach is the use of in vitro equivalents. An excellent example of this strategy in the context of lignin modeling is the work by Wang et al. [9], who constructed a dynamic model based on kinetic reaction and inhibition parameters of pathway enzymes in the black cottonwood, Populus trichocarpa. The authors derived 189 kinetic parameters associated with generalized Michaelis-Menten mechanisms, primarily in the form of K cat , K m , and K i of the 21 enzymes involved in monolignol biosynthesis. They also measured absolute enzyme quantities using mass spectrometry. Furthermore, the authors used a measured S/G ratio to quantify the input flux with a customized optimization algorithm. Such optimization methods are often needed in large-scale metabolic modeling, because the number of fluxes is typically greater than the number of metabolites, which creates a mathematical situation that cannot be directly solved. The information from their experiments allowed Wang's team to construct a fully parameterized model with estimated input flux, which they formulated as ordinary differential equations (ODEs). They were able to obtain the steady-state flux distribution and to investigate the effects of enzyme perturbations on lignin content and composition.
In principle, the well-established strategy used by Wang's team is excellent, as it leads to a fully dynamic model that permits explanations and predictions. The somewhat disconcerting issue is the use of in vitro data, which at present seems unavoidable, but leads to the following questions: (1) To what extent are in vitro data accurate and representative of the pathway behavior in vivo, and does enough in vivo information exist to validate the results of such models? In other words, it is unclear how to assess the reliability of these models. (2) It is clear that no biomathematical modeling effort can presently claim to have taken all components and modulators of a pathway into account. Thus, is it possible to ensure that all relevant information is present quantitatively to reproduce and explain in vivo observations? Or is it simply not feasible to reconstruct the complex in vivo cell environment with sufficient reliability from in vitro information? For example, Wang et al. did not include the enzyme caffeoyl shikimate esterase [10] in their poplar lignin model [9]; this enzyme was discovered as a new component in the lignin pathway while their studies were ongoing.
These concerns are not exaggerated and can even be found in a very detailed microbial investigation by Teusink et al. [11], which provides a good perspective in this regard based on the much simpler pathway system of glycolysis in baker's yeast, Saccharomyces cerevisiae. Specifically, these authors compared in vivo flux and concentration profiles with the results of a computational model that had been constructed based on the best available kinetic parameters obtained in vitro. Despite the authors' dedicated efforts to use the same yeast source and obtain measurements under the same assay conditions, the discrepancies between the model results and the observed in vivo behavior were alarming. For possible explanations, Teusink et al. pointed to potential factors that may be active in vivo and cause uncertainties that are almost impossible to implement in in vitro models. Some of these uncertainties are apparently not adjustable by tuning of rate constants or through modifications in the model structure, but may be due to complicated combinations of molecular interactions between the pathway metabolites and enzymes or agents outside the investigated metabolic pathway. The authors proffered that these small details might have caused drastic differences during the integration of in vivo information into systemic models.
Similar concerns about in vitro-in vivo extrapolations were voiced some while ago by Savageau and others [12,13], while Albe and Wright and others came to the conclusion that such an extrapolation is, by and large, justified in many cases [14,15]. In any case, these problems are disconcerting, as in vivo measurements are incomparably more difficult to perform than experiments in vitro. Then again, if it is only possible to obtain in vitro data, are there means of in vivo validation? It appears that a direct validation of individual process representations will be difficult. Thus, one must hope that different types of in vivo data may fill the gap as they are combined with models designed from in vitro data for more reliable results.

Use of limited in vivo data
Lignin synthesis in poplar At present, the total body of in vivo data is dwarfed by information obtained in vitro, and this situation is not likely to change any time soon. As a case in point here, the lignin pathway simply does not permit many concentration or flux measurements in vivo. Instead, the typical dataset that can reasonably be expected today consists of lignin content and composition under different conditions, possibly augmented with a few metabolite concentrations. Although this situation might seem to be quite dire for modeling, mathematical and computational approaches can still offer interesting results.
As a pertinent example, Lee and Voit [16] investigated the lignin biosynthesis pathway in Populus xylem based on a relatively small set of data consisting of the S/G ratios and down-regulation levels of enzymes in five transgenic plants. In addition, the authors used information regarding the pathway stoichiometry, regulatory information of five enzymes of the pathway, and an enzyme capacity measurement for COMT. Utilizing a predetermined lignin monomer composition as input and maximum lignin production as the cell's alleged objective, the authors were able to generate steady-state flux distributions through FBA methods. Furthermore, to convert this information into a dynamic model, they employed a strategy derived from biochemical systems theory (BST) [17][18][19][20][21]. In this modeling framework, all fluxes are represented with power-law functions, so that the parameters can be coarsely estimated without knowledge of direct measurements.
Once the model was fully parameterized, the authors were able to run simulations that ultimately reproduced the lignin composition in all measured transgenics. In fact, an entire ensemble of models was generated, rather than a single model with a unique set of parameter values. This ensemble of models was validated against two transgenics that had not been used to set up the model. Upon validation, an indirect optimization method [22] was implemented to propose enzyme profiles that were expected to lead to a minimal S/G ratio in order to minimize recalcitrance. Single, double, and triple enzyme alterations were conducted to give insights and to determine the most effective perturbations. An interesting detail to note is that the best triple mutation did not contain the double mutation plus an additional mutation, but a different set. Specifically, in comparison to the wild-type S/G ratio of about 1.8, the model predicted a minimal S/G ratio of about 1.3 for two modifications, namely reduction of COMT and CAld5H activities, but a minimum of about 1.1 for three modifications in which the activities of C4H, CAD, and CAld5H were somewhat increased. These computational predictions have not been tested in actual plants.
The fact that a model is able to predict the results of perturbations is intriguing, especially because the powerlaw representation does not explicitly model specific reaction mechanisms, but only the overall effect of a metabolite or regulator on a given process. Then again, the in vivo data used to formulate and instantiate the model encapsulate in some sense everything occurring in the plant, which is not the case for in vitro models. Encouragingly, the estimated parameters in Lee's analysis are in agreement with biochemical knowledge of the pathway and provide new insights into the dynamics of the pathway (see results in [16]). Similarly, the predictive capacity of the model to characterize the best candidates for gene alterations is interesting, but it remains to be seen whether explanations and hypotheses obtained with the model are comparable with those obtained with a model like Wang's [9], which was based on experimentally laborious in vitro data.
Lignin synthesis in alfalfa The structure of the lignin biosynthesis pathway and its regulation in alfalfa (Medicago sativa L.) are fairly well known, but some observations on transgenics were confusing as they seemed to contradict the pathway structure. In particular, some gene knock-downs led to different S/G ratios even though they occurred before the branch point where the pathways toward S-and G-monolignols diverge. Lee et al. [8] set out to investigate this situation, using an in vivo dataset of lignin content and composition in eight stem internodes in wild-type and seven transgenic lines (with reduced PAL, C4H, HCT, C3H, CCoAOMT, F5H, or COMT activity). The internode classification in this case provided the opportunity to characterize the differential biosynthesis of lignin during the maturation of stem tissue.
Without formal computation, an analysis of the logic of the pathway topology mandated the reversibility of the enzymatic steps catalyzed by HCT and C3H (Fig. 1), which had not been considered before. Taking this reversibility into account did not resolve the puzzle regarding S/G ratios though. Thus, the authors constructed a computational model of the pathway by first using FBA to compute the steady-state flux distribution in wild type, and then applying the method of MOMA [7] to analyze the redistribution of fluxes in transgenics. This analysis revealed that the results regarding S/G ratios in transgenics could not be explained unless functional channels were active to partition the pathway flux into dedicated S-and G-pathways.
Using statistical analysis, the authors showed that there was a strong correlation between the flux catalyzed by CCR1 and the flux of the consecutive reaction catalyzed by CAD in all strains except for the CCoAOMT-deficient line. This curious result indicated a lack of product exchange between coniferyl aldehyde produced by either COMT or by CCR1. To examine this situation more carefully, the authors tested the possibility of kinetic regulation by the CCR2-COMT and CCoAOMT-CCR1 routes ( Fig. 1), but extensive Monte-Carlo simulations indicated only a very remote possibility of kinetic regulation by substrate/product interactions. Instead, the analysis suggested regulation by one or more distant metabolites. The authors proposed that salicylic acid (SA) could act as the potential regulator of the pathway leading to S-lignin synthesis. Indeed, experimental data characterizing the correlation between SA and lignin content supported the computational hypothesis. Moreover, additional in vivo data, demonstrating the co-localization of COMT and F5H [23,24], provided further evidence supporting the channeling hypothesis.
Wang et al. [9] criticized Lee's approach on grounds that the method was rather indirect and, in particular, suggested that a complete kinetic model would be able to capture the experimental data without the need for channeling. While the existence of channels awaits further validation with direct experimental means, it is unclear whether a bottom-up kinetic approach would have led to the crisply targeted hypothesis of differentially regulated channels directing flux toward either S-or G-lignin.
In a different study, Lee et al. [25] investigated the channeling hypothesis in Medicago by setting up an ensemble of dynamic kinetic models in 19 pathway configuration variants. Each of these variants preserved mass conservation, while allowing alternative routes including one or two metabolic channels across coniferaldehyde (Fig. 2). The models also examined the presence or absence of putative regulatory mechanisms. Extensive Monte-Carlo simulations over a biologically meaningful range of kinetic values identified only 6 among the 19 plausible configurations as feasible and demonstrated that only 4 out of 16 combinations of plausible regulatory mechanisms could match the experimental data. A graph analysis of these six configurations showed that they were topologically closely related and corresponded to a closed network, if closeness between two configurations was defined as a difference in only one enzymatic reaction. Interestingly, all six feasible configurations in the analysis included one or both proposed metabolic channels.
While the computational results strongly suggest the existence of channels, and independent experimental evidence supports these results [8,23,24], it is of course imaginable that other explanations could be found for the counterintuitive data in alfalfa, because even the best model fit to data can never offer a guarantee that the model is in some sense correct or that there could not be other models satisfying the same data in a similar manner. It is interesting though, that the computational results were inferred directly from actual data from these same species and with a minimum of assumptions, whereas models based on in vitro data, obtained from bacteria, should be validated in the target species in situ, before they can be considered true. Furthermore, while the power-law formulation used by Lee is mathematically guaranteed to be correct at an operating point of choice, there is no such guarantee for Michaelis-Menten functions; in fact, it is clear that their underlying assumptions and prerequisites are violated in situ [2,12,26].
Lignin synthesis in switchgrass Similar to the investigations on poplar and alfalfa, a limited dataset characterizing lignin content and composition was available for switchgrass (Panicum virgatum) [27], one of the most promising plants in bioenergy research. This dataset was used to set up a model of lignin biosynthesis and to examine for this species the hypothesis of channeling at a diverging branch point, leading to either S-or G-lignin. Specifically, wild-type and four transgenic (4CL, CCR, CAD, and COMT) lignin profiles were analyzed with FBA methods to compute steady-state flux distributions. The stoichiometric model included three variants permitting alternative, slightly differing pathways with and without a hypothetical metabolic channel comprising CCR and CAD. Extensive Monte-Carlo simulations generated thousands of random kinetic parameters to test whether any of the three configurations could reproduce the experimental  [25]. The two crossing channels are associated with coniferaldehyde (see Fig. 1) data in a dynamical manner. Surprisingly, none of the configurations was able to capture the increase in H-lignin in 4CL-transgenics. Instead, the computational results suggested the necessity to include product inhibition by downstream pathway metabolites, as well as substrate competition between CCR substrates. These computational suggestions identified p-coumaroyl-CoA and feruloyl-CoA as possible regulators that were arguably necessary to reproduce the observed increases in H-lignin. The model also revealed that the reaction catalyzed by 4CL, which converts ferulic acid into feruloyl-CoA, constitutes an impediment for explaining the counterintuitive accumulation of ferulic acid in COMT transgenics.
Further computational analysis suggested the accumulation of some so-far unidentified metabolite as an inhibitor of 4CL and as the mechanism by which ferulic acid increased. Revisiting the experimental data indicated a slight accumulation of p-coumaric acid and caffeic acid, which was shown to suffice to support the model-based hypothesis. Taken together, the pathway configuration including both the CCR-CAD channel and two independent CCR and CAD reactions, along with the deduced regulatory mechanisms, turned out to be the only structure capable of matching the in vivo data. The authors validated the model to some degree by testing the responses to an enzyme expression profile in an independent transgenic PvMYB4 line that had not been used at all to set up the model. Overall, the analysis produced satisfactory results with respect to lignin content and composition, as well as the concentration profiles of several of the pathway intermediates [27,28].

Use of pathway data and 13 C-labeling data in Brachypodium distachyon
This case study describes new results that have not been published so far. For this reason, a short description of methods is provided in a later section.
While the results of analyzing in vivo alfalfa and switchgrass transgenics data in a somewhat indirect manner were interesting and could be validated to some degree, the data themselves constitute a rather thin base for model development. This base becomes more solid if it is combined with other types of data. An example for such a merging of heterogeneous data types is the lignin biosynthetic pathway in Brachypodium distachyon. In contrast to dicots, monocot grasses use both phenylalanine and tyrosine as the initial substrate for monolignol production (Fig. 3). One puzzling aspect of this apparent redundancy is that, despite the nearly equal contribution of both precursors to the total lignin content, phenylalanine is preferentially incorporated into G-lignin, and tyrosine into S-lignin, although both pathways converge at the same intermediate metabolite, p-coumaric acid [29]. This result is surprising and cannot easily be explained with putative structure of the lignin pathway in Brachypodium. Beyond the existence of this intermediate, where the two pathways converge, the G-and S-lignin pathways appear to be the same until they split at the coniferaldehyde node.
A computational model directly corresponding to the alleged pathway structure (Fig. 3) confirms the logicbased analysis: the pathway, as currently alleged, cannot reproduce key observations, such as the differential channeling of phenylalanine and tyrosine toward G-and S-lignin. Specifically, model simulations demonstrate that the pathway scheme in Fig. 3 is unable simultaneously to satisfy the following observed requirements: • Match the amount of 13 C-labeled H-lignin in experiments with [U-13 C 9 ]phenylalanine; • Match the observed 13 C incorporation into ERbound ferulic acid in the same experiment; • Capture the differential 13 C incorporation levels from [U-13 C 9 ]phenylalanine and [U- 13 C 9 ]tyrosine in lignin units.
One great advantage of a modeling approach is the relative ease with which it is possible to test different hypotheses and variations of the pathway structure in order to obtain possible explanations. As a specific example, it was reported that the three enzymes C4H, C3′H and F5H of the lignin biosynthesis pathway in B. distachyon are bound to the outer surface of the ER, while the remaining enzymes are located freely in the cytosol ( [29]; unpubl. data). This finding led to the hypothesis that the spatial localization of enzymes might be a reason for the preferential incorporation of phenylalanine and tyrosine into different monolignols. This hypothesis was readily tested with a computational model that distinguishes the two locations (see below). These two locations, or compartments, are physically not strictly separated, but allow the handing over of metabolites through diffusion.
To test the hypothesis of two distinct locations, we set up a refined model scheme by assigning the reactions catalyzed by the ER-bound enzymes, C4H, C3′H and F5H, to the ER compartment, and all others to the cytosol compartment (Fig. 4). While there is no strict spatial separation between ER and cytosol, we assumed preferential enzyme activity within each compartment and slower diffusion between compartments. As a note, only the net diffusion fluxes are shown in the pathway model, but both forward and reverse diffusions are considered explicitly in the computational model (see later section). Specifically, we took the following steps for our model design. In the current scheme (Fig. 4), the only means for incorporation of 13 C into H-lignin is through the diffusion flux D 2 , and this flux is diluted with the influx from unlabeled tyrosine. To increase 13 C incorporation into H-lignin, a second diffusion flux, D 9 , is added between the ER compartment downstream of D 2 , and this flux compensates for the dilution of tyrosine (Fig. 5). This diffusion flux D 9 can be interpreted as partial activity of 4CL in the ER compartment.
Beyond the inconsistent amount of H-lignin, low incorporation of 13 C in ferulic acid in the [U-13 C 9 ]phenylalanine labeling experiment is an indication for dilution by unlabeled tyrosine through caffeic acid. Therefore, a downstream influx, D 10 , from the ER compartment is postulated to compensate for tyrosine dilution and to increase 13 C incorporation in wall-bound ferulic acid. Again, this flux corresponds to partial activity of C3H in the ER compartment (Fig. 5).
Closer inspection of the pathway reveals that the key site for preferential incorporation of [U-13 C 9 ]phenylalanine and [U-13 C 9 ]tyrosine into different lignin units is the branch point where the pathways toward G-and S-lignin diverge; this divergence happens at the coniferaldehyde node. The original scheme in Fig. 4 dictates the same level of 13 C-labeling into both G and S units, due to dilution in both compartments at the coniferaldehyde node into the free cytosol. To explain the actually observed higher incorporation of 13 C into G-lignin in the phenylalanine labeling experiment, an undiluted upstream flux from the ER is necessary to compensate for the dilution from the cytosol influx (D 5 and D 6 ) into the immediate G-lignin precursors coniferaldehyde and/or coniferyl alcohol. We first modeled this hypothesis by simply adding a suspected direct flux from p-coumaroyl-CoA into coniferyl alcohol (Fig. 5, thick blue arrow).
Simulations with this amended model showed that the scheme in Fig. 5 is able to capture the levels 13 C incorporation in H-lignin and ferulic acid from [U-13 C 9 ]phenylalanine experiments. Also, by acting as a metabolic channel, the direct flux from p-coumaroyl-CoA into coniferyl alcohol shields the flow within the ER compartment from strong dilution by diffusion from the cytosol compartment, and thereby enables the preferential  While the long metabolic channel in Fig. 5 is able to simulate the preferential incorporation of precursors into lignin units, it is intriguing to determine whether fewer enzymes in such channel could still reproduce the data. Therefore, we examined the scheme in Fig. 5 toward the shortest channel possible (Fig. 6). This analysis suggested that the critical point to shield the ER compartment from strong dilution by cytosolic diffusion fluxes is coniferaldehyde. Without this compound protected, Gand S-lignins cannot attain different 13 C 9 -labeling levels. If this conjecture can be validated, the simplest scheme consists merely of a CCR/CAD channel.
Simulations of the scheme in Fig. 6 resulted in steadystate flux distributions that capture the experimental 13 C-labeling data (Fig. 7). Phenylalanine and tyrosine contribute nearly equally to the resulting lignin content: in the [U-13 C]phenylalanine experiment, 35% of phenylalanine is labeled and tyrosine is unlabeled (natural abundance), while in the [U-13 C 9 ]tyrosine experiment, 35% of tyrosine is labeled and phenylalanine is unlabeled (natural abundance). The labeled fluxes in Fig. 7 compare the contributions of phenylalanine and tyrosine in each pathway flux. Figure 8 exhibits the total flux values, which combine the values of labeled and unlabeled fluxes. Since the magnitude of the input flux is unknown, we normalized the input to a base value of 100 units of mass per unit of time.
Because the system is mathematically underdetermined, its steady-state solution is not unique (see Steady-state analysis section). Therefore, a range of admissible steady-state values is possible for each flux. It is worth emphasizing in this context that all solutions in the resulting ensemble are consistent with all pertinent observations; namely: • Each model in the ensemble captures the experimental data with respect to the label distribution in steady-state fluxes. For instance, V 24 shows a higher labeled portion than V 28 when phenylalanine contains the feeding label; • The lignin compositions and S/G ratios in all scenarios are compatible with experimental data; • The labeled lignin composition is compatible with 13 C 9 -phenylalanine and 13 C 9 -tyrosine experimental data; and • The labeled ferulic acid and p-coumaric acid match with experimental data.
Further details are presented in Table 1.
The boxplots in Figs. 7 and 8 reflect the distributions of admissible values. As can be seen, V 8 admits small values in comparison to its parallel reactions in cytosol compartment, i.e., V 22 and V 23 . This result demonstrates that, while the main pathway for the reactions catalyzed by CCR and CAD resides in the cytosol, a relatively small and undisturbed flux through CCR/  CAD at the ER is sufficient to establish the metabolic channel necessary for preferential incorporation. In fact, considering the wrinkled environment of the ER surface, it is not hard to imagine that localized pools would keep a small fraction of the pathway undisturbed from exchanges of metabolites with the cytosol.

Generic model formulation
In a kinetic systems model, the dynamics of the pathway is represented by a system of ordinary differential equations (ODEs) in which the metabolites are the states. The rate of change in each metabolite is determined by sums and differences of all fluxes that directly affect this metabolite. Each flux is a mathematical function of the metabolites and other variables of the system that needs to be selected. Although the fluxes are usually nonlinear functions, the collection of fluxes itself forms a linear system, which can be represented as a matrix equation of the type Here, X is the vector of metabolites, Ẋ is its derivative with respect to time, S is the stoichiometric matrix, and V is the vector of fluxes. The stoichiometric matrix S defines the pathway structure. An element S i,j of this matrix equals 1 if flux V j is directed toward metabolite X i . It is − 1, if flux V j removes material from metabolite X i , and it is equal to 0, if flux V j has no direct effect on metabolite X i . In long form, the matrix equation can be rewritten for each equation as where n is the total number of fluxes.

Steady-state analysis
The steady state of a system is important for two reasons. First, many biological systems tend to operate close to such a state, where the overall concentrations of metabolites do not change, even though flux is running through the system. Second, from a mathematical point of view, many analyses at a steady state are much simpler than for the differential equations themselves, because now one has, by definition, Ẋ = 0 , so that all differential equations become explicit algebraic equations that can be analyzed with methods of linear algebra. If all fluxes are known, it is usually not difficult to compute the steady-state of a system. However, the reverse is not true: if only the metabolite concentrations at the steady state are known, it is not easy to compute the corresponding flux distribution, because metabolic systems almost always contain more reactions than variables. In this case, optimization methods like FBA or MOMA need to be employed.
In the Brachypodium study, we chose an alternative to FBA and MOMA. Namely, we intended to obtain the most likely solution without specifying an objective function for the FBA optimization. Because the degrees of freedom of a solution to our system are directly associated with diverging branch points, we focused on the flux split

Table 1 Computational model results compared to experimental data.The model results demonstrate a good match with experimental data in terms of total lignin (A) and the incorporation of label (B)
a Label incorporation in H-lignin was not considered as a criterion during the model calibration. The recorded experimental value in the [U-13 C 9 ]phenylalanine feeding experiment is greater than the reported label level in phe, which is 35% [29]. As a consequence, we deemed the measurement unreliable and did not use  ratios (FSRs) at these points. In cases where these FSRs were known, we used their values; otherwise, we performed large-scale Monte-Carlo simulations with thousands of combinations of FSRs and retained only those solutions where all fluxes were positive at all time points of an experiment. This strategy led to the most likely flux profiles. Details of this method are discussed in [30]. As a simplified example, consider the hypothetical pathway in Fig. 9, which has two FSRs, and hence two degrees of freedom.
The system of differential equations corresponding to pathway in Fig. 9 is Given a set of metabolite concentrations over time, the pathway can be driven by infinitely many flux distributions [31]. To determine the most likely, the system in Eq. (3) It is theoretically possible that the estimation strategy based solely on split ratios does not converge to an acceptable solution, and we have discussed means of addressing this situation elsewhere [30]. Here the splitratio method succeeded without the need for alternative methods.

Modeling 13 C-labeling experiments
The diffusion flux between two pools of the same metabolite in different locations is comprised of two directions (Fig. 10). Although the two opposing fluxes have a net value, as shown in Figs. 6 and 9, it is necessary to consider them individually when modeling a labeling experiment. The reason is that the labeling content of each pool affects the flow of label, but the net diffusion alone would not reflect the free passing of label in both directions. For instance, the illustration scheme in Fig. 9 does not allow flow of label from X 4 to X 2 through D when labeled metabolite is fed to the pathway through X 3 , but due to the bidirectional nature of diffusion fluxes, it is evident that flow would happen in reality. As Fig. 10 illustrates, these bidirectional diffusion fluxes form cycles and don't allow the direct computation of steady-state fluxes.
To tackle this issue, we first consider only the net diffusion flux as shown in Fig. 9 and compute the steady states. Then we consider the bidirectional model in Fig. 10, employ Eq. (2), and use conservation of mass for labeled and total fluxes at each metabolite. For a given labeling percentage L i , where L i represents the labeled portion of the pool of metabolite X i , we obtain for the pathway system in Fig. 10: which can be rewritten as Fig. 9 Material flow through net fluxes in an illustration example. Without labeling, it is sufficient to model diffusion fluxes as net fluxes. However, this is not the case for labeling experiments (Fig. 10)   Fig. 10 Illustration of the flow of label in the same example as Fig. 9, but with explicit flux directions. In contrast to the scenario in Fig. 9, labeling experiments mandate the modeling of diffusion fluxes in both directions. Specifically, the simpler model in Fig. 9 does not allow flow of label from X 4 to X 2 through D when labeled metabolite is fed to the pathway through X 3 , but the figure here demonstrates that such flow is clearly possible Similar calculations for pool X 4 yield By equating D f from Eqs. (6) and (7) one obtains Assuming that L 1 and L 3 are known from inputs of the pathway, V in , we can compute L 4 from the computed steady-state fluxes in the previous step and an estimated L 2 . Therefore, L 2 is the only unknown to be estimated, and L 4 , D f and D r can consequently be computed. In fact, we only need to estimate the label level of one of the parallel metabolite pools in the cytosol and ER compartments.
Similar to the use of split ratios, vector L is generated randomly by Monte-Carlo sampling, and the labeled fluxes can then be computed. The labeled fluxes corresponding to Fig. 10 are Closer inspection demonstrates that the model in Fig. 9, which considers only net diffusion, computes the labeled portion of D as L 2 D, which is equal to L 2 (D f − D r ), whereas Eq. (9) quantifies the net labeled flux as The fluxes for the Brachypodium example were defined in this manner. Labeled fluxes that satisfied the model criteria for labeling experiments were recorded. The criteria were • The labeled lignin composition is compatible with 13 C 9 -phenylalanine and 13 C 9 -tyrosine experimental data; and • The labeled ferulic acid and p-coumaric acid levels match the experimental data.
The recorded flux vectors were plotted using boxplots, which offer a visual representation of the distribution of most likely flux values within their admissible ranges.

Discussion
Mathematical modeling in biology is still in its infancy.
Especially within the realm of plant and crop science, the number of modeling articles is negligible in comparison to experimental papers. As a consequence, the collective experience with plant and crop modeling approaches is still limited, and much more practice and many more case studies will be needed to gain a glimpse into the systemic responses of plants to interventions and manipulations. It may even be, as some experts claim (Leroy Hood, pers. comm.), that a "new math" is needed that allows us to combine different data and heterogeneous information in a more efficacious manner than is possible today. Ultimately, a deeper understanding of such responses would allow us to answer questions like "how does 'a' plant react to natural or artificial changes?" or "why does plant (or plant species) A respond differently to a perturbation than plant (or plant species) B?" To obtain more practice and experience of this type, experimentalists and modelers should collaborate more closely. On the one hand, modelers will need experiments specifically performed for some modeling aspects. At present, many data are available, and the data flow from -omics experiments can be overwhelming. However, not all data are useful for the type of modeling outlined in this article, and modelers will be dependent on experimentalists to perform other types of experiments [32]. On the other hand, experimentalists will want to see genuinely new results coming out of models, especially if they had contributed data to the modeling effort. They will benefit from new, integrative interpretations of their data and from reliable modeling results and computationally achieved hypotheses guiding the "next steps" in their research programs. The generic differences between laboratory or field experiments and computational approaches render it evident that this type of collaboration has true potential, but that it will take time and patience on both sides to make progress toward reaching some of this potential.
As a tangible target, experimentalists and modelers should explore together to what degree metabolic responses can be predicted (qualitatively or quantitatively) from the existence of genes and enzymes (as, for instance, TAL in the case of Brachypodium) or from quantitative transcriptomics, where one would expect to find similarities between gene expression and changes in enzyme activities, which however do not always materialize in reality, due to post-transcriptional alterations. It would also benefit both sides to obtain and computationally analyze data describing the same process in different species, as we demonstrated here with the different models for lignin biosynthesis. At first, these comparative analyses could shed light on questions such as: whether apparent differences between species are experimental or modeling errors; whether different designs have evolved in line with the general phylogeny of these species or whether they are due to other factors; and whether different natural designs are dictated by different environmental needs or demands. Together, these combined analyses would have the potential of revealing design principles that govern these processes and could provide deep explanations for why certain species solve a task in the observed fashion and not in a different fashion.
As an example, we discussed lignin biosynthesis. The lignin heteropolymer and its monomers are quite similar among different species and, indeed, the intermediates of their biosynthetic pathway are essentially the same. However, the enzymatic reactions of the pathway exhibit striking differences, not just in terms of their kinetic features or parameters, but even in their existence. These differences raise the question why one pathway design is favored in one case, but not in another. The answer to this question is not only academically interesting but is of immediate pertinence to the biofuel scientist and metabolic engineering in general, because ignorance of the true reasons mandates a new conceptual model for every untested species. By contrast, knowledge of a general design or operating principle would allow predictions regarding the pathway topology of a new species based on the criteria on which the principle is founded.
The different pathway models for lignin synthesis in a number of plant species [8,9,16,25,27,28] have revealed commonality, but also differences in regulatory features and, as suggested here, compartmentalization. One could thus come to the conclusion that every species manages its lignin production differently. However, the fact that some distinctive features of one model are not part of the other models should not be over-interpreted, at least not quite yet. It is well possible that the pathway in alfalfa and switchgrass is as compartmentalized as the one in Brachypodium, but the dictum of simplicity in modeling, and the data that these models were based upon, suggested that specific compartments were not needed to match the data in these species. The same is true for other apparently distinguishing features, such as product inhibition, which we found necessary in switchgrass, and which may well be in effect in other species. These features were needed to make the models consistent with specific data typesets, and if one re-analyzed the models with other types of data types, the same features could well be suggested for other species. As it stands, the collective experimental database is sparse, and the published models are minimalistic special cases of the same "master model, " which can even account for the fact that some species seem to be missing certain pathway metabolites or enzymatic reactions. It remains to be the subject of further data generation and analysis to determine whether these differences disappear toward one common, complex model, whether they are immaterial byproducts of evolution that did not exert strong selective pressure, or whether they evolved for reasons that are germane to these species and their environments. The cooperation between experimentalists and modelers has led to early successes. Some of these are narrowly focused by explaining observations that had been puzzling before. We described some of these in the context of lignin synthesis and recalcitrance. To study some of the in vivo complexity in an in vitro system, the lignifying cell suspension cultures reported in several species (Arabidopsis [33], poplar [34], switchgrass [35]) could be useful for modeling purposes of the lignin pathway. These systems can be studied along a time course when lignin deposition and cell differentiation occur, allowing the evaluation of different parameters such as pH or temperature and the use of dynamic models that could be proposed as potential in vivo validation systems.
Others studies have attempted to connect several scales of the biological hierarchy of processes, both in time and size (e.g., [36][37][38][39][40][41]). Not surprisingly, such much larger models cannot account for every detail at every lower level. Nonetheless, it might be useful, for instance, to use agent-based models at the highest scale considered, such as overall plant growth, and to anchor into them detailed models of key sub-systems, such as photosynthesis, respiration, and stress responses. Instead of an agent-based model, the highest level could also be a dynamic FBA model [42], or it is even possible to use low-level models as constraints in genome-wide metabolic models [43].
The epitome of such models is whole-plant-plus-environment models that have been developed in recent years and capture governing processes and responses quite well. For instance, SOYSIM [44] is a computational model that simulates soybean growth on a daily basis throughout its lifecycle. It permits reasonably accurate explorations of water use, additional irrigation, and potential yield under different conditions. Another example is the WIMOVAC simulation model, which allows investigations of the carbon balance in plants and permits predictions of crop responses to changes in climate [45][46][47][48]. It is applicable to different plant and soil types and can be used by researchers, managers and students as an exploratory tool. These organismal models are now to be coupled more comprehensively to environmental and agricultural models, an effort that has recently been coined Crops in silico [49][50][51]. To be successful, this effort will require the close collaboration not only between experimental plant scientists and mathematical modelers, but also involve experts in biophysics, hydrogeochemistry, meteorology, high-performance computing, visualization, and many other fields. The challenges are grand indeed, but a solid foundation is being built by the collaboration of several communities, such as BESC [52], and it will only be a matter of time and collective willpower to increase momentum allowing us to achieve some of the set goals.


extfieldsofstudy: Chemistry Medicine

text:
Prohibitin 2 localizes in nucleolus to regulate ribosomal RNA transcription and facilitate cell proliferation in RD cells

Prohibitin 2 (PHB2), as a conserved multifunctional protein, is traditionally localized in the mitochondrial inner membrane and essential for maintenance of mitochondrial function. Here, we investigated the role of PHB2 in human rhabdomyosarcoma (RMS) RD cells and found substantial localization of PHB2 in the nucleolus. We demonstrated that PHB2 knockdown inhibited RD cell proliferation through inducing cell cycle arrest and suppressing DNA synthesis. Meanwhile, down-regulation of PHB2 also induced apoptosis and promoted differentiation in fractions of RD cells. In addition, PHB2 silencing led to altered nucleolar morphology, as observed by transmission electron microscopy, and impaired nucleolar function, as evidenced by down-regulation of 45S and 18S ribosomal RNA synthesis. Consistently, upon PHB2 knockdown, occupancy of c-Myc at the ribosomal DNA (rDNA) promoter was attenuated, while more myoblast determination protein 1 (MyoD) molecules bound to the rDNA promoter. In conclusion, our findings suggest that nucleolar PHB2 is involved in maintaining nucleolar morphology and function in RD cells by regulating a variety of transcription factors, which is likely to be one of the underlying mechanisms by which PHB2 promotes tumor proliferation and represses differentiation. Our study provides new insight into the pathogenesis of RMS and novel characterizations of the highly conserved PHB2 protein.

Rhabdomyosarcoma (RMS) is the most prevalent soft tissue sarcoma in children and adolescents, accounting for 5% of all pediatric tumors and over 40% of pediatric soft tissue sarcomas 1 . RMS can be grouped into two major histological subtypes, embryonal (ERMS) and alveolar (ARMS) rhabdomyosarcoma 2 . ERMS represents approximately 70% of all childhood RMS, mainly affecting the head and neck and genitourinary regions 3 . Advances in the treatment of RMS have promoted the 5-year survival rate from 25% up to approximately 70% 4 . However, a high rate of metastases, drug resistance and tumor recurrence remain to be overcome. The lack of directed therapies against RMS emphasizes the need to further illuminate the molecular underpinnings of the disease.
Myogenic differentiation arrest is a main characteristic in RMS 5 . Accompanied by incomplete differentiation, cell proliferation is no longer arrested in these tumor cells 2 . The skeletal muscle differentiation process is conserved in evolution and involves an orderly cascade of events dependent on the activities of two families of transcription factors, myogenic regulatory factors (MRFs) and the myocyte enhancer factor-2 (MEF2) 6 . MRFs, which include MyoD, myogenin, myogenic factor 5 (Myf5) and MRF4, synergistically cooperate with MEF2 to induce muscle-specific gene transcription and ultimately the onset of myogenesis. Among those factors, MyoD is considered to be muscle-determining, and the expression of myogenin is a typical early myogenic differentiation marker. Although RMS cells express MyoD and myogenin, they fail to achieve complete myogenesis under differentiation conditions via mechanisms that are not fully understood. Previous studies have associated higher expression of myogenic markers with better prognosis in children with RMS 7 . Therefore, gaining a better understanding of pathogenesis that contribute to the proliferation and growth of RMS is critical.
PHB2, also designated as B-cell receptor-associated protein 37 (BAP37) and repressor of estrogen receptor activity (REA), is a highly conserved protein found in fungi, plants and humans 8,9 . Mainly, PHB2, together with PHB1, forms the eukaryotic mitochondrial prohibitin complex which is essential for mitochondrial morphogenesis and genome stability 10,11 . In addition, PHB2 has been implicated in diverse mitochondria-related functions, such as cell proliferation, cell death and aging 12 . However, its diverse functions in other cellular compartments have also been widely reported, such as in the nucleus and on the membrane 13,14 . We previously demonstrated that PHB2, as a repressor, could inhibit muscle differentiation by recruiting histone deacetylase 1 (HDAC1) to repress the transcriptional activity of both MyoD and MEF2 15 . However, whether PHB2 is involved in the development and progression of RMS has not been reported.
Our current study aimed to explore the role of PHB2 in RMS by using the RNA interference (RNAi) technique. We found that PHB2 knockdown inhibited the proliferation of the (human RMS) RD cell line, as evidenced by arrest of cell cycle and reduction of DNA synthesis. Meanwhile, PHB2 knockdown also induced apoptosis in a fraction of the cells and enhanced the differentiation of RD cells to some extent. More importantly, we observed that a substantial amount of PHB2 was localized in the nucleolus which has never been reported before. We further showed that nucleolar PHB2 might be involved in maintaining the structure and function of the nucleolus in RD cells partly by modulating the transcription of rDNA through regulating the accessibility of c-Myc and MyoD to the rDNA promoter. In summary, our results suggest a pro-tumorigenic role for PHB2 in human RMS at least partly through its specific localization in the nucleolus.

PHB2 knockdown by siRNA inhibits RD cell proliferation.
To explore the functional role of PHB2 in RMS, we first knocked down endogenous PHB2 in RD cells, an ERMS cell line, using siRNA. As shown in Fig. 1a, two separate siRNAs targeting PHB2 could both effectively down-regulate the expression of PHB2 of RD cells. Thereafter, we determined the effect of PHB2 siRNAs on RD cell proliferation in a time course experiment. Results of the MTT assay showed that RD cells transfected with PHB2 siRNAs displayed a reduced proliferation rate compared with cells transfected with the negative control siRNA (Fig. 1b). We further demonstrated that the expression of antigen Ki67, a nuclear protein associated with cell proliferation, also substantially decreased in RD cells transfected with PHB2 siRNAs (Fig. 1c). As DNA synthesis is directly related to the ability of cells to proliferate, we also examined DNA synthesis in siRNA-transfected RD cells by using a BrdU incorporation assay. Similarly, BrdU incorporation was decreased in RD cells treated with PHB2 siRNAs compared with control cells (Fig. 1d). Next, we investigated whether PHB2 silencing would also impact the clonogenic ability of RD cells. The results showed that PHB2 knockdown significantly reduced the colony forming efficiency of RD cells (Fig. 1e). Collectively, these findings indicated that down-regulation of PHB2 would inhibit the proliferation of RD cells, suggesting that PHB2 is required for RD cell proliferation. Interestingly, knockdown of PHB2 expression in normal rat myoblast L6 cells and normal liver L02 cells did not significantly affect the cell proliferation (see Supplementary Figs S3 and S4).
PHB2 knockdown results in cell cycle arrest and apoptosis in RD cells. The cell cycle plays a key role in cell proliferation. Thus, we further analyzed the cell cycle of siRNA-treated RD cells by using flow cytometry in order to evaluate whether the inhibition of RD cell proliferation by PHB2 siRNAs was achieved by arresting the cell cycle. As shown in Fig. 2a, a higher proportion of cells transfected with PHB2 siRNAs was arrested in the G1 phase compared with negative control cells. Upon transfection with siPHB2-1 and siPHB2-2, the percentage of G1 phase cells increased to 50.90% and 52.01%, respectively (Fig. 2a). To investigate the underlying molecular mechanisms, we examined various cyclins which are key cell cycle regulators. The results showed that PHB2 knockdown reduced the expression of cyclin E (Fig. 2b), while it had no obvious effect on the levels of cyclin A and cyclin B1. Cyclin E plays an essential role in the G1-S transition. Therefore, these results validated that PHB2 silencing could indeed induce cell cycle arrest at the G1 phase. Next, we tested the possibility of inducing apoptosis by PHB2 siRNAs. We found that although apoptosis was not induced in a large number of RD cells, transfection with PHB2 siRNAs resulted in a slight but stable increase in the proportion of those cells undergoing apoptosis compared with negative control cells (3.7% vs 8.5% or 12.5%) (Fig. 2c). Thus, the decreased expression of PHB2 inhibited the proliferation of RD cells possibly by impairing cell cycle and inducing cell apoptosis.
PHB2 knockdown promotes RD cell differentiation. The clinical outcome of RMS is associated with the degree of muscle differentiation, and greater expression of myogenic transcriptional factors indicates a better prognosis. We wondered whether PHB2, as an inhibitor of myogenic transcription factor and muscle differentiation, is linked to differentiation defects in RMS. Therefore, we investigated the effect of PHB2 knockdown on the differentiation of RD cells. Myogenin is a transcriptional factor critical for differentiation of skeletal muscle and also designated as an early myogenic differentiation marker. The lowered expression of PHB2 was observed to lead to the higher expression of myogenin to some extent (Fig. 3a). We further detected myogenin by immunofluorescence staining in RD cells. Consistent with the immunoblotting results, the percentage of myogenin-positive cells increased in PHB2-silenced RD cells, compared with the negative control ( Fig. 3b and c). Altogether, these results indicated that PHB2 might participate in maintaining the low degree of differentiation in RD cells.

PHB2 is partially localized in the nucleolus in RD cells.
To further delineate the biological function of PHB2 in RD cells, we evaluated the subcellular localization of PHB2 in RD cells. By using immunofluorescence, we observed that aside from the expected mitochondrial localization, a substantial amount of PHB2 seemed to be localized in the nucleolus of RD cells. To further confirm this new observation, we examined the cellular ScIENtIfIc RepoRts | (2018) 8:1479 | DOI:10.1038/s41598-018-19917-7 co-localization of PHB2 and nucleophosmin (NPM, also known as B23, a nucleolar protein) by immunofluorescence staining and confocal microscopy. PHB2 showed an obvious co-localization with NPM in RD cells. NPM is a phosphoprotein localized in granular regions of the nucleolus, which is involved in diverse biological functions including genomic stability, tumorigenesis and ribosome biogenesis [16][17][18] . This result confirmed that PHB2 is located in the nucleolus of RD cells (Fig. 4a). Consistently, partial nucleolar localization of PHB2 was also observed in two other ERMS cell lines, A204 and TE671, and one ARMS cell line RH30 ( Supplementary  Fig. S5). Interestingly, in normal cells, such as L6 and L02 cells, PHB2 was mainly localized in the cytosol/mitochondria and hardly observed in the nucleolus (Fig. 4b and c). We then extended the co-localization analysis to other human cancer cell lines, such as HepG2, MDA-MB-231 (MM231) and A2780. The nucleolar localization of PHB2 was also observed in HepG2 cells (Fig. 4d), but not in MM231 and A2780 cells ( Fig. 4e and f). These data indicated that PHB2 might elicit unique regulatory functions in RD cells via its specific subcellular localization to the nucleolus.
PHB2 silencing alters morphology of the nucleolus in RD cells. Since a considerable amount of PHB2 appeared in the nucleolus of RD cells, we speculated that it may be involved in regulation of nucleolar morphology and function. To test this hypothesis, we observed the subcellular structures of siRNA-transfected RD cells by using transmission electron microscopy (TEM). As expected, a series of alterations could be found in the morphology and architecture of the nucleolus in PHB2 knockdown RD cells compared with negative control cells. The nucleolus in mammalian cells has three primary structures with morphologically and functionally distinct components: fibrillar centers (FCs), dense fibrillar component (DFC) and granular component (GC) 19 . FCs  RD cells were transfected with siRNA and the growth medium was replaced 6 hours later. The next day, the medium was replaced with the differentiation medium (2% calf serum, H-DMEM). Cell extracts were subjected to immunoblotting with indicated antibodies, and GAPDH was detected as a loading control. Full-length blots are presented in Supplementary Fig. S1. (b) Immunofluorescence staining of myogenin in RD cells. Forty-eight hours after transfection, RD cells were fixed and subjected to immunostaining for myogenin (red). The nuclei of cells were counterstained with DAPI (blue). Scale bar 200 μm. Results are representative of three independent experiments. (c) Quantitative analysis of myogenin in panel B. Data are presented as the mean ± S.D. (n = 3), *P < 0.05, **P < 0.01, ***P < 0.001, two-tailed Student's t test.
house the bulk of the rDNA. The synthesis of pre-ribosomal RNA (rRNA) executed by RNA polymerase takes place at the interface between FCs and DFCs 20,21 . Early pre-rRNA cleavages occur in DFC territory 22 , and further processing is at the GC 23 . In our study, we observed that PHB2 knockdown led to a series of alterations in the nucleolus, i.e. increased nucleolar density, disappearance of sponge-like nucleolus, irregular nucleoli structure, fuzzy nucleolus boundaries and nucleolus contraction (Fig. 5a). These changes could partly explain the inhibitory effect of PHB2 silencing on RD cell proliferation. Aside from these alterations in the nucleolus of RD cells, PHB2 knockdown also led to the damage of mitochondria, as previously reported, such as loss of structural integrity of mitochondria, and cristae becoming intumescent and intricately stacked (Fig. 5b). In addition, we also found that PHB2 silencing induced a more uniform cell boundary and fewer cell microvilli (Fig. 5c).
Collectively, these data indicate that PHB2 deficiency by siRNA treatment in RD cells impaired the ultrastructure of the nucleolus and mitochondria, likely contributing to the inhibition of RD cell proliferation and subsequent inhibition of the tumorigenicity of these cancer cells.
PHB2 silencing suppresses rRNA transcription, reduces occupancy of c-Myc at rDNA and facilitates MyoD binding to rDNA. rRNA is the RNA component of the ribosome, participating in ribosomal biogenesis. Transcription and processing of most rRNAs occur in the nucleolus. To investigate whether PHB2 is involved in regulating nucleolar function, we examined the effect of PHB2 knockdown on transcription of rDNA in the nucleolus. We detected the level of 45S rRNA, a pre-rRNA transcribed by RNA polymerase I, and its mature 18S rRNA product by quantitative real-time RT-PCR (RT-qPCR) in siRNA-treated RD cells. The results showed that both 45S pre-rRNA and 18S rRNA were down-regulated in siPHB2-transfected RD cells compared with control cells (Fig. 6a). These data suggest that rRNA gene transcription was suppressed in siPHB2-transfected RD cells.
The above results indicated that nucleolus-localized PHB2 might participate in regulation of the rRNA transcription in RD cells. To further clarify the possible underlying mechanisms, we investigated the effect of PHB2 knockdown on the binding of c-Myc to the rDNA promoter. The c-Myc oncogene is one of the well-known transcription factors that activates rDNA transcription in response to mitogenic signals during proliferation 24,25 . The rDNA repeats contain numerous canonical and non-canonical E-boxes [element recognized by  24 . We therefore used several sets of primer pairs 25 to perform a chromatin cross-linking and immunoprecipitation (ChIP) assay. As shown in Fig. 6b, PHB2 silencing decreased the enrichment of rDNA immunoprecipitated with anti-c-Myc antibodies, which may have contributed to the reduced rRNA transcription. In addition, our previous study showed that PHB2 could interact with the bHLH transcription factor MyoD and repress MyoD-dependent gene transcription in myoblasts 15 . However, aside from catalyzing the RNA polymerase II-mediated muscle-specific gene transcription, Ali et al. also demonstrated that MyoD could occupy rDNA loci and suppress rRNA expression during skeletal muscle lineage progression, concomitant with decreased rRNA expression and reciprocal loss of occupancy by c-Myc 26 . Therefore, we further investigated the occupancy of MyoD on rDNA in siRNA-transfected RD cells. Results of the ChIP assay with an anti-MyoD antibody showed that the occupancy of MyoD on rDNA repeats exhibited a significant increase when PHB2 was silenced (Fig. 6c). Taken together, our data indicate a potential PHB2-dependent mechanism which inhibits myogenic differentiation and facilitates proliferation through promoting c-Myc activation of rDNA transcription, while also inhibiting MyoD binding to rDNA at least in RD cells.

Discussion
RMS is the most common soft-tissue sarcoma in children and adolescents, which is a mesenchymal malignancy arising as a consequence of regulatory disruption of skeletal muscle progenitor cell growth and differentiation. As a result, RMS cells are unable to undergo full-scale terminal myogenic differentiation despite the expression of MyoD and myogenin, two important factors required for myogenic cell differentiation [27][28][29] . PHB2, a highly evolutionarily conserved protein, has been demonstrated to possess multiple biological activities, including transcriptional regulation, nuclear signaling, mitochondrial structural integrity maintenance and cell division modulation in various cellular compartments 30 . Here, using the RD cell line, we first showed the localization of PHB2, a repressor of skeletal muscle differentiation, in the nucleolus where it participated in regulation of Myc-mediated rRNA synthesis, which may serve as a mechanism to facilitate cell proliferation while inhibiting differentiation.
The involvement of PHB2 in the regulation of cancer cell proliferation has been reported. However, its effects are debated in cancer cells of different origins. In estrogen receptor (ER)-positive breast cancer, PHB2 acts as a tumor inhibitor to suppress the growth of breast cancer cells by repressing the estrogen/ER signaling pathway 31 . Meanwhile, in hepatocellular carcinoma, PHB2 is required for cancer cell proliferation 32 . In the present study, we found that silencing of PHB2 could obviously inhibit RD cell proliferation and survival (Fig. 1). Furthermore, we observed that upon PHB2 silencing, RD cells showed a stable tendency of cell cycle G1 phase arrest. These results are consistent with the study by Equilibrina et al., which demonstrated that PHB2 was required for G1/S progression since the number of cells undergoing extended G1 phase increased significantly under the PHB2 knockdown condition 33 . We also found that PHB2 silencing could induce apoptosis in a small proportion of RD cells. This result is in agreement with a previous report showing that PHB2 knockdown in Hela cells resulted in apoptosis 34 . However, we did not find that PHB2 silencing could affect the migration and invasion of RD cells (data not shown), which disagrees with previous findings 13 . Therefore, the precise activities of PHB2 in these cellular processes are worth further study.
Aside from supporting RD cell proliferation, PHB2 might also play a role in the differentiation state of these cells. Our results showed that PHB2 silencing resulted in a higher level of myogenin expression in RD cells (Fig. 3) and an increased proportion of myogenin-positive RD cells. These results are expected since we previously demonstrated that PHB2 could act as a repressor for muscle differentiation by inhibiting the transcriptional activities of MyoD and MEF2 15 . Although RMS cells express multiple essential myogenic transcription factors such as MyoD, MEF2 and myogenin, they cannot fully differentiate 35 . The uncontrolled repressive activity of PHB2 may partly contribute to the defective differentiation state of RMS cells. In other tumors, the role of PHB2 in maintaining the poorly differentiated state has not been reported thus far. Notably, the differentiation enhancement by PHB2 silencing on RD cells was not as significant as expected. We speculate that the main role of PHB2 in RD cells is to regulate cell proliferation, although we cannot exclude the possibility that the results were affected by low siRNA transfection efficiency.
PHB2 is located in different cellular compartments to elicit different functions 36 . In inner mitochondrial membranes, PHB2 (by interacting with PHB1) maintains the structure of mitochondria and the integrity of mitochondrial respiration machinery 37,38 . In the nucleus, PHB2 is involved in transcription regulation by interacting with multiple transcription factors, including estrogen receptor, MEF2, MyoD 15,34,39 . In the plasma membrane, PHB2, as a receptor-associated protein, is required for insulin-like growth factor binding protein (IGFBP)-6 induced cell migration of RMS cells 13 . PHB2 in the mitochondria is found dominantly and ubiquitously in almost in all tested cells, while other subcellular localizations may be dependent on the type or given stage of cells and are usually restricted to only a small fraction of PHB2 molecules. Here, we identified a new pattern of subcellular localization in RD cells by observing quite a few PHB2 molecules in the nucleolus (Fig. 4). Meanwhile, we further found that PHB2 also localized in the nucleolus in some other examined tumor cells, but a clear nucleolar localization of PHB2 was not observed in normal myoblasts L6 cells, normal liver L02 cells and several types of tumor cells (MM231 and A2780). Thus, mechanisms that mediate PHB2 localization in the nucleolus or other cellular compartments are worth deeper investigation.
The nucleolus, as the most prominent nuclear substructure, is the center of rDNA transcription and facilitates ribosome synthesis. rDNA transcription is important for malignant cell growth and proliferation because the cell must increase ribosome biogenesis to meet the high demand of protein synthesis. Increased sizes and numbers of nucleoli have been recognized as hallmark features of most of tumor types 40 . Recently, increased ribosome biogenesis, which is an important event of tumorigenesis and cancer progression 41 accompanied by an increased rate of rRNA synthesis in the nucleolus, has become a hallmark of cellular proliferation 42 . Therefore, aberrant nucleolar substructure and function are closely related to cancer progression. Our data showed that PHB2 silencing led to smaller and more dense nucleoli compared with the negative control (Fig. 5). Accordingly, upon PHB2 knockdown, the transcription of rDNA was repressed, and the rRNA level clearly declined (Fig. 6). These results suggest that PHB2 directly or indirectly participates in the regulation of the nucleolus function, which may be one of the mechanisms by which PHB2 is involved in regulating cell proliferation.
Myc is an important and direct regulator of rRNA transcription 43 . In this study, we found that PHB2 silencing in RD cells significantly decreased Myc occupancy at the rDNA promoter, which may account for the reduction in rRNA level upon PHB2 knockdown by siRNA. This data indicates that nucleolar PHB2 is required for Myc to regulate the transcription of rRNA at least in RD cells. On the other hand, MyoD, as a phenotypic transcription factor, has been shown to control differentiation of mesenchymal progenitors into myoblasts by suppressing rRNA gene transcription through selectively interacting with rRNA genes 26 . Our previous studies have demonstrated that PHB2 is a binding partner and repressor of MyoD and inhibits the transcriptional activity of MyoD for skeletal muscle specific structural gene expression during myogenesis 15 . Interestingly, here we also found that down-regulation of PHB2 in RD cells enhanced the binding of MyoD to the promoter of rDNA (Fig. 6), suggesting that PHB2 not only represses the transcriptional activity of MyoD but also alleviates the repressive effect of MyoD on rRNA transcription. Therefore, the present data indicate that the repression of MyoD binding to rDNA by nucleolar PHB2 is one of the mechanisms underlying the differentiation defect of RD cells even though MyoD is expressed. Altogether, these results illustrate the functional significance of the nucleolar localization of PHB2 in RD cells. In view of the above results, we hypothesize that PHB2 plays multiple pathophysiological roles in rhabdomyosarcoma. On one hand, it supports proliferation by enhancing the activity of Myc in nucleolus and/ or maintaining normal mitochondrial processes. On the other hand, PHB2 inhibits differentiation by repressing the activity of MyoD.
Cellular localization is closely related to their biological functions of proteins. Both mitochondria and nucleolus localized PHB2 can possibly play a role in regulating RD cell proliferation, survival and differentiation. However, our present study did not provide direct evidence to demonstrate the localization-function correlation of PHB2 in RD cells since we could only silence global PHB2 expression and not specifically down-regulate either mitochondrial or nucleolar PHB2. Although a considerable amount of PHB2 does exist in nucleolus of RD cells and the morphology and function of nucleolus were indeed destroyed by PHB2 silencing in RD cells, we could not exclude the possibility that alterations in the PHB2 level of the mitochondria may indirectly influence rDNA expression and nuclear stability. Therefore, future studies are needed to clarify how PHB2 localizes in different cellular compartments and its functional roles in those different compartments. Notably, some of our preliminary data showed diverse regulatory roles of PHB2 in various cell lines. Upon PHB2 silencing, the human embryo liver cell line L02 neither showed a decrease in cell proliferation, nor exhibited cell cycle G1 phase arrest and reduced rDNA expression (Supplementary Fig. S3). Similarly, the proliferation of another normal cell line, rat myoblast L6 cells, was also not affected by PHB2 silencing (Supplementary Fig. S4). Both of these normal cell lines did not show obvious nucleolar PHB2 staining in our study. Meanwhile, cell proliferation in human hepatocarcinoma HepG2 cells and RH30 cells, which are positive for nucleolar PHB2 staining, was significantly inhibited by PHB2 down-regulation to different extents ( Supplementary Figs S2a and S5d). However, compared to RD cells, HepG2 cells transfected with PHB2 siRNAs demonstrated no obvious cell cycle arrest, and rDNA expression was only decreased slightly ( Supplementary Fig. S2b and c). These preliminary data point to the complicated working mechanisms of PHB2 in cellular processes, which may differ between normal and malignant cells or vary among cells with different tissue origins or in distinct stages. All of these phenomena are worth further investigation.
In conclusion, for the first time our study identified a unique localization for PHB2 in RD and other tumor cells, which further deepens our understanding of this conserved protein and highlights the significance of its function. PHB2 may play a different biological role in different cells or different physiological and pathological conditions. Novel and complex roles of PHB2 are continually brought to light. Wei et al. recently found that PHB2 is a crucial mitophagy receptor involved in targeting mitochondria for autophagic degradation 44 . In our study, ectopic nucleolar-localized PHB2 in RMS, as a regulator of rDNA transcription, was found to support Myc-dependent rDNA transcription and maintain nucleolar structure which participates in sustaining the survival and proliferation of RD cells. In addition, PHB2 via the regulation of the access of MyoD to rDNA adopts an additional mechanism to repress MyoD activities and thus keep RD cells less differentiated. Altogether, the findings of our study provide new insight into the pathogenesis of RMS. RNA extraction and RT-qPCR. Total RNA was prepared from cultured cells using Trizol reagent (Invitrogen) following the manufacturer's instructions. Total RNA was reverse transcribed into cDNA using a real-time PCR kit (TransGen Biotech, Beijing, China). RT-qPCR was performed using 2 × FastStart Universal SYBR Green Master (Rox) (Roche Diagnostics, Mannheim, Germany) according to the manufacturer's protocol. Cycling conditions were 95 °C for 10 min, followed by 40 cycles of 95 °C for 10 s and 60 °C for 1 min. Sequences for primers used for RT-qPCR are as follows: 18S rRNA (sense: 5′-GCAATTATTCCCCATGAACG-3′, antisense: 5′-GGGACTTAATCAACGCAAGC-3′) and 45S rRNA (sense: 5′-TGTCAGGCGTTCTCGTCTC-3′, antisense: 5′-GAGAGCACGACGTCACCAC-3′). Relative quantification of the expression of each gene was calculated using the △△CT method. All real-time PCR reactions were performed in triplicate.
BrdU incorporation assay. RD cells were incubated with 10 µl/well of 10 µM BrdU labeling solution for 6 h at 37 °C. The cells were then fixed with 200 µl FixDenat for 20 min and incubated with 100 µl/well anti-BrdU-POD working solution for 90 min. Subsequently, substrate solution (100 µl/well) was added and incubated for 5-30 min. BrdU uptake was measured at 450 nm using a microplate reader (Bio-Rad, Hercules, CA, USA). Immunofluorescence staining. Cells were first fixed in 4% paraformaldehyde for 10 min and then permeabilized with 0.2% TritonX-100 in PBS for 10 min before being blocked with PBS containing 5% bovine serum albumin (BSA) for 30 min. Thereafter, the cells were sequentially incubated with anti-PHB2 and anti-NPM primary antibodies and then with anti-rabbit IgG (H + L) (Alexa Fluor ® 647 Conjugate, Cat # 4414, Cell Signaling Technology) and anti-mouse IgG (H + L) (Alexa Fluor ® 488 Conjugate, Cat # 4408, Cell Signaling Technology), respectively. The images were acquired using an Olympus FV-1000 confocal laser scanning microscope (Olympus, Tokyo, Japan).

Flow cytometry analysis.
MTT assay. RD cells transfected with siRNA were seeded in 96-well plates at an initial cell density of 5000 cells per well. After 24 h (defined as the first day), 20 µl of MTT solution (5 mg/ml) was added to each well and incubated for an additional 4 h. The absorbance was measured at 570 nm using a microplate spectrophotometer (Thermo Electron Corporation, Shanghai, China). The experiment was performed on three replicates and repeated three times.
Clonogenic assay. RD cells transfected with siRNA were plated in 6-well plates at 1500 cells per well in 2 ml of media. After 3 days, cells were washed with fresh media and cultured in 2 ml of fresh media. Ten days later, colonies were fixed in 80% ethanol, stained with 0.1% crystal violet and counted.
Electron microscopy analysis. RD cells were transfected with indicated siRNAs. Forty-eight hours later, cells were harvested and washed with PBS. The cells were then processed as previously described 47 . Ultrathin sections were examined with the FEI TECNAI SPIRIT electron microscope (Brno, The Czech Republic).