Computational Complementation: A Modelling Approach to Study Signalling Mechanisms during Legume Autoregulation of Nodulation

Autoregulation of nodulation (AON) is a long-distance signalling regulatory system maintaining the balance of symbiotic nodulation in legume plants. However, the intricacy of internal signalling and absence of flux and biochemical data, are a bottleneck for investigation of AON. To address this, a new computational modelling approach called “Computational Complementation” has been developed. The main idea is to use functional-structural modelling to complement the deficiency of an empirical model of a loss-of-function (non-AON) mutant with hypothetical AON mechanisms. If computational complementation demonstrates a phenotype similar to the wild-type plant, the signalling hypothesis would be suggested as “reasonable”. Our initial case for application of this approach was to test whether or not wild-type soybean cotyledons provide the shoot-derived inhibitor (SDI) to regulate nodule progression. We predicted by computational complementation that the cotyledon is part of the shoot in terms of AON and that it produces the SDI signal, a result that was confirmed by reciprocal epicotyl-and-hypocotyl grafting in a real-plant experiment. This application demonstrates the feasibility of computational complementation and shows its usefulness for applications where real-plant experimentation is either difficult or impossible.


Introduction
Legumes are one of the largest families of flowering plants that occupy about 15% of Earth's arable surface; yet they provide 27% of the world's primary crop production and more than 35% of the world's processed vegetable oil [1], signifying their cropping potential. Legumes are also the major natural nitrogen-provider to the ecosystem, contributing roughly 200 million tons of nitrogen each year [2] equivalent to over 200 billion dollars worth of fertiliser replacement value. Underlying this powerful fixation capability is a plant developmental process termed ''nodulation'', which results from the symbiosis of legume roots and soil-living bacteria broadly called rhizobia. Yet for a legume plant itself, excessive nodulation may cause over-consumption of metabolic resources and disproportional distribution of internal growth regulators [3], and may interfere with developmentally related lateral root inception and function.
Legume plants have evolved a long-distance systemic signalling regulatory system, known as autoregulation of nodulation (AON), to maintain the balance of nodule formation [3][4][5][6][7]. It has been hypothesised that the induction of the nodule primordium produces a translocatable signal Q, which moves through a root-shoot xylem pathway to the leaves. This Q signal, or an intermediate, is detected in the phloem parenchyma of leaf vascular tissue by a transmembrane leucine-rich repeat (LRR) receptor kinase [8] related in structure to CLAVATA1 in Arabidopsis. This kinase is referred to as GmNARK in soybean [9,10], HAR1 in Lotus [11], and SUNN in Medicago [12]. Q is presumed to be a CLV3/ESR-related (CLE) peptide [13,14]. The perception of the Q signal by the LRR receptor kinase triggers production of a hypothetical shoot-derived inhibitor (SDI) that is transported to the root to inhibit further nodule initiation. SDI can be extracted from wild-type leaves, re-fed via petiole feeding into loss-of-function mutants, resulting in restoration of the wild-type phenotypes [15]. It is a small, water-soluble, heat-stable and inoculation-dependent molecule. However, other mechanisms involved in AON signalling remain largely unknown, though the pre-NARK events (those setting up the signal transmission and then Q signal transduction) as well as the post-NARK events (firstly KAPP phosphorylation, ensuing transcriptional changes, and then SDI production) are being investigated [10,15,16].
To help understand such biological complexities, system modelling has been broadly applied [17][18][19]. From a systematic view, behind the signalling mechanisms is a network of components connected by intricate interfaces, with activities such as ''assembly, translocation, degradation, and channelling of chemical reactions'' occurring simultaneously [20]. These components and their interactions -also responding to the temporally and spatially changing environment -frame dynamic and complex systems at multiple scales to orchestrate plant development and behaviour. As a full understanding of system properties emerging from component interactions cannot be achieved only by ''drawing diagrams of their interconnections'' [17], computational techniques become indispensable for processing massive datasets and simulating complex mechanisms [21].
Although computational approaches have been progressing rapidly for modelling plant signalling, such as for signal transport [22,23], canalization [24] and signalling network [25], most efforts have focused on cellular or tissue levels. Since AON is in essence a long-distance inter-organ regulatory network, our investigation required modelling at the whole-plant scale. Functional-structural plant models [26], such as those developed for resource allocation [27][28][29] and shoot signalling [30][31][32][33][34], can take inter-organ communication into account and use plant architecture as a direct reporter of underlying processes. Functional-structural modelling allowed us to simulate the hypothesised AON signalling and integrate it with nodulation. Yet the major difficulty was not how to model the hypotheses but how to test them through modelling. To meet this challenge, we have developed a new approach -Computational Complementation -for AON study.
Following description of the computational complementation method, we will present its first application in investigating whether wild-type cotyledons participate as an SDI producer in the AON system. Previous studies have indicated that mRNA for GmNARK, which, if translated, is responsible for perceiving the Q signal and triggering the SDI signal, exists in wild-type unifoliate and trifoliate leaves. It is expressed in all vascular tissue [8] of the plant (including the root), but its product is functional only as a nodulation control receptor in the leaf [35]. Thus the RNA expression pattern does not match biological function in AON. Relevant to the investigation here, the vasculature of the cotyledon also expresses RNA for GmNARK; whether this is functional in AON signalling was unclear. Therefore we used computational complementation to test two opposing hypotheses: (a) cotyledons function as part of the root, incapable of perceiving Q and producing SDI; or (b) cotyledons function as part of the shoot, involved in regulating root nodules.

Methods
Genetic complementation [36] is a classical approach to define genetic cause-and-effect relations. For example, assuming two mutant organisms exhibit the same phenotype caused by loss-offunction (recessive) mutations, then their hybrid will be wild-type, if the mutations are in different genes (called cistrons); conversely the hybrid will be mutant if the mutations are in the same cistron. In other words, the wild-type (functional) allele complements the deficiencies of the mutant. Genetic complementation is also used in transgenic analysis of organisms, as a loss-of function mutation in a candidate wild-type gene is deemed causal for a mutant phenotype if that mutant is effectively complemented by the transfer of a dominant wild-type allele. The complementation approach introduced here does not cross one genotype with another, but will use computational modelling to complement the deficiency (in an empirical model) of a mutant to determine if this recovers the virtual wild-type phenotype.
We use two well-characterized soybean (Glycine max L. Merrill) genotypes: the wild-type soybean Bragg and its loss-of-function mutant nts1116 [37]. Wild-type soybean Bragg performs AON to keep its nodulation balance well-maintained ( Fig. 1A and C), leading to characteristic crown nodulation in upper root portions. In its near-isogenic mutant nts1116, the Q signal generated from early nodule proliferation cannot induce SDI due to the lack of GmNARK activity in leaves (Fig. 1B). Reduced SDI in GmNARK-

Author Summary
Endogenous signals, such as phytohormones, play a vital role in plant development and function, controlling processes such as flowering, branching, disease response, and nodulation. However, the signalling mechanisms are so subtle and so complex that details about them remain largely unknown. In this study, we develop a ''Computational Complementation'' approach for the investigation of long-distance signalling networks during legume autoregulation of nodulation (AON). The key idea is to use computational modelling to complement the deficiency of an empirical model of an AON deficient mutant with hypothesised AON components. If the complementation restores a wild-type nodulation phenotype, the modelled hypotheses would be supported as reasonable. To evaluate the feasibility of this approach, we tested whether wild-type soybean cotyledons participate in AON, commonly controlled by ''real'' leaves. The test gave an affirmative result (i.e., cotyledons do have AON activity), which was subsequently confirmed by a graft experiment on real plants. Future applications of this approach may be to test candidate AON signals such as auxins, flavones, and CLE peptides, and other plant signalling networks. This results in a phenotype with a normal number of nodules (C). In the mutant nts1116, GmNARK is not functional in the leaves, leading to the lack of SDI production (B) and consequently a supernodulation phenotype (D) with many more nodules than the wild type. doi:10.1371/journal.pcbi.1000685.g001 deficient plants leads to a phenotype with many more nodules than wild-type, called ''supernodulation'' or ''hypernodulation'' ( Fig. 1D) [5]. Compared with Bragg, the only deficiency of nts1116 plants is the significantly reduced capacity of producing SDI.
The key idea of our complementation approach comes from this point. We ''add'' hypothetical components of AON signalling, including those of signal production, transport, perception and function (see also Text S3), into the empirical model that depicts the growth behaviours of nts1116 plants to see if a wild-type phenotype can be restored. The flowchart of methodology for this approach is given in Fig. 2, including the following steps: Build empirical models to simulate architectural development of Bragg and nts1116 plants based on biometric growth data collected from cultivation of the two genotypes under the same conditions. The empirical data include architectural information such as internode length and  Compare the new phenotype generated by the nts1116+AON model in step (iii) and the nodulation pattern produced by the Bragg architectural model in step (i). If they are same or similar, the hypotheses will be supported as reasonable. Otherwise, the hypotheses need to be modified and tested again from step (iii).
If the hypotheses supported in step (iv) are testable by realworld experiments, the virtual-experiment process can suggest appropriate real-experiment methods to further evaluate them. The mechanisms further supported by realplant experiments will then be used as ''confirmed mechanisms'' in step (iii) to serve the testing of remaining hypotheses.
(vi) If the hypotheses supported in step (iv) are not suitable or possible for evaluation through real-world experiments (due to limitation of current biological techniques), unknown attributes or characteristics about AON signalling can be predicted by virtual experiments.
The architectural and functional-structural models mentioned in steps (i) and (ii) have been built with context-sensitive L-systems [31]. The empirical data used for building architectural models of Bragg and nts1116 plants were collected every second day from growth experiment under the same conditions until the 16 th day post-sowing (all plants were inoculated on the 2 nd day). Materials and methods for this glasshouse experiment are given in supporting Text S1. The growth data, algorithms and techniques used for model construction are described in supporting Text S2. The remaining steps of the flowchart, including (iii), (iv), (v) and (vi), are implemented for hypotheses testing and prediction.

Application of the Approach through Virtual Experiments
In this initial application of our computational complementation approach, two opposing hypotheses were tested: (a) cotyledons function as part of the root, incapable of perceiving Q and producing SDI (abbreviated as ''cotyledon-root'' hypothesis); (b) cotyledons function as part of the shoot, involved in regulating root nodules (abbreviated as ''cotyledon-shoot'' hypothesis). Since GmNARK is expressed in all organs [8] (including cotyledons) and since cotyledons are short-term terminal organs (as they are degraded 7-14 days after germination), neither the cotyledon-root nor the cotyledon-shoot hypothesis was favoured a priori.
Theoretically speaking, if all other AON mechanisms (such as signal production, transport, perception and function) had been confirmed and used as basis for this application, the tested hypothesis leading to a wild-type nodulation pattern could be the correct one. However, the actions of many other signalling components also remain unclear. One or two virtual experiments are obviously insufficient to allow conclusions. Implementing too many experiments (to test all mechanisms together), however, would miss the emphasis and undermine efficiency. With these concerns, our strategy was to adjust parameters for signal production, transport, perception and function within a limited range, and use them as different conditions for different virtual experiments. Among all these experiments, if the complementa-  tion results (nts1116+AON) based on the cotyledon-root hypothesis are always or in most cases closer to Bragg than those based on the cotyledon-shoot hypothesis, then the cotyledon-root hypothesis would be considered plausible; otherwise, the cotyledons are more likely to function as general-sense leaves to regulate root nodulation.
According to this specific strategy, 27 virtual experiments (varying three rates of transport for both Q and SDI and three levels of nodulation inhibitory threshold) were designed for each of the two hypotheses: CRH_1,CRH_27 for cotyledon-root testing and CSH_1,CSH_27 for cotyledon-shoot testing. The only difference between CRH_i and CSH_j, if i = j, is whether cotyledons can function for AON signalling or not. Details of the virtual-experiment assumptions and conditions are described in the supporting Text S3.
To quantify the comparison between complementation results and Bragg phenotype, we define their similarity degree S cp as where N nt , N br and N cp are the nodule numbers generated respectively by the architectural model of nts1116 plants, the architectural model of Bragg, and the functional-structural model of nts1116+AON. This can be understood as the ratio of the number of nodules inhibited by the virtual experiment to the number of nodules inhibited by a real Bragg plant. The similarity degrees of overall nodule number produced by virtual experiments on the 10 th and the 16 th day after sowing are listed in Fig. 3 and Fig. 4, where R q and R sdi represent the transport rates of Q and SDI signals (mm/day). These data indicated that the similarity degrees resulting from cotyledon-shoot hypothesis were generally much higher than those from cotyledon-root hypothesis, supporting the former hypothesis. Considering that values of S cp greater than 100% may mean over-regulation and might not be optimal, the criterion for further evaluating S cp is defined in Fig. 5. According to this criterion, the virtual experiments based on cotyledon-root hypothesis produced unsatisfactory results on the 10 th day (Fig. 3, left-hand column), in sharp contrast to the cotyledon-shoot experiments (Fig. 3, right-hand column). Although there were good results derived from virtual experiments CRH_1, CRH_2, CRH_11 and CRH_13 on the 16 th day (Fig. 4, left-hand column) in terms of nodule number, the nodule size and density from these experiments were all far from similar with the Bragg pattern (Fig. 6). In comparison, the nodule distribution generated by CSH_1 (Fig. 6D) -the opposite of CRH_1 -was quite close to that of the Bragg architectural model.  We predicted from these complementation experiments that the cotyledons should be part of the shoot and participate as an SDI producer in wild-type soybean plants.

Confirmation of the Virtual-Experiment Result
To confirm the above prediction and also to evaluate the effectiveness of this approach, a ''real-plant'' grafting experiment was conducted. The critical experiment was to graft -between Bragg and nts1116 plants -the shoot of one genotype with cotyledons to the root of the other genotype without cotyledons, and also to graft the shoot of one genotype without cotyledons to the root of the other genotype with cotyledons, forming four graft combinations: Ns+Nc+Br, Ns+Bc+Br, Bs+Bc+Nr and Bs+Nc+Nr ( Table 1). Materials and methods for this graft experiment are given in the supporting Text S1. The collected empirical data for nodule number were not only classified by each graft type but also according to each plant's cotyledon retention status ( Table 2).
According to the experimental results, the nodule number from the Ns+Nc+Br graft type was much higher than that from the Ns+Bc+Br (Fig. 7A). For the Ns+Bc+Br graft type alone, its plants with fallen cotyledons had more nodules than those with persisting cotyledons, and the plants with yellow cotyledons had more nodules than those with green cotyledons (Fig. 7C). These differences suggest Bragg cotyledons were the only leaves to regulate nodulation in Ns+Bc+Br plants, because unifoliate and trifoliate leaves of nts1116 plants were unable to do so.
Data of another graft type with Bragg cotyledons -the Bs+Bc+Nr (Fig. 7D) also suggested that the Bragg cotyledons participated in providing SDI. However, more nodules were found in the Bs+Bc+Nr plants than in the Bs+Nc+Nr plants that had no Bragg cotyledons (Fig. 7B). An explanation for this observation is that the Bs+Nc+Nr allowed more nodules to be formed at early stages than the Bs+Bc+Nr, leading to more Q signal moving from root to shoot. As the cotyledon biomass declined greatly at later stages of seedling growth (resources are unloaded for plant growth and the ''spent'' cotyledon is eventually discarded), the difference in shoot between Bs+Bc+Nr and Bs+Nc+Nr became insignificant. Therefore larger amounts of Q triggered more SDI, which finally inhibited more nodules in Bs+Nc+Nr.
To better understand this nonlinear characteristic brought out by real-plant experiments, we returned to the virtual-experiment models and visualised the dynamic signal allocation during CRH_1 and CSH_1 (Fig. 8). As demonstrated by the visualisation, the SDI concentration (in the root) of CRH_1 was lower than that of CSH_1 on the 5 th day but became higher from the 10 th day on, in agreement with the above analysis of the nodulation difference between Bs+Bc+Nr and Bs+Nc+Nr. Thus, we conclude that the testing result from our initial application of computational complementation is confirmed: the cotyledons ''belong'' to the shoot and function as a source of the nodulation regulator in wildtype soybeans.

Discussion
The computational complementation approach introduced here is an original contribution to the study of legume autoregulation of nodulation. Compared with conventional biological technologies with broader implications to plant development, one of the major advantages of this approach is its capability to complement the deficiency of a mutant plant at an organ scale with totally hypothetical and concept-derived physiological components. It is also able to make hypothetical signalling details manipulable and visible. For example, as demonstrated in the above case, signal transport rates can be modified as hypothesised and the allocation of signal can be dynamically visualised. These functionalities not only enable AON researchers to test hypotheses or make predictions using time-and resource-saving virtual experiments, but also bring out possible underlying details that are unobservable through real-plant experiments. Moreover, the application of this approach is not only limited to AON research, but also potential to other plant signalling studies such as those on branching regulation (e.g., [38]), flowering control (e.g., [39]) and lateral root initiation (e.g., [40]). This approach contributes a new idea to the domain of computational plant modelling -computational complementation. From a classic modelling point of view, one can formulate a model based on empirical data and then verify the model against the data, which has been used for development of crop (e.g., [41]) and architectural (e.g., [42]) models. However, what we investigate is a largely unclear internal signalling system -most of the detailed mechanisms remain unknown, which determines there is no direct parameterisation-and-verification data to evaluate the modelled signalling hypotheses. Using an indirect strategy, functionalstructural modelling allows us to use the observable structure as a reporter for estimation of the unobservable function. But for this study, we have to link the structure of one genotype with the function of another genotype. The reason for this is: the wild-type Bragg nodulation has already been regulated, thus incorporating AON to Bragg architecture would double the regulation and have no reasonable comparison target for validation; in contrast, the nts1116 is a non-AON plant and this is its only difference with Bragg, therefore activating AON in nts1116 plant could result in system behaviours comparable with the wild type.
Another feature of this approach resides in the level of complexity for simulation of structural and signalling processes. We captured root details for studying shoot-root signalling rather than oversimplifying the root system. And the signalling pathways are constructed with sub-modules of which the size and number can be manipulated without limitation, which allows future modelling work to be extended to lower-scale mechanisms (such as tissue and cellular scale). We also created a synchronisation algorithm for coordination of multi-rate procedures to enhance the precision of signalling-development interactions. A description of these modelling techniques is given in the supporting Text S2.
The approach also has some limitations. For example, due to the nature of complementation, it can only be used for a single mutation at a time, though leaky mutants can be handled by parameter optimization. Another drawback is that it cannot distinguish between different mutations in the same pathway that result in the same phenotype in the first instance. In other words, if the hypothesised mechanisms used to complement the mutant are the same in both cases, and so is the phenotype of the two mutants, computational complementation cannot be used to say which gene component of the regulatory network has been mutated.
Our first application of this approach was to test whether wildtype soybean cotyledons are involved in production of SDI. Also but more importantly, we expected this application to evaluate whether the computational complementation idea is effective. The virtual-experiment results suggested the wild-type cotyledons can produce SDI, which was further confirmed by a graft experiment on real plants. This demonstrates the feasibility of computational complementation and shows its usefulness for future applications. The next step is to apply this approach to support research for the identification of Q and SDI. Candidate signals, such as CLE peptide for Q [13,14] and auxin for SDI [43], will be tested to see if they play the roles in AON as hypothesised. In addition, environmental factors, such as soil nitrogen status, that have effects on the process could also be tested with this approach. Furthermore, the finding that wild-type soybean cotyledons act as an SDI producer in AON opens the door for testing physiological transgenerational effects, such as altered nodulation patterns influenced by the Bradyrhizobium infection status of mother plant through presence of SDI in cotyledons.

Supporting Information
Text S1 Materials and methods for glasshouse experiments