text:
How Many Parameters Does It Take to Describe Disease Tolerance?

The study of infectious disease has been aided by model organisms, which have helped to elucidate molecular mechanisms and contributed to the development of new treatments; however, the lack of a conceptual framework for unifying findings across models, combined with host variability, has impeded progress and translation. Here, we fill this gap with a simple graphical and mathematical framework to study disease tolerance, the dose response curve relating health to microbe load; this approach helped uncover parameters that were previously overlooked. Using a model experimental system in which we challenged Drosophila melanogaster with the pathogen Listeria monocytogenes, we tested this framework, finding that microbe growth, the immune response, and disease tolerance were all well represented by sigmoid models. As we altered the system by varying host or pathogen genetics, disease tolerance varied, as we would expect if it was indeed governed by parameters controlling the sensitivity of the system (the number of bacteria required to trigger a response) and maximal effect size according to a logistic equation. Though either the pathogen or host immune response or both together could theoretically be the proximal cause of pathology that killed the flies, we found that the pathogen, but not the immune response, drove damage in this model. With this new understanding of the circuitry controlling disease tolerance, we can now propose better ways of choosing, combining, and developing treatments.

Introduction
The clinical goal of treating infectious diseases is to reduce the levels of sickness experienced by infected hosts. One approach to studying this problem is to quantitate illness by correlating the dose response of pathology to pathogen load; graphs like these are called "disease tolerance curves" [1][2][3][4][5][6]. We argue that the shape of the dose response curves and the underlying mathematical functions producing these curves can teach us how to alter the health-by-microbe relationship. For example, if disease tolerance curves are linear, we need only discover the molecular mechanisms that control the curves' slopes. If these curves have a more complex shape-for example, a sigmoid shape-then we will need to measure more variables.
Since there are not universally agreed upon definitions for many of the words we use to describe the immune response and health, we start by defining the following terms: vigor, resistance, resilience, and tolerance. "Vigor" is the health of an uninfected individual, and "resistance" describes the host's ability to limit microbe load [1][2][3][4][5][6]. We use "resilience" to describe the ability of infected hosts to return to their original healthy condition, much in the same way the word is used to describe the way a perturbed ecological system returns to its origin [7,8]. For example, a resilient host may suffer from an infection but would easily bounce back to health, whereas infection in a non-resilient host might lead to permanent disability or death. We apply the term resilience to an infected individual because we can follow the path that individual takes from health through sickness and back. "Disease tolerance" is the dose response curve summarizing the pathogen loads that are required to produce a range of host responses. Disease tolerance differs from resilience in that tolerance is an emergent property of populations, while resilience is a property of an individual; disease tolerance measures how resilience changes across a population as pathogen load is varied. Disease tolerance, by definition, cannot be measured in a single individual [9]. The model described in this paper shows how these distinctions between resilience of individuals and tolerance of populations disappear once we understand infection dynamics to the point that we can predict how a population would behave knowing the properties of an individual.
To study disease tolerance we used a simple infection system in which we injected a pathogen (Listeria monocytogenes) into a fruit fly (Drosophila melanogaster) [10][11][12][13][14]. These L. monocytogenes-infected flies suffer from a variety of ailments, including alterations in their circadian rhythm, cold-coma recovery, climbing ability, feeding behavior, metabolism, and death [12,[15][16][17]. To determine the disease tolerance curve, we plot the survival rates of the infected flies against microbe loads measured 2 days post infection (DPI). We used this system because it allowed us to measure large numbers of samples and to take advantage of the genetic tools and the understanding of innate immunity available for this model organism.
Our immediate goal is to define the mechanisms controlling tolerance. We start with a simple mathematical model representing feedbacks between microbe growth, the immune response, damage, and health, modeling each using logistic equations. We then measured the ground-truth of these models using a D. melanogaster/L. monocytogenes infection that can lead to lethal outcomes. We experimentally measured microbe growth, the immune response, and disease tolerance and found that each of them was well described by sigmoid curves, suggesting that we need at least three variables to describe each. In the case of tolerance, this meant that we had to measure host vigor, the number of bacteria required to injure the host, and the maximum death rate achieved by the fly. In the case of antimicrobial peptide (AMP) transcript production, we measured basal levels, the number of bacteria required to induce transcripts to 50% of their maximum value, and the maximum transcript level. We examined how alteration of the values of the parameters in the model changed the output of the model to distinguish between changes that we might see if microbes or the immune response were the principal mediators of damage. Upon testing a variety of D. melanogaster natural variants and mutants along with some L. monocytogenes mutants, we concluded that, in this system, the bacteria and not the immune response is responsible for the damage caused by the infection. When we monitored how disease tolerance changed as we varied host or pathogen properties, we found the tolerance curve changed shape in the manner predicted by the model; for example, the number of bacteria required to damage health changed, while the vigor and maximum death rate remained constant. We demonstrate that this dual approach of measuring and modeling full-length disease tolerance curves can reveal previously unrecognized parameters controlling disease tolerance.

Results
We started by making a graphical model to explain four parameters: microbe load, the immune response, microbe-induced damage, and health (Fig 1). We modeled the immune response so that it could both limit microbe growth and kill microbes. We previously found that Mycobacterium marinum infected flies waste upon infection and predicted that the day they died they would have exhausted their glycogen and fat stores [18,19]. We also found that L. monocytogenes infected flies waste during infection [15]; to explain death, we used these empirical data to imagine that there was a store of "health" that could be depleted by the infection. In our model, both the immune response and the microbes can cause damage, which depletes health and increases the death rate. Immunopathology directly affects health, while microbes secrete damage effectors that impact health.
We built a simple mathematical model based on the graphical model largely using sigmoid functions to define each parameter (see Experimental Procedures) to help define our assumptions about the mechanisms governing disease tolerance. This was an obvious choice for microbe growth and the immune response, as these equations are well described in classical mathematical models of infections [20]. We used the same mathematical function to describe microbe-induced damage, since we anticipate that the shape of the curve is due to some enzymatic rate-limiting step that will likely follow sigmoid Michaelis-Menten kinetics [21,22].
To determine the ground-truth of this model, we measured the behavior of microbes, the immune response, and health in a model system in which we inject L. monocytogenes into the body cavity of D. melanogaster. We chose this model infection because it is simple to obtain the measurements we need to define resistance and disease tolerance and because we have a collection of host mutants that we know affect resistance and tolerance to L. monocytogenes [10][11][12]. We can readily measure microbe loads by homogenizing whole flies at any point in the infection to plate them to determine colony-forming units (CFU). We can determine the health of the infected flies by measuring the median time to death following infection.
When L. monocytogenes is injected into flies that are incubated at 29°C, the resulting bacterial growth is best described as logistic (Fig 2A, logistic adjusted r 2 = 0.594, adjusted linear and adjusted exponential r 2 < 0.1) with a maximum growth rate of 0.1656 log10 (CFU)/h. L. monocytogenes loads reach a plateau at 24 h post infection.
We wanted to define the relationship between AMP levels to a given microbe load to understand how many parameters we need to use to measure the immune response. Because bacteria grow rapidly to reach their plateau, and wounding flies to introduce bacteria causes an immune response on its own, we reasoned it would be difficult to measure the immune responses at intermediate microbe loads by simply injecting flies with different amounts of microbes. Instead, we took flies that had high levels of bacteria and slowly cured them of the infection so that we could measure the relationship between AMP transcripts and microbe load essentially at steady state. The fly immune response cannot clear L. monocytogenes on its own; therefore, to clear the bacteria we transferred infected flies to fly media containing 1 mg/ml of ampicillin   DPI. We observed that L. monocytogenes clearance followed exponential decay dynamics ( Fig  2B, adjusted r 2 = 0.7140). L. monocytogenes infection in the fly decreases survival, but ampicillin-treated flies lived just long as uninfected flies (Fig 2C).
Having developed a method to slowly decrease microbe loads, we performed the following experiment: We injected one set of flies with phosphate buffered saline (PBS), fed them with ampicillin, and followed them regularly to serve as wounding, aging, and antibiotic treatment controls. We took a second set of flies and injected L. monocytogenes. This group was split in two; one portion was fed ampicillin to clear the microbes, and the other was not. We collected flies for plating to determine microbe loads, and for RNA preparation to run a microarray at regular intervals (S1 Fig).
We analyzed the expression of the known and presumed AMPs with respect to microbe loads and found that AMPs are expressed according to a sigmoid relationship with respect to microbe load. We plotted the expression of AMPs to CFU levels in the fly and found the relationship between AMP transcription and microbe load is best fit by a four-parameter sigmoid ( Fig 3A). Accordingly, each AMP transcript has a baseline (bottom plateau), maximal concentration (top plateau), effective bacterial concentration producing a half maximal effective concentration (EC 50 ), and slope. We observe that the AMPs expressed during infection differ in their EC 50 ( Fig 3B) over a range of about 100-fold (20-2,000 CFUs). AMP expression does not continue to increase as microbe number increases past a certain point, and each AMP transcript reaches a characteristic plateau. Each AMP transcript can thus be characterized by its maximum expression in relation to its baseline, and this relationship varies over an 8-fold range between AMPs. We observed no significant relationship between EC 50 and maximum expression, suggesting that EC 50 and maximum expression can vary independently of each other. We found an exponential relationship between Hill slope and EC 50 (Fig 3D), suggesting that AMPs requiring large numbers of bacteria for induction turn on in a switch-like manner, while those turned on at low microbe loads are turned on gradually. This gives us insight into the tuning of the fly's immune response; overall, AMPs will be turned on gradually at microbe loads below 100 bacteria per fly but will rapidly increase as loads pass 1,000 bacteria per fly.
Though AMPs had a rather simple sigmoid relationship with microbe load, this was not the case for many other genes modulated during infection. Some of these other gene expression patterns showed hysteresis, by which we mean that gene expression depends upon past conditions and not only the immediate conditions. For example, gene expression for these hysteretic transcripts had one pattern as the fly moved from health to sickness and another as the fly returned from sickness to health. This hysteresis can be useful because it helps us define different states of the infection. To visualize the "disease space" traversed by infected flies and to identify the different states of the infection process, we used topological data analysis (TDA). It is common practice when clustering data to plot the data using some sort of tree structure, as is done with hierarchical clustering or spade analyses [23,24]. This makes good sense when dealing with data for a developing system in which the animal or cell enters a new state, but it is problematic when dealing with a system that returns to its original state; resilient systems should not be fit by a tree and are better described by loops. Topological data analysis is sensitive to the topology of the data and will not arbitrarily linearize a loop and then force it to fit a tree [25,26]. Instead, the analysis simply clusters related data points, represented as nodes on a network graph, and the shape of that graph reveals the connections between the time points flies demonstrate that antibiotic treatment led to recovery of the flies, which lived as long as wounding control flies. Data are reported in S1 Data.  [27][28][29]). In the case of a resilient system, such as flies recovering from an infection, we find this graph forms a loop. Fig 4 shows the TDA graph we built to describe our data, and 4A describes the three basic treatment groups in the dataset (uninfected, infected, and antibiotic-treated). In Fig 4B, the phase curve is colored by CFU (blue indicates low and red indicate high microbe loads). The green arrow indicates progression from health to infection and back, moving clockwise around the figure. Uninfected flies are at the top of the figure, and the figure progresses through acutely wounded flies, sick flies, and on to recovering flies, which link back to the original uninfected flies. Infection was initiated by the injection of 100 CFU, and L. monocytogenes levels increase following this. The path followed by moribund flies deviates from recovering flies to form a spur on the loop (marked as section iv). Flies that received ampicillin decrease in CFUs, and flies that have low CFUs loop back to overlap with uninfected controls.
To determine which genes show altered expression in different parts of the disease map, we performed fuzzy c-means clustering on this dataset to identify genes with similar expression patterns (S2 Fig, S2 Table). We then performed gene ontology analysis on each cluster to identify characteristic physiological changes within the cluster (S3 Table) and then asked about the expression patterns of these genes on the TDA graph. Acute stress response genes, such as heat shock protein 26, were activated immediately after injection and then quickly returned to  baseline levels (Fig 4B and 4C). AMP expression follows CFU (Fig 4B and 4D). We identified groups of genes that are highly expressed by flies in the spur approaching death (Fig 4B and  4E) and recovering from infection (Fig 4B and 4F). The transcripts for these "morbidity genes" (clusters 1,2,11,12,22,35,43) continue to increase even when microbe levels plateau. Gene ontology analysis of fuzzy c-means clusters containing these morbidity genes did not yield significantly enriched terms. Among the top 50 differentially expressed morbidity genes, Bteb2 and run are predicted to encode proteins with DNA binding activity. CG3117 has predicted protease activity, and diedel is a putative negative regulator of JAK/STAT signaling. By contrast, the recovery genes (clusters 18,20,23,25,29,30) are enriched for sugar and fatty acid metabolism genes that are turned up above baseline upon recovery.

The Shape of a Disease Tolerance Curve
To measure a disease tolerance curve, we recorded the response of the host to a broad range of initial pathogen doses. We did this by injecting L. monocytogenes into flies over a range of ten to 100,000 bacteria and allowing the infected flies to die. We injected L. monocytogenes into the hemocoel of flies and monitored bacterial numbers 2 DPI to measure the ability of the fly to resist microbe growth when challenged with a range of infection intensities (Fig 5). We determined the median time to death (MTD) for each inoculum and used this time as a measurement of health ( Fig 5B). Plotting microbe load versus MTD produced a curve that was readily fit by a four-parameter logistic sigmoid model (r 2 > 0.96) (Fig 5C).
Occasionally, tolerance curves are fit with mathematical functions, but the reason these functions are chosen is that they fit the data and not that they are dissected for further biological insights [5,[30][31][32]. More typically, tolerance is visualized as a linear system, which requires just two parameters, vigor and slope [5,11]. In contrast, our sigmoid model provides four parameters ( Fig 5D). In addition to vigor and slope used to describe a line, the sigmoid model adds the parameters of EC 50 and maximum effect. The EC 50 of the system is the number of microbes present at day 2 that cause a 50% change in MTD. Maximum effect is defined by the sigmoid's asymptotic tail at high microbe loads. This defines a maximum death rate, suggesting that there is a previously unrecognized rate-limiting step for death.
In the model shown in Fig 1A, depletion of health could be induced either by bacterial damage effectors or indirectly by self-harm caused by the resistance response. We modeled the two possibilities by observing how the shape of the curve changed as we altered the model such that either the resistance response or bacteria were the sole cause of damage. We concentrated on changes in the rate that the immune response was turned on (τ), and the inflection points for the relationships between microbe density and the rates of immune induction and microbeinduced damage induction (σ I and σ M ) (Fig 6A-6D, S4 Table). In the case in which resistance mechanisms drove damage (Fig 6A and 6C) and bacterial damage was set to zero, changes in τ or the inflection point for microbe-induced immunity caused shifts in both the EC 50 and microbe loads. The effect is more extreme when the inoculum is far below the microbial carrying capacity, as this gives the immune response an opportunity to control microbe loads. This results in a reciprocal mechanistic link between resistance and tolerance in which one is always high when the other is low. In the model in which bacteria drive damage, a loss of resistance results in high microbe loads and health pegged at maximum severity ( Fig 6B). Shifts in the EC 50 of the bacterial damage-driven system, in which resistance-induced damage is set to zero, are caused by changes in the inflection point for bacterial-induced damage, as might be expected for differentially virulent strains of microbes or hosts that are better at neutralizing bacterial toxins (Fig 6D). In this second case, there is no mechanistic reciprocal link between resistance and tolerance, and the two vary independently of each other.

Testing the Mathematical Model for Disease Tolerance
To test whether the shape of the disease tolerance curve in this host-pathogen pairing was driven by immunological damage or microbial damage, and to determine how the parameters of a sigmoid disease-tolerance curve varied, we measured the tolerance curves for a collection of Drosophila mutants that we previously showed differ in their resistance and tolerance defenses to L. monocytogenes [10][11][12]. We also tested natural variant flies from the Drosophila Genetics Reference Panel that were preselected because they showed extreme changes in their response to L. monocytogenes infections [33,34]. In addition, we tested L. monocytogenes mutants that altered pathogenesis in the fly [11]. These data supported the model in which bacteria caused damage in the system and did not support the immunological damage model.
A mutation in the D. melanogaster gene CG2247 was found previously to reduce the fly's ability to both survive an infection and control L. monocytogenes growth [11]; here, we found that at 2 DPI, even CG2247 mutant flies injected with just ten bacteria had the same high microbe loads as flies injected with 100,000 bacteria, supporting the idea that these mutants alter resistance. This is visible in the tolerance curves, in which all of the points from infected mutant flies clustered on the bottom asymptote of the parental tolerance curve (Fig 6E-6H, S5 Table). An additional example of a fly strain with poor resistance is shown in the supplementary data (S3 Fig and S5 Table). Since flies suffering a loss of resistance showed no change in the maximum death rate of the infection, these results support the model in which bacteria are responsible for damage.
RNAseq analysis of CG2247 mutant flies suggests a molecular mechanism for the observed phenotype. We observed that in wild-type flies, the transcripts encoding the enzymes required to generate a reactive oxygen response against L. monocytogenes dropped over the course of the infection but recovered upon antibiotic treatment (S6 Fig) [11,35]. In contrast, these transcripts dropped to 10-fold lower levels in CG2247 mutants and did not recover to their original levels. We observed a reduction in the characteristic dark spots resulting from the action of this immune response (S6 Fig). We conclude that mutations in this gene disrupt the melanization immune response.
The disease-tolerance curve for the natural variant D. melanogaster strain, RAL 359, showed a reduction in the EC 50 of the system without a change in resistance [33]. RAL 359 was as capable as our lab control strain w 1118 in controlling L. monocytogenes growth; however, low doses of the microbe had a larger effect on health in this strain than in the control. This shifted the EC 50 from 922 to 83 colony-forming units per fly (Fig 6I-6L and S5 Table). The vigor and maximum death rate of these two strains was similar. This shift in EC 50 was a common phenotype, as shown in S4 and S5 Figs.
Changing the virulence of the L. monocytogenes also altered the EC 50 of the system. Mutant L. monocytogenes lacking the virulence factors actin assembly-inducing protein (ActA) or listeriolysin O (LLO) had previously been shown to kill flies more slowly than wild-type L. monocytogenes (Fig 5C versus Fig 6O) [11]. The tolerance curves for flies infected with an ΔactA mutant showed a shift in the EC 50 from 922 to 5,753 (Fig 6M-6P and S5 Table). We observed a ceiling for the number of L. monocytogenes that can be maintained in a wild-type infected fly (Figs 2A and 5A) and hypothesize that this is defined by the number of phagocytes in the fly that provide a niche for L. monocytogenes growth [13]. The ΔactA mutant bacteria approached this limit before they induced maximal pathology and, thus, we did not observe a low asymptotic tail. Δhly mutants produced a similar result to ΔactA mutants, only more extreme (S4

Discussion
As a field, we can generate complex network diagrams showing how the parts of the immune system communicate with each other and identify key signaling pathways and nodes [36,37]. Though this lets us predict that the system will fail when we remove a necessary component of the network, evolution and medicine work by subtler means, for example, by changing rate constants in addition to deleting entire pathways. To predict how evolution or medicine will change an immune response, we need to identify mathematical functions that accurately describe the behavior of that process so that we can understand all of the parameters that control the behavior.
L. monocytogenes growth in the fly follows logistic kinetics; the variation of these curves observed in natural fly variants suggests that microbes are controlled by separable parameters defining the growth rate and carrying capacity of the host [34]. It seems reasonable that this is not a property limited to L. monocytogenes growing in the fly and that it will be true for many pathogens, especially considering that this is a near textbook description of microbial growth [20,38,39]. Both of these properties should be under natural selection, but studies of infections in D. melanogaster typically do not measure microbe growth at all and, when they do, they measure only one time point, conflating growth rates with carrying capacities. The fly is studied because it is easily manipulated and has a relatively simple immune response, but despite decades of work, the field has been underestimating the number of dimensions in which natural variation can change the immune response. Thinking about this generally, it raises the question of whether there might be separate host mechanisms that control microbe growth rates versus ceilings and, thus, diseases might result from a defect in either of these mechanisms. To study these different disease mechanisms, we must ensure that our experiments let us observe the phenomenon.
The AMP response in the fly is under the control of two signaling pathways, Toll and imd [36,40]. The classical descriptions of these pathways argue that the pathways show specificity for different elicitors; a problem with the experiments supporting this model is that the experiments did not correlate response to elicitor levels, and most experiments tested just one arbitrarily determined time point [41,42]. We find that AMP induction is described well by a fourparameter sigmoid function. This means that there are parameters describing the basal level, EC 50 , maximum expression level, and slope, which may all vary independently of each other; that we don't see a correlation between maximal gene induction and EC 50 supports this hypothesis. If the field's experiments measure gene expression at just one time point post infection and do not correlate this to microbe loads, we cannot learn which of the three induction parameters are changing. Thus, though it is clear that the Toll and imd signaling pathways, as well as sex, age, the environment, and evolution, can change the immune response, we don't have a good definition of specificity in this system. Now that we understand that there are new variables to explain, we need to explain this from an evolutionary perspective; for example, what selective pressures cause changes in the EC 50 of gene expression of an AMP rather than the maximal gene expression or slope? Will a host benefit from a quick induction of antimicrobials at low microbe loads, or should it invest in an enormous response as microbe loads climb higher?
We found that health was correlated with microbe load according to a four-parameter sigmoid function. As described above for microbe growth and immune induction, past experiments looking at health outputs in model infections are problematic because they don't measure the full shape of the response curve and don't report where an experiment is performed on the curve. Perhaps it should have been obvious that microbes would have a growth ceiling and AMPs would have a maximum expression level; classic infectious disease modeling texts start with chapters describing these functions [20,43]. However, one thing these books don't do is link pathogen growth and the immune response to health. When we examine this relationship, we find that we can make simple models that predict the presence of more parameters than we've been studying. For example, the implication of a sigmoid response for health suggests that there is a rate-limiting step for death and that death has a maximum rate. This means that it should be possible not only to change the EC 50 of the system to reduce the rate at which damage occurs, but we should also be able to increase the health reservoir of a host so that it has high stamina and will be able to suffer damage for longer periods of time.
The underlying problem with the field's experimental characterization of immune responses is that we focus on when responses are induced in terms of time and not why they are induced in terms of the quantity of inducer. The result is that we systematically overlook and conflate control mechanisms. By correlating effects to effectors, we can better understand microbial pathogenesis.

Disease Spaces
The disease space analysis of recovering and dying flies reveals that far more is going on in these sick animals than AMP gene expression, and it lets us organize these events with a simple map. A commonly used time point for gene expression studies in Drosophila is 6 h post infection [44,45]; this is a particularly problematic time point because it combines transient expression of heat shock genes with expression of antimicrobials. Flies fated to die have a progressing gene signature in which a set of transcripts is increased without a corresponding increase in microbe load. We anticipate that one will need to measure the course of pathogenic infections to identify these genes, rather than follow microbes that simply elicit an immune response but do not kill wild-type flies (for example Escherichia coli or Micrococcus luteus [40,46]. Just as each pathogen causes a different disease in humans, we expect that there will be a range of pathologies caused by insect pathogens and do not assume that these L. monocytogenesinduced genes are the only morbidity genes in the fly. The nature of this particular death response is unclear, as the identity of the genes does not provide obvious clues. In contrast, recovering flies have a unique gene expression signature that suggests function; for example, enzymes in biosynthesis and energy metabolism are repressed during infection but pop up hysteretically as microbes are removed, suggesting the induction of a recovery response that rebuilds damaged tissues.

Why Does Shape Matter?
The shapes of tolerance curves are important because they will define which medications can be used to improve the health of an infected patient depending upon their position on the sigmoid curve. Here is a concrete example: If we think about the dose response curve to sepsis in humans as a linear response, then any drug that limits damage will help every patient. We come to a different conclusion if we consider a sigmoid relationship. Humans suffer from sepsis when relatively small numbers of microbes enter the blood, and very sick patients can have much higher levels of microbes. This suggests the sickest patients will be found on the asymptote of a sigmoid curve, where they will be suffering maximum pathological effects. Drugs that change the basal level of a response, EC 50 , or slope will not help these patients; they will only respond to drugs that manipulate maximal response levels. Thus, we need to understand the shape of these tolerance curves and where the patients lie on the curves to select appropriate treatments.
The discussion of tolerance in patients raises a tricky point. We've shown it is possible to measure a disease tolerance curve in a model system; this demonstrated the nonlinear nature of tolerance and suggested the existence of new dimensions we should follow when measuring disease outcomes. The problem is that this approach requires us to perform many more experiments than we do currently when analyzing phenotypes, and this approach is likely impossible to apply to most human infections for the following two reasons: First, disease tolerance is plotted using summary characteristics of an infection. For example, one might measure the maximum parasite load and minimum health for an infection [5]. This is accessible experimentally, but it is unethical to gather these data from a patient, as you need to treat the patient when they enter the clinic. In one rare case of HIV, infected patients' tolerance curves have been recorded because these chronically infected patients are followed for years, but that approach isn't going to be useful for acute infections [47]. Second, disease tolerance is a measurement of the behavior of a population and not an individual. The information we gather about an individual allows us to place a datum on a health by microbe graph, but we don't know the shape of the curve that should be fit through that individual [9]; thus, we can't easily determine which parameters need fixing in the sick patient. What we learn from model systems is that there are underlying rules that can be used to explain disease outcomes, and that these rules are nonlinear. By studying these rules in models, we can identify the molecular mechanisms corresponding to the mathematical parameters and then apply this knowledge to human biology by analogy. If we want to measure changes directly in human immune infections, we need to rely on descriptions of disease space, which are described in the accompanying paper [48].

Injection
Bacteria were injected into flies essentially as described previously [10,12,14,49]. Flies were anaesthetized with carbon dioxide. A drawn glass needle carrying L. monocytogenes was used to pierce the cuticle on the ventrolateral side of the abdomen. A picospritzer III was used to inject 50 nl of liquid into the fly. Bacteria were delivered at different concentrations to produce injections of approximately 10, 100, 1,000, 10,000 or 100,000 CFU. Infectious doses were determined for each experiment by plating a subset of flies at time zero. Approximately 200-400 flies were used for each dose in the experiment to measure survival and to count colonies.

Bacterial Strains and Culturing Conditions
All L. monocytogenes stocks: Wild type/mutant parental strain (10403S), Δhly (DP-L2161) [50], and ΔactA (DP-L3078) [51] were stored at -80°C in brain and heart infusion (BHI) broth containing 25% glycerol. To prepare L. monocytogenes for injection, bacteria were streaked onto Luria Bertani (LB) agar plates containing 100 ug/mL streptomycin and incubated at 37°C overnight. Single colonies of L. monocytogenes from the LB agar plate were used to inoculate 4 mL of brain and heart infusion (BHI) broth and incubated overnight at 37°C without shaking. Bacteria were removed from the incubator at log growth phase. Prior to injection, L. monocytogenes cultures were diluted to the desired optical density (OD) 600 in phosphate buffered saline (PBS) and stored on ice.

Bacterial Plating
Single flies were homogenized in PBS using a motorized plastic pestle in 1.5 ml tubes. The supernatants were plated using an Autoplate spiral plater and counted using a Qcount automated counter. At least six samples were counted to determine the median number of bacteria for each inoculum. Bacteria were plated onto LB medium and incubated overnight at 37°C before counting.

Survival Curve Analysis
200-400 flies were injected and checked daily to measure mortality for each inoculum. Flies were housed in vials containing 20-25 flies each.

Repetitions
Each set of conditions was repeated at least three times; for example, an experiment for a mutant fly line would be set up independently on three different days to gather microbe load and survival data. Pairs of plating and survival data from these multiple experiments were all plotted on the same tolerance curve.

RNA Isolation for Microarray
Flies were injected with 50 nl of Listeria (OD 600 = 0.001, or approximately 100 CFUs) or left manipulated. Following injection, flies were placed in vials containing dextrose fly media and incubated at 29°C. Samples were separated into three groups: Moribund, recovering, and uninfected control. Moribund flies were infected, and samples were collected on the following DPI: 0.25, 1, 2, 2.25, 3, 4, 5, and 6. Recovering flies were infected and flipped onto dextrose fly media containing 1 mg/ml ampicillin 2 DPI. Recovering samples were collected on the following DPI: 2. 25, 3, 4, 5, 6, 7, 9, and 16. Uninfected control flies were not infected but were flipped onto dextrose fly media containing 1 mg/ml ampicillin 2 dpi. Uninfected control samples were collected on the following dpi: 2.25, 3, and 9. At each of the indicated time points, groups of 20 flies were homogenized in TRIzol, and RNA was isolated using a standard TRIzol preparation. Additional flies were used to determine CFUs and monitor survival. Biological triplicates were obtained from three independent experiments. Quality of RNA was determined using a BioAnalyzer 2100. Samples were labeled using the Quick Amp Labeling Kit, One-Color (Agilent), and hybridized to 4x44K (V2) Drosophila Gene Expression Microarray (Agilent) using the Agilent Gene Expression Hybridization kit following the manufacturer's protocol. Microarrays were scanned using an Agilent Technologies Scanner, and processed signal intensities were determined using Agilent's Feature Extraction software. RNA quality assessment, labeling, hybridization, and microarray feature extraction were performed at the Stanford Functional Genomics Facility.

Microarray Analysis
The microarray data were analyzed using Genespring v12.1 (Agilent). Microarrays were normalized to the 75th percentile. The median expression level on 0 dpi was set as the baseline. Prior to statistical analysis, low-quality spots were removed based on flag calls. To determine differential gene expression, one-way ANOVA was performed comparing all samples to 0 dpi. P-values were corrected by the Benjamini-Hochberg method. Genes with p-values <0.05 at any time point with a fold-change greater than 2 were categorized as differentially expressed. Differentially expressed genes were then clustered using Mfuzz [52].

Topological Data Analysis
Only genes that were differentially expressed in moribund and recovering but not in uninfected control were used in the topological data analysis using the Ayasdi 3.0 software platform (ayasdi.com, Ayasdi Inc., Menlo Park, California). Nodes in the network represent clusters of samples of infected flies, and edges connect nodes that contain samples in common. Nodes are colored by the average value of their samples for the variables listed in the figure legends. TDA was used to map the way hosts loop through the disease space in an unsupervised fashion. Two types of parameters are needed to generate a topological analysis. First is a measurement/ notion of similarity, called a metric, which measures the distance between two points in some space (usually between rows in the data). Second are lenses, which are real valued functions on the data points. Lenses are used to create overlapping bins in the dataset. Overlapping families of intervals are used to create overlapping bins in the data. Metrics are used with lenses to construct the Ayasdi 3.0 output. Multiple lenses can be used in each analysis. There are two parameters used in defining the bins. One is resolution, which determines the number of bins; higher resolution means more bins. The second is gain, which determines the degree of overlap of the intervals. Once the bins are constructed, we perform a clustering step on each bin, using single linkage clustering with a fixed heuristic for the choice of the scale parameter. This heuristic is described in [26]. This gives a family of clusters within the data, which may overlap. We built a network with one node for each such cluster and where we connect two nodes, if the corresponding clusters contain a data point in common. We used two types of lenses. The first type was lenses based on dimension reduction algorithms such as multidimensional scaling and nearest neighbor analyses; these helped analyze the data in an unsupervised manner. The second type was lenses based on the data alone, for example, the levels of antimicrobial gene expression or recovery gene expression. The gene expression markers helped us dissect the graphs using knowledge about the biology of the system. To build our map, we analyzed our dataset with samples as rows and genes as columns. We used the Variance Normalized Euclidean metric. We used the PCA coord1 and PCA coord2 lens as well as two data lenses: CG32373 and Fuca (Resolution = 19, Gain = 6, and equalized was used for all lenses).

Model Rationale and Methods
To evaluate the impact of within-host infection dynamics on host health curves, we created a compartmental model consisting of four ordinary differential equations for microbes (M), immune effectors (I), microbial damage effectors (D), and health (H).
Microbes grow in a sigmoid fashion at rate r as they deplete the host substrate, and immune effectors can either kill microbes outright at rate ε per effector or inhibit microbial growth at rate η. In these simulations, growth rate limitation was used as the main immune response. Two mechanisms contribute to immune effector production. Reflecting the sigmoid relationship between microbe density and AMP levels, immune effectors can increase in proportion to microbe density up to a maximal rate of τ Ã k M until effector levels saturate immune pathway machinery (k I ). The parameter σ I reflects the microbe concentration at the inflection point for immune induction. Immune effectors decay at rate μ. Microbes secrete virulence factors and toxins at a maximal rate φ Ã k D , modulated by a sigmoidal relationship between damage factor production and microbe density to reflect a process like quorum sensing controls on virulence factor production. As with immunity, there is a carrying capacity imposed on damage effectors. This could reflect either negative feedbacks on further effector production at high effector densities or a bottleneck limiting flux through the system. The parameter σ D reflects the microbe density that produces a half-maximal induction rate. Damage effectors degrade at rate ρ. Immunopathology and damage effectors deplete health (H) at rates γ and α, respectively, while hosts can recover health according to a sigmoidal relationship with current health (reflecting the difficulty of achieving recovery at low health due to organ failure and other catastrophe) at rate z. "Death" is called when health dips below 10% of maximum (k H ). All simulations were conducted within a parameter space outlined in S4 Table and run in Matlab (v. 7.11.0) using the ode45 solver.

Curve Fitting
The sigmoidal curves were fit using the four-parameter method in Prism. We tested several models (linear, exponential, logistic, and sigmoid) and picked the model that gave the best adjusted r 2 . We used the adjusted r 2 to account for overfitting. When the curve-fitting program suggested a clearly erroneous result, such as an extremely high top or low bottom, we fixed the top or bottom at the average level for the vigor or for the average of the last three points for the bottom.

RNAseq
Flies were injected with 50 nL of 0.01 OD600 (1000 CFUs) of L. monocytogenes resuspended in PBS. Flies were placed in vials containing dextrose fly media and incubated at 29°C. 3 DPI flies were flipped onto dextrose fly media containing 1 mg/mL ampicillin. Groups of 20 flies were homogenized in TRIzol for RNA extraction on the following DPI: 0, 1, 2, 3, 4, 6, 8, 10, 12, 14, 16, and 18. Additional flies were used to determine CFU and monitor survival. We tested w 1118 , CG2247, Pcmt, and RAL359 strains. As a control, additional sets of w 1118 were left uninfected and flipped onto dextrose fly media containing 1 mg/mL ampicillin 3 dpi, and samples for RNAseq were collected at the indicated time points. RNA was isolated using standard TRIzol preparation. Library preparation and sequencing was performed by the Duke Center for Genomic and Computational Biology. Briefly, quality of RNA was assessed by Bioanalyzer 2100. Polyadenylated RNA was enriched from total RNA. Barcoded TruSeq cDNA Libraries were constructed and quality was assessed by Qubit and Agilent Tapestation. 50 bp single end reads were obtained by sequencing the libraries on an Illumina HiSeq 2000 using a full flow cell. Reads were mapped using STAR, and RPKM for each gene was determined by Cufflinks.

Melanization Assay
Groups of 20 female flies were anesthetized with carbon dioxide and injected with 50 nl of 0.01 OD600 (~1,000 CFUs) of L. monocytogenes resuspended in PBS. After injection, flies were flipped onto dextrose fly media and incubated at 29°C. At 2 dpi, the percent melanized for each group was determined by anesthetizing the flies with carbon dioxide and scoring each abdomen for presence of melanization spots. For each genotype, eight groups of flies were scored.  Fig. 50 fuzzy c-means clusters. All differentially expressed genes were clustered into 50 fuzzy c-means clusters. Prior to clustering, data were standardized to have a mean value of zero and a standard deviation of one. The 50 fuzzy c-means groups were then manually curated into groups of similar expression patterns. Each line represents mean gene expression. Blue represents uninfected control. Red represents infected and no treatment. Green represents infected and ampicillin-treated. Cluster numbers highlighted in red contain mortality genes, and cluster numbers highlighted in green contain recovery genes.


extfieldsofstudy: Biology Medicine

text:
Tomato HAIRY MERISTEM genes are involved in meristem maintenance and compound leaf morphogenesis

Highlight Functional analysis of the tomato HAM genes suggest that they are involved in critical meristem functions and play homologous roles in compound leaf development.

Introduction
The shoot apical meristem (SAM) produces all of a plant's aboveground structures throughout its lifespan and it is also the site for the formation of new leaf and flower primordia. The SAM can be subdivided into a central zone (CZ) which consists of the stem cell niche at the SAM summit and the underlying organizing center (OC), a surrounding peripheral zone (PZ) of rapidly dividing cells, and an underlying rib zone (RZ). During organ formation, PZ and RZ cells are recruited into differentiating lateral organ primordia and the elongating stem, respectively (Steeves and Sussex, 1989). The ability of the SAM to maintain its size and continuously produce new organs depends on a small population of slowly dividing pluripotent stem cells in the CZ that constantly renews itself while providing daughter cells to the surrounding PZ and RZ (Tucker and Laux, 2007). WUSCHEL (WUS) is a central regulator of stem cell homeostasis and is specifically expressed in the meristem OC (Laux et al., 1996;Mayer et al., 1998). The floral meristem (FM) is a modified SAM, but, unlike the SAM, it is not indeterminate as its stem cells give rise to a precise number of floral organs and are abolished once the carpel primordia, the final organs to be made from the FM, form (Sablowski, 2007). Thus, a delicate balance between the two opposing activities of maintaining the undifferentiated stem cell population and promoting differentiation is precisely and carefully regulated in primary meristems.
The petunia hairy meristem mutant prematurely terminates both lateral organ and stem production due to precocious SAM termination. Terminating SAMs had a histological structure reminiscent of the radial pattern of the stem tissue subtending the wild-type SAM. This included the development of a differentiated epidermis with trichomes, hence the gene associated with the mutation designated HAIRY MERISTEM (HAM). Nevertheless, at termination, hairy meristem SAMs displayed normal PhWUS expression, suggesting that PhHAM is required for PhWUS responsiveness. PhHAM codes for a member of the GAI, RGA, and SCR (GRAS) family of transcription regulators (Bolle, 2004). It was concluded that PhHAM is required for meristem maintenance by protecting meristematic cells that initiate differentiation from developing by default into stem (Stuurman et al., 2002).
The Arabidopsis thaliana (Arabidopsis) genome encodes four HAM (AtHAM)/LOST MERISTEM (LOM) homologs, of which AtHAM1-AtHAM3 are targeted by miR171 (Llave et al., 2002;Wang et al., 2010). Atham1,2 and Atham1,2,3 mutants showed abrupt shoot termination and disturbed axillary bud formation, similar to AtMIR171c overexpressors and the Phham mutant (Stuurman et al., 2002;Engstrom et al., 2010;Schulze et al., 2010;Wang et al., 2010). The abrupt shoot termination was associated with an abnormal broader and flatter SAM that contained differentiating cells, indicating a reduced capacity to maintain itself and a loss of its polar organization (Schulze et al., 2010). Occasionally, Atham1,2,3 and Atham1,2 mutants developed bulges or symmetrical organs in their leaf axils that contained ectopic meristematic cell clusters (EMCCs). EMCCs were also evident in the differentiated young stem (Schulze et al., 2010). The blocked organogenesis despite ongoing cell division suggested that AtHAM1 and AtHAM2 are involved in the promotion of differentiation at the periphery of the SAM.
Recently, it was found that AtHAM1-AtHAM4 interact with WUS. The functionality of this interaction is supported by their co-localization in the SAM and the observation that the Atham1,2,3 loss-of-function mutant and the Atham4-pWUS:MIR171 mutant display defects similar to those observed in wus7 and wus1 mutants, respectively. In addition, AtHAM2 has been found to bind genomic regions similar to those reported to associate with WUS and to enhance WUS transcriptional activities. These data suggest that HAM proteins also function as cofactors for WUS-mediated stem cell niche maintenance (Zhou et al., 2015).
In contrast to the significant effect of ham mutations on meristems, in Petunia and pepper ham mutants, leaf morphology could not be distinguished from that of the wild type (Stuurman et al., 2002;David-Schwartz et al., 2013).
In Arabidopsis, loss of HAM functions caused narrower and more serrated rosette leaves and epinastic cauline leaves (Engstrom et al., 2010;Wang et al., 2010). These subtle leaf phenotypes suggest that as opposed to their essential function in meristems, HAM genes may play a minor role in the development of a simple leaf, which is a determinate organ.
Although determinate, specific regions of the leaf maintain transient intermediate growth. These include regions at the margin of the leaf primordium that may be patterned into a blade or any other marginal structures such as leaflets. Accordingly, these regions are active for short periods in the simple leaf primordium and remain active for longer periods in the compound leaf primordia that are elaborated by leaflets (Hagemann and Gleissberg, 1996). The mechanism that regulates leaf margin maintenance is not completely understood. To investigate whether HAM genes are involved in compound leaf development, we functionally analyzed SlHAM and SlHAM2 from the compound leafed Solanum lycopersicum (tomato) via reverse genetics. This analysis suggests that the function of SlHAMs in meristem maintenance was also conserved in tomato. However, in contrast to the minor role of HAM genes in simple leaves, our analysis uncovers an essential role for SlHAM genes in compound leaf morphogenesis.

Plant material and growth conditions
The tomato cv. M82 driver lines 35S:LhG4, FIL:LhG4, AP1:LhG4, and OP:AtCKX3 are described elsewhere (Shani et al., 2010;Hendelman et al., 2012). Tomato plants were grown under greenhouse conditions with temperatures ranging between 15 °C and 30 °C in a turf-peat mix with nutrients, using 4 liter pots. Germination and seedling growth took place in a growth chamber under a 16 h light, 8 h dark photoperiod (photosynthetic photon flux density: 50-70 μmol m −2 s −1 ) at a constant temperature of 24 °C. For crosses, the closed flowers of the corresponding homozygous responder line were emasculated by removal of the petals and stamens, and hand-pollinated with the pollen of the respective homozygous driver line.

Sly-miR171 cleavage site mapping
The sly-miR171 cleavage site was mapped in mixed mRNA extracted from leaves, flower buds, and flowers at anthesis using a modified procedure of RLM-RACE (RNA ligase-mediated rapid amplification of cDNA ends) as described before (Hendelman et al., 2012) with corresponding RACE and RACE-nested primers (for primer sequences, see Supplementary Table S1 at JXB online).

Plasmid construction
For sly-MIR171a and sly-MIR171b reporter constructs, the sequences flanking pre-miR171a (SL2.50ch02:46752720..46753003) and pre-miR71b (SL2.50ch02:44783423..44783704) were PCRamplified from tomato M82 genomic DNA with the primer pairs SlMIR171a_SalI_F or SlMIR171b_SalI_F and SlMIR171a_ BamHI_R or SlMIR171b_BamHI_R which contained SalI and BamHI sites at their 5' ends. After sequence verification, the amplified fragments were cloned into the XhoI/BamHI sites of the pOp-TATA-BJ36 shuttle vector between an OP array (Moore et al., 1998) and Agrobacterium tumefaciens octopine synthase terminator (OCS) to generate OP:SlMIR171a and OP:SlMIR171b constructs. The NotI fragments of these constructs were then mobilized into the binary vector pART27 to generate pART27-OP:SlMIR171a and pART27-OP:SlMIR171b.

Transformation of tomato plants
The binary vectors pART27-OP:SlMIR171a and pART27-OP:SlMIR171b were transformed into tomato cv. M82 by cocultivation of cotyledons derived from 14-day-old seedlings with A. tumefaciens strain GV3101 as described previously (Hendelman et al., 2013).
Total RNA extraction and sly-miR171 gel-blot analysis Total RNA was extracted from different tomato tissues with Bio-Tri RNA reagent (Bio-Lab, Jerusalem, Israel) according to the manufacturer's protocol. Sly-miR171 gel-blot analysis of total RNA was performed as described previously (Hendelman et al., 2013) using a complementary radiolabeled oligonucleotide as a probe.
Quantitative RT-PCR assays First-strand cDNA was synthesized from 2 µg of total RNA with a Maxima first strand cDNA synthesis kit (Thermo Scientific, Vilnius, Lithuania) following the manufacturer's instructions. A negative control (-RT) was used to ensure the absence of genomic DNA template in the samples. Three independent biological replicates were used for each sample, and quantification was performed in triplicate. PCR was performed in StepOnePlus (Thermo Scientific) following the manufacturer's instructions. PCR products were analyzed using StepOne software version 2.2.2 (Thermo Scientific). Primers were designed around the corresponding sly-miR171 complementary site. The relative expression levels were calculated using the two-standard-curve method normalized to SlTIP41 as the reference gene. Unless otherwise mentioned, statistically significant differences between samples were determined by Tukey-Kramer multiple-comparison test.

SEM
For scanning electron microscopy (SEM) analysis, different plant tissues were collected and placed in FAA (3.7% formaldehyde, 5% acetic acid, 50% EtOH, v/v/v) solution until use. Then the FAA was removed and the tissues were washed in an increasing gradient of ethanol (up to 100%). Fixed samples were critical point dried, mounted on a copper plate, and gold coated. Samples were viewed in a Jeol 5410 LV microscope (Jeol, Tokyo, Japan).

Histological analyses
For apices and leaves, tissues were collected and fixed with PFA [4% paraformaldehyde; 1× phosphate-buffered saline (PBS); 0.2% Tween-20 (v/v)] for 16 h. Then the samples were washed twice with 1× PBS and dehydrated in increasing concentrations of ethanol (10, 30, 50, 70, 80, 90, 95, and 100%). The samples were embedded with JB-4 (Electron Microscopy Sciences, Hatfield, PA, USA) according to the manufacturer's protocol with slight modifications: after dehydration, the samples were incubated in the infiltration solution and placed at 4 °C in the dark under vacuum for 6 h (the infiltration solution was replaced twice). Polymerization was performed under anaerobic conditions. For flower analyses, tissues were collected and fixed in FAA until use, then dehydrated in increasing concentrations of ethanol, cleared with K-clear (Kaltek, Padova, Italy), and embedded in paraffin. Microtome-cut sections (10 µm and 2.5 µm thick for plastic and wax, respectively) were spread on microscope slides and stained with 0.1% (w/v) Safranin followed by 0.2% (w/v) Fast green (flowers) or 0.1% (w/v) Toluidine blue O (apices and leaves). Slides were examined under bright-field using an Olympus DP73 microscope equipped with a digital camera.

The tomato miR171 (sly-miR171) guides the cleavage of three GRAS-like genes
Four members of the tomato miR171 family have been cloned to date, of which sly-miR171a and sly-miR171b are offset by three nucleotides relative to each other (Moxon et al., 2008). Using the psRNATarget web server (Dai and Zhao, 2011) and applying a mismatch score of ≤2.0, we predicted four putative sly-miR171 targets from the available tomato gene models (ITAG release 2.30) ( Supplementary Fig. S1A). Analysis of RNA from leaves, flower buds, and flowers by RLM-RACE confirmed sly-miR171-directed endonucleolytic cleavage of Solyc08g078800 (SlHAM), Solyc01g090950, and Solyc11g013150. Sequencing of the amplified products indicated that SlHAM and Solyc01g090950 are cleaved at positions that derived from targeting by sly-miR171a and sly-miR171b, whereas Solyc11g013150 was cleaved at positions that indicated targeting by sly-miR171b only ( Supplementary  Fig. S1A). In contrast, the Solyc02g085600 cleavage product was not recovered by us, consistent with its absence in published tomato degradome data (Karlova et al., 2013), suggesting that it is not subject to significant sly-miR171-mediated cleavage.
Analysis of the predicted protein sequences of SlHAM, Solyc01g090950, Solyc11g013150, and Solyc02g085600 revealed a similar overall architecture, including a variable N-terminus, which was much longer in SlHAM and Solyc01g090950, followed by a conserved GRAS domain and its characteristic VHIID, PFYRE, and SAW motifs (Pysh et al., 1999) (Supplementary Fig. S1B). Phylogenetic reconstruction of their GRAS domains ( Supplementary Fig. S1C) revealed that SlHAM, Solyc01g090950, and Solyc02g085600 belong to the HAM branch of GRAS (Bolle, 2004). The Solyc01g090950 protein was found to belong to a subclade that includes Arabidopsis HAM1, 2, 3, Petunia HAM, and SlHAM proteins (Stuurman et al., 2002;Engstrom et al., 2010) and, accordingly, was named SlHAM2. The Solyc02g085600 protein was found in the HAM subclade but shares the highest similarity with AtHAM4, which also does not undergo miR171-guided cleavage (Engstrom et al., 2010), and accordingly was named SlHAM4. The Solyc11g013150 protein was not found in the HAM subclade but in a subclade of GRAS that contains the Medicago truncatula GRASlike transcription factor NODULATION SIGNALING PATHWAY 2 (MtNSP2), which is targeted by mt-miR171h (Lauressergues et al., 2012). In line with this, amino acid sequence alignment between Solyc11g013150 and MtNSP2 showed significant overall homology (62%/77% identity/similarity) and therefore it was named SlNSP2L ( Supplementary  Fig. S1C).

Silencing of SlHAMs in meristems is associated with the formation of EMCCs
It was shown that HAMs from petunia, Arabidopsis, and pepper play a role in meristem maintenance (Stuurman et al., 2002;Schulze et al., 2010;David-Schwartz et al., 2013;Zhou et al., 2015). Consistent with this, published tomato RNA-sequencing data (Park et al., 2012) have indicated that SlHAM and SlHAM2 are relatively abundant in the SAM and FM. In contrast, in both meristems, SlNSP2L was expressed at ~11-to 15-fold lower levels and SlHAM4 expression was almost negligible (~50-to 150-fold lower), suggesting that they are less important for meristem function ( Supplementary Fig. S1D). To study the role of miR171targeted SlHAM and SlHAM2 (collectively referred to as SlHAMs) in tomato meristems, we utilized the OP/LhG4 transactivation system to silence them by overexpressing sly-miR171a, which directs the cleavage of both genes, or sly-miR171b that additionally cleaves SlNSP2L ( Supplementary  Fig. S2A). Following transformation into M82 tomato, 12 OP:MIR171a and five OP:MIR171b responder plants were obtained and their F 1 progeny characterized for corresponding mature miR171 overexpression following a cross with the 35S:LhG4 driver line. This analysis identified several responder lines that strongly expressed the respective sly-miR171 upon transactivation ( Supplementary Fig. S2B) from which we selected OP:MIR171a-4 (hereafter OP:MIR171a) and OP:MIR171b-20 (hereafter OP:MIR171b) for further analysis ( Supplementary Fig. S2B). Constitutively activated 35S>>MIR171a and 35S>>MIR171b F 1 progeny seedlings accumulated ~3-fold higher levels of corresponding sly-miR171 than the control (Fig. 1A) and that was accompanied by significantly lower levels of SlHAMs transcripts (Fig. 1B). A similar reduction in the levels of SlNSP2L was also observed, but only in 35S>>MIR171b seedlings (Fig.  1B), in accordance with its specific cleavage by sly-miR171b ( Supplementary Fig. S1A). Conversely, SlHAM4 showed an insignificant reduction in its transcript levels, thus confirming its resistance to sly-miR171a-and sly-miR171b-guided cleavage in vivo (Fig. 1B). In contrast to whole seedlings, examination of SlHAMs levels in shoot apices of seedlings showed more efficient silencing in 35S>>MIR171b than in 35S>>MIR171a. In addition, SlNSP2L silencing was weaker than that of SlHAMs in 35S>>MIR171b shoot apices. Whereas SlHAMs levels were down-regulated by ~55%, SlNSP2L levels were down-regulated by only 36% (Fig. 1C). Compared with the control (35S:LhG4), the 35S>>MIR171a seedlings displayed growth arrest after producing a few pairs of leaves, suggesting that SlHAMs are required for normal meristem function (Fig. 1D). Consistent with the stronger silencing of SlHAMs in 35S>>MIR171b seedlings, they displayed a similar, albeit more severe phenotype, and usually arrested growth after producing a single pair of abnormal leaves (Fig. 1D). In addition, they developed trichomes on the adaxial side of their arrested P2 leaf primordia, indicating precocious differentiation ( Fig. 1F-H). However, unlike the Petunia ham mutant (Stuurman et al., 2002), their SAMs did not display trichomes ( Fig. 1I-K). Moreover, in two-thirds of the 35S>>MIR171b seedlings, the shoot apex already appeared slightly swollen at 12 days after germination (DAG) and the swollen area enlarged over time to form a bulge between the cotyledons (Fig. 1E), which grew until it tore the surrounding tissue ( Fig. 2A, 25 DAG). A similar phenotype was occasionally observed in 35S>>MIR171a apices ( Supplementary Fig. S3). The etiology and ultrastructure of the bulge were investigated by histology of young 35S>>MIR171b seedlings and histone H4 distribution, as a marker for cell division (Brandstädter et al., 1993). In control tomato apices, high meristematic activity, as indicated by small densely Toluidine blue-stained cytoplasmic cells ( Fig.  2A, 35S:LhG4) and positive histone H4 signal (Fig. 2B sections 4-5, and C), was restricted to the SAM and the adaxial side of the leaf primordia. In contrast, in 35S>>MIR171b apices, ectopic clusters of small densely staining ( Fig. 2A) and histone H4-positive ( Fig. 2B sections 2-3, and D) cells were apparent between the cotyledons, indicating abnormal cell proliferation. With time, this led to the build up of a cell mass in the form of a bulge between the hypocotyl and the SAM ( Fig. 2A, 35S>>MIR171b). The apical part of the bulge consisted of small meristematic cell clusters that produced an inner mass of large and vacuolated cells, suggesting that as cells were displaced to the periphery they lost their meristematic nature without adopting a specific fate ( Fig. 2A, 25 DAG).
Since SlHAM genes are also relatively abundant in the FM ( Supplementary Fig. S1D), we examined their roles by reducing their expression levels specifically in the FM through crossing the OP:MIR171a and OP:MIR171b responder lines with the flower-specific AP1:LhG4 driver line. AP1:LhG4 drives expression throughout the FM and, when meristem activity is abolished, its expression is confined to sepals and petals (Hendelman et al., 2012). Quantitation of sly-miR171targeted genes in the young inflorescence with flower buds <1 mm in size revealed a significant ~65% reduction of SlHAM and SlHAM2 in AP1>>MIR171a and AP1>>MIR171b, whereas no significant reduction in SlNSP2L levels was detected in these plants, which is consistent with its less efficient silencing in shoot apices (Fig. 3A). Phenotypic analysis of developing floral buds showed a similar abnormal morphology in AP1>>MIR171a and AP1>>MIR171b. The observed phenotype, which was more pronounced in AP1>>MIR171b, was already visible in 2 mm buds that displayed swollen proximal ends compared with AP1:LhG4 control buds. Like the 35S>>MIR171b seedling bulge ( Fig. 2A), the swelling increased with time (Fig. 3B). In contrast to AP1>>MIR171b flowers, which senesce and could not bear fruit (Fig. 3B), the less severely affected AP1>>MIR171a buds reached anthesis and could set fruit. Moreover, the swollen area in the flower buds of AP1>>MIR171a plants continued to expand during fruit development, suggesting its indeterminate nature (Fig.  3C). Accordingly, histology of very young AP1>>MIR171a and AP1>>MIR171b flower buds revealed that bud swelling was caused by excess proliferation of cells above the receptacle, which was more dramatic in the AP1>>MIR171b buds ( Fig. 4A-C). In contrast to the characteristic large vacuolated cells of the receptacle, the excess growth contained a mixed cell population that included numerous foci of relatively small densely stained cells (Fig. 4D-F), reminiscent of the actively dividing cell foci observed in the 35S>>MIR171b apex ( Fig.  2A, B) and the EMCCs observed in the Atham1,2 mutant (Schulze et al., 2010). Similar analysis of AP1>>MIR171a and AP1>>MIR171b flowers at anthesis revealed a larger tumor-like cell mass intertwined with numerous vascular bundles, demonstrating that as the flower matures, the respective cells continue to proliferate (Fig. 4G-I). Taken together, our data indicate that depletion of SlHAMs from the SAM and FM causes the formation of EMCCs in their periphery.

SlHAMs are involved in compound leaf morphogenesis
We found that SlHAMs, which are relatively abundant in meristems, are also abundant in the compound leaf primordia, whereas SlNSP2L and SlHAM4, which were weakly expressed in the SAM and FM, maintained the same weak expression trend in the leaf primordia ( Supplementary Fig. S1D). This led us to ask whether SlHAMs have a similar function in compound leaves. To investigate this, the OP:MIR171b responder line was crossed with the FIL:LhG4 line driving transgene expression in lateral organ primordia but not in meristems  LhG4) and sly-miR171-overexpressing (35S>>MIR171a, 35S>>MIR171b) seedlings. Total RNA was extracted from pooled independent seedlings (n=6 per biological replicate). Sly-miR171 expression levels were determined relative to the control after normalization to the U6 snRNA and are indicated below each panel. (B, C) Quantification of sly-miR171 target mRNA levels in the RNA samples analyzed in A (B) and in apices of seedlings at 1.5 days after germination (DAG) (C) by qRT-PCR. Error bars indicate the SD over three biological replicates. The levels are expressed relative to the control, which was set to 1 ±SD. Different letters indicate statistically significant differences at P<0.01. (D) Seedlings from the indicated genotypes at 26 DAG. An arrow marks the bulged tissue. Scale bar=1 cm. (E) A close-up view of vegetative apices of 26 DAG seedlings from indicated genotypes. Two representative 35S>>MIR171b apices with different bulges are shown. C, cotyledon; L1, first leaf; L2, second leaf. Scale bar=1 cm. (F-K) Scanning electron micrographs of shoot apices of 11 DAG seedlings from the indicated genotypes. The first leaf pair (F, G) and advanced leaf primordia (F-H) were removed. Arrows indicate representative trichomes on the adaxial side of the leaf primordium. (I-K) Magnified view of the meristem epidermis. Scale bars (F-H)=100 µm; (I-K)=25 µm. (Shani et al., 2009). Quantitation of sly-miR171-targeted genes in FIL>>MIR171b leaf primordia confirmed a significant reduction in the expression levels of SlHAMs compared with that of the control (Supplementary Fig. S4A). In the early developing leaves of FIL>>MIR171b plants, we observed relatively subtle phenotypes. Compared with the control FIL:LhG4 line, these leaves were simpler and exhibited epinastic curling (Fig. 5A-D; Supplementary Fig. S4B). Histology of the fifth leaf indicated that the terminal leaflet blade is thicker than the control and contains more spongy mesophyll layers (Fig. 5E, F). In addition, both the leaflet and rachis displayed an atypical vasculature (Fig. 5E-H). Surprisingly, later developing leaves showed a dramatic increase in phenotype severity in FIL>>MIR171b plants. By leaf 10, the produced leaves were simpler, twisted, and their leaflet growth angle was distorted (Fig. 5I, J; Supplementary Fig. S4C). Later developing leaves were even more twisted and their leaflet growth angle was severely distorted. In addition, compared with the control, their rachis and attached petiolules were thicker and parts of their adaxial domain were abnormally brownish (Fig. 5K, L). As the leaves grew, a shiny light green tissue mass appeared on the adaxial side of the brownish rachis and petiolules (Fig. 5M,   N). SEM of a cross-section of the FIL>>MIR171b brownish rachis revealed abnormal morphology and a bulge composed of relatively small cells at its adaxial side (Fig. 5O-Q). Moreover, in contrast to the hairy epidermis of control rachis, the epidermis above the bulge completely lacked trichomes (Fig. 5Q). The lack of trichomes and the presence of small cells instead of large cells suggest an undifferentiated meristematic nature of the ectopic bulging tissue. Microscopic examination of the adaxial side of the FIL>>MIR171b rachis confirmed its abnormal morphology and revealed a mixture of vegetative and reproductive lateral organ primordia that may have differentiated from the meristematic tissues (Fig. 5R). These observations suggest that in the organogenic compound leaf rachis, as in meristems, silencing of SlHAMs provokes overproliferation of meristematic cells.

The EMCCs in the FIL>>MIR171b leaves are suppressed by expression of a CK-degrading enzyme
In tomato, elevation of cytokinin (CK) levels or response enhancement led to the proliferation of ectopic shoot and inflorescence meristems on the adaxial side of the compound leaf rachis which is highly reminiscent of the phenotype observed in FIL>>MIR171b later developed leaves (Janssen et al., 1998;Kim et al., 2003;Shani et al., 2010;Steiner et al., 2012). This suggested that an increased CK level or response may be the cause of the abnormal cell proliferation observed in later developed leaf rachises of FIL>>MIR171b plants. To test this idea, we reduced the CK levels in the FIL>>MIR171b leaves by expressing the Arabidopsis CK degradation gene CYTOKININ OXIDASE3 (AtCKX3) (Werner et al., 2003). To do this, OP:MIR171b was crossed to the previously characterized OP:AtCKX3 responder line, which expresses AtCKX3 upon transactivation (Shani et al., 2010), and their double responder line progeny were crossed to the homozygous FIL:LhG4 driver line. The leaf rachises of FIL>>AtCKX3 tomato plants displayed a wild-type phenotype (Fig. 6B,F). In contrast to later developed leaves of FIL>>MIR171b which develop ectopic organ primordia on the adaxial side of their rachises due to abnormal proliferation of meristematic tissue (Fig. 6G, K), the leaves of FIL>>MIR171b>>AtCKX3-coexpressing MIR171b and AtCKX3 under the FIL promoter displayed rachises that are morphologically identical to control and FIL>>AtCKX3 leaves (Fig. 6D, H). In addition, histology of the FIL>>MIR171b>>AtCKX3 leaf rachises confirmed that they possess wild-type morphology and are devoid of any abnormal clusters of densely stained meristematic cells, suggesting their suppression by a reduction of CK content (Fig. 6I-L). Still, similarly to FIL>>MIR17b leaves, FIL>>MIR171b>>AtCKX3 leaflets exhibited epinastic curling (Fig. 6D), suggesting that this phenotype is not affected by CK content reduction.

SlWUS is misexpressed in the leaf rachis, SAM, and FM silenced for SlHAMs
In Arabidopsis, CK is an important regulator of meristem establishment and maintenance by promoting WUS expression (Gordon et al., 2009;Kuroha et al., 2009;Bartrina et al., 2011). Thus, to understand further how CK is involved in the overproliferation of meristematic cells in the FIL>>MIR17b rachis, we compared SlWUS expression in the 10th leaf rachis, which does not harbor visible ectopic meristematic tissue, between control, FIL>>AtCKX3, FIL>>MIR17b, and FIL>>MIR171b>>AtCKX3. In control and FIL>>AtCKX3 rachises, we could not detect SlWUS by quantitative PCR (data not shown) or semi-quantitative RT-PCR, and similarly SlWUS was not detected by either in the rachises of FIL>>MIR171b>>AtCKX3 leaves, corroborating their wildtype phenotype. Conversely, SlWUS was readily detected in FIL>>MIR171b rachises (Fig. 7A). Moreover, in situ hybridization detected SlWUS in the FIL>>MIR171b leaf rachis, but not in control leaf rachis that expressed SlHAM (Fig. 7B-E). This suggested that SlHAMs silencing in the leaf rachis caused SlWUS misexpression probably via CK. To test if SlWUS misexpression also occurs upon SlHAMs silencing in meristems, we compared the spatial distribution of SlHAM and SlWUS between control and SlHAMs-silenced SAM and FM. In situ hybridization on control vegetative apices showed that SlHAM transcripts accumulate in a pattern similar to that of AtHAM1 (Schulze et al., 2010). SlHAM accumulated mainly in the PZs from the L2 layer downward and in the rib meristem, but it was absent from the CZ that hosts the OC and stem cell niche (Fig. 7F). Instead of the characteristic OC-confined expression of SlWUS in the SAM of control plants (Fig. 7G), in 35S>>MIR171b plants, the SlWUS expression domain expanded into the SAM periphery (Fig. 7I), probably causing the elevated SlWUS transcript levels detected by qRT-PCR (Fig. 7J). As in the SAM, non-overlapping expression domains of SlHAM and SlWUS were evident in the FM of control plants ( Supplementary Fig. S5C, E) while SlWUS expression domains expanded in the FM of FIL>>MIR171b plants ( Supplementary Fig. S5D, F). This suggests that similar to the leaf rachis, silencing of SlHAMs in the SAM and FM is correlated with SlWUS misexpression. In agreement with WUS being a positive regulator of stem cells and CLV3 being a stem cell marker (Fletcher et al., 1999), an ectopic SlCLV3 expression was observed in the leaf rachis, SAM, and FM silenced for SlHAMs (Supplementary Fig. S6). Thus, in the absence of SlHAMs, misexpression of SlWUS led to the formation of extra stem cells that most probably caused the EMCCs observed in the SAM, FM, and the leaf rachis. The same RNA was used as in Fig. 1C. Error bars indicate the SD over three biological replicates (n=20 per biological replicate). Different letters indicate statistically significant differences at P<0.01.

Discussion
The HAM genes encode GRAS transcriptional regulators that have been implicated in meristem maintenance but were found to play only minor roles in simple leaf development (Stuurman et al., 2002;Schulze et al., 2010;Engstrom et al., 2010;David-Schwartz et al., 2013;Zhou et al., 2015). In this study, we used global and tissue-specific trans-silencing to analyze SlHAMs functions in tomato that have compound leaves. A common phenotype of trans-silenced meristems and compound leaves was the formation of an ectopic cell mass as a result of overproliferation of meristematic cells.
Does SlNSP2L contribute to the observed overproliferation phenotypes? Despite SlNSP2L targeting by sly-miR171b, our analyses indicate that its silencing in vegetative and reproductive apices, which overexpress sly-miR171b, was weaker than that of SlHAMs (35S>>MIR171b; Fig. 1C) or not significant (AP1>>MIR171b; Fig. 3A), suggesting a minor or no contribution, respectively, to their abnormal phenotypes. The reduced silencing efficiency of SlNSP2L compared with SlHAMs may be caused due to the presence of an additional mismatch with in the seed sequence of sly-miR171b ( Supplementary Fig. S1A). Moreover, cell proliferation phenotypes were also observed in 35S>>MIR171a and in AP1>>MIR171a apices, which contained wild-type SlNSP2L levels (Figs 1C, 3A), indicating that its presence cannot rescue these phenotypes, probably because its function is not redundant with that of the silenced SlHAMs. This conclusion is also supported by the relatively weak expression of SlNSP2L in the SAM, FM, and leaf primordium ( Supplementary Fig. S1D) and the absence of abnormal cell proliferation phenotypes in known nsp2 mutants (Kaló et al., 2005;Liu et al., 2011). In contrast to SlNSP2L, SlHAMs were abundant in the SAM and FM, and SlHAM transcript accumulated mainly in the PZ and was highly reduced or excluded from the CZ. This pattern of transcript accumulation in meristems was very similar to that of HAM genes from Petunia, pepper, and Arabidopsis, suggesting a conserved function for SlHAMs in tomato meristems (Stuurman et al., 2002;Schulze et al., 2010;David-Schwartz et al., 2013).
The ability of the SAM to produce new organs requires a delicate balance between maintenance of the indeterminacy potential and organogenesis. All known ham mutants ceased organ formation after producing a variable number of leaves as a result of precocious termination of the SAM (Stuurman et al., 2002;Engstrom et al., 2010;Schulze et al., 2010;David-Schwartz et al., 2013). This indicated that HAM genes are required to maintain SAM activity, and several roles have been proposed. In the Petunia ham mutant, the SAM differentiated into a stem two weeks after termination, suggesting that PhHAM is required for meristem maintenance by protecting meristematic cells that initiate differentiation from developing by default into stem (Stuurman et al., 2002). In Arabidopsis, AtHAMs were recently demonstrated to be important for stem cell production by physically interacting with WUS and facilitating its related transcriptional activities. Accordingly, similar to wus mutants, Atham mutants show an arrested growth phenotype probably due to precocious termination of their stem cell-depleted SAMs. In tomato, physical interaction between SlHAM and SlWUS suggests a similar function in stem cell homeostasis for the former (Zhou et al., 2015). This is supported by early vegetative growth arrest of SlHAMs-silenced seedlings and the direct correlation between the silencing efficiency and growth arrest severity. While the 35S>>MIR171b apices arrested just after producing two leaves, the 35S>>MIR171a apices, which contained higher SlHAMs levels, produced more leaves (Fig.  1C, D). However, since SlHAMs-silenced seedlings contained EMCCs already at 1.5 DAG, it is hard to determine the cause and effect of their growth arrest.
In shoot meristems, AtHAM1 and AtHAM2 have also been suggested to promote cellular differentiation of CZ descendants. This is supported by their preferred accumulation in the SAM PZ, the blocked organogenesis, and the development of EMCCs in the axils and young stems of corresponding mutants. Moreover, the EMCCs accumulated the transcripts of WUS and the stem cell marker CLV3, and occasionally they proliferated to produce an axillary bulge instead of an organ (Schulze et al., 2010). In tomato, SlHAMs accumulate in the meristem PZ, and their silencing in the SAM and FM was also associated with EMCCs that produced bulges beneath the corresponding meristem. In addition, formation of EMCCs was associated with expansion of SlWUS and SlCLV3 expression domains. Taken together, these observations support a similar role for SlHAMs as promoters of cellular differentiation in meristems.
Interestingly, we also observed the formation of EMCCs associated with SlWUS and SlCLV3 misexpression in the adaxial side of the compound leaf rachis, the site from which leaflet primordia emerge, and meristematic potential persists due to the activity of a pool of undifferentiated cells at the leaf margin (Hagemann and Gleissberg, 1996;Kim et al., 2003). It was suggested that once lobes and leaflets are initiated in the tomato compound leaf, this meristematic activity becomes restricted to those structures and down-regulated along the rachis (Kang and Sinha, 2010). This raises the possibility that as in meristems, SlHAMs may promote the differentiation of the rachis cells in order to prevent excess cell proliferation that could lead to the development of ectopic lateral organs. In tomato, leaf complexity increases acropetally, suggesting that the early developed leaf rachis maintains its meristematic potential for a shorter period than the later developed one (Busch et al., 2011). This provides a possible explanation for why EMCCs were never observed in the rachises of early developing leaves. In these rachises, the meristematic activity probably arrested before their cells could proliferate to significant numbers to generate an EMCC. Accordingly, meristematic cell proliferation has never been reported in the previously characterized ham mutants (Stuurman et al., 2002;Engstrom et al., 2010;Wang et al., 2010;David-Schwartz et al., 2013), probably due to the relatively rapid differentiation of their simple leaves. Instead, their simple leaves showed no phenotype (Stuurman et al., 2002;David-Schwartz et al., 2013) or relatively subtle phenotypes, including an increase in cell number along the adaxial-abaxial axis of the leaf and epinastic curling (Engstrom et al., 2010;Wang et al., 2010). The latter are reminiscent of the phenotypes observed in FIL>>MIR171b leaflets, suggesting that HAM genes also play a role in blade development.
We have shown that leaf-specific expression of AtCKX3, a CK-degrading enzyme, did not rescue blade development but suppressed the formation of EMCCs in FIL>>MIR171b leaf rachises, suggesting that EMCCs are CK dependent. Consistent with that, in tomato, increased CK levels or signaling cause the formation of ectopic meristems on the adaxial side of the leaf rachis and leaflet petiolules (Shani et al., 2010;Steiner et al., 2012). Since CK was shown to promote WUS expression in meristems (Gordon et al., 2009), the expansion of the WUS domain in meristems silenced for SlHAMs may imply that CK status increased in these meristems. This may hold for SlHAM-silenced leaf rachises since SlWUS was also misexpressed in them. Taken together, our results raise the possibility that by negative regulation of CK levels in the SAM, FM, and the compound leaf margin, either directly or indirectly, SlHAMs may promote cellular differentiation in them. Further research is required to confirm this hypothesis.

Supplementary data
Supplementary data are available at JXB online. Table S1. Primers and probes used in this study. Figure S1. Characterization of sly-miR171 predicted target genes. Figure S2. Generation and screening of transgenic tomato responder lines. Figure S3. Phenotype of 35S>>MIR171a seedlings 25 DAG. Figure S4. Characterization of FIL>>MIR171b leaves. Figure S5. In situ hybridization of SlHAM and SlWUS in FIL>>MIR171b floral buds. Figure S6. In situ hybridization of SlCLV3 in SlHAMsilenced plants.