The size and shape of organs is species specific, and even in species in which organ size is strongly influenced by environmental cues, such as nutrition or temperature, it follows defined rules. Therefore, mechanisms must exist to ensure a tight control of organ size within a given species, while being flexible enough to allow for the evolution of different organ sizes in different species. We combined computational modeling and quantitative measurements to analyze growth control in the Drosophila eye disc. We find that the area growth rate declines inversely proportional to the increasing total eye disc area. We identify two growth laws that are consistent with the growth data and that would explain the extraordinary robustness and evolutionary plasticity of the growth process and thus of the final adult eye size. These two growth laws correspond to very different control mechanisms and we discuss how each of these laws constrains the set of candidate biological mechanisms for growth control in the Drosophila eye disc.
The size and shape of organs is species specific to an extent that they represent major taxonomical traits. Even in species where organ size is strongly influenced by environmental cues, such as nutrition or temperature, it follows defined rules, called reaction norms (Moczek et al., 2011). Therefore, mechanisms must exist to ensure a tight control of organ size within a given species, while at the same time they should be flexible enough to allow for the evolution of organs that are different sizes in different species. These mechanisms operate during the development of the organism, i.e. they are developmental mechanisms. Still, the nature of these mechanisms is far from understood.
Much of what is currently known about how final organ size is controlled has been learnt from insects and particularly from studies in Drosophila (Shingleton, 2010). The external organs of the Drosophila adult (or imago) develop from larval primordia formed by epithelial sacs, called imaginal discs (Cohen, 1993). Imaginal discs have an interesting property: they stop growth when their normal size has been attained even if given extra developmental time (Bryant and Levinson, 1985; Bryant and Simpson, 1984; Garcia-Bellido, 1965; Martín and Morata, 2006). Therefore, imaginal discs have an autonomous mechanism that surveys organ size and terminates growth once a target size has been reached. Further evidence suggests that intrinsic mechanisms of growth termination also exist in vertebrates (Dittmer et al., 1974; Metcalf, 1964; Twitty and Schwind, 1931). A number of experiments rule out a simple cell-counting mechanism, by which cell proliferation would stop after a prefixed number of cell divisions. Thus, increasing or decreasing cell proliferation by modifying cyclin/cyclin-dependent kinase activity results in smaller or larger cells, respectively, but normal-sized discs (Neufeld et al., 1998; Weigmann et al., 1997).
In imaginal discs, growth rates are controlled by patterning signals (Shingleton, 2010). In the wing disc – where most studies have been carried out – the pace of growth is modulated by the combined action of at least two signaling pathways: that of the BMP2-like morphogen Decapentaplegic (Dpp), which is secreted from a stripe in the middle of the disc, and signaling of the protocadherin Fat, so that, together, uniform growth rates are attained throughout the disc (Schwank et al., 2012). Two sets of models have been proposed to explain how organ growth terminates. One of them posits that the dynamics of Dpp signaling itself is responsible for growth termination (Wartlick et al., 2011). In this model, cells measure the relative increase in the Dpp signal during development. Cell division is triggered every time a fixed relative change in the local Dpp concentration has been reached and uniform tissue growth is achieved despite localized Dpp production and the resulting Dpp gradient, because the Dpp gradient scales with organ growth. According to this model, the observed decline in the cell division rate with time is explained by the almost linear increase in the Dpp concentration with time, such that it takes progressively longer to reach that fixed relative concentration increase. Wing disc growth has, however, been reported to be normal in clones that lack the only Dpp transducer Mad and its downstream target Brk (Schwank et al., 2012), and recent studies show that growth depends on Dpp only during the first half of larval development (Akiyama and Gibson, 2015) and that lateral cells divide at normal speed if the Dpp gradient is abolished (Harmansa et al., 2015).
Another set of models includes the potential role of mechanical feedback during growth. According to these models, growth is limited by the mounting pressure in the center of the growing disc domain and stimulated by the increasing strain in the outer parts of the domain, thus enabling uniform growth and growth termination (Aegerter-Wilmsen et al., 2007, 2012; Hufnagel et al., 2007; Shraiman, 2005). When linked with a signaling model for the disc's key patterning signaling system [Dpp, Wingless (Wg), Notch] as well as for Vestigal, Four-jointed, Dachsous and components of the Hippo pathway, the mechanical feedback model can reproduce many additional experimental observations beyond uniform growth and growth termination, but still fails to yield the experimentally observed growth kinetics and final disc size (Aegerter-Wilmsen et al., 2012). Growth and tension anisotropies have been associated with the regulation of the direction of cell division, thus linking growth to organ shape (Aegerter-Wilmsen et al., 2012; Baena-López et al., 2005; Benham-Pyle et al., 2015; González-Gaitán et al., 1994).
Complementary to the wing imaginal disc, the eye imaginal disc offers an attractive system to study growth control (Fig. 1A). In the eye disc, dpp is also expressed in a stripe and plays an important role in eye patterning (Chanut and Heberlein, 1997; Pignoni and Zipursky, 1997). However, the relative position of the Dpp-producing stripe is not stationary as in the wing disc. Instead, the Dpp stripe is initiated at the posterior pole, and then sweeps across the disc like a wave towards its anterior pole (Heberlein et al., 1993). Moreover, proliferation in the eye disc is not uniform: the Dpp stripe coincides with an indentation of the epithelium called the morphogenetic furrow (MF). Anterior to the MF, progenitor cells proliferate, with a peak in the first mitotic wave, followed by a synchronous entry in the G1 phase of the cell cycle at the MF. Posterior to the MF, cells differentiate, either directly or after a terminal mitotic round (Wolff and Ready, 1993). Therefore, most proliferation occurs anterior to the MF. Quantitative studies of eye disc growth revealed a correlation between cell proliferation rates and the relative change in Dpp signaling, similar to the previously reported correlation for the wing disc (Wartlick et al., 2014). In the eye disc, a local increase in Dpp concentration is the result of movement of the Dpp gradient as the MF progresses. However, much as in the wing disc (Schwank et al., 2012), clones that lack both the Dpp signal transducer Mad and its downstream target and pathway repressor Brk grow at the same rate as wild-type cells (Wartlick et al., 2014). Therefore, the mechanism that controls growth termination still remains elusive.
To define alternative candidate mechanisms for growth control we studied the eye disc growth process quantitatively. To this end, we measured the anterior and posterior areas and lengths at different developmental stages. We then analyzed the growth process computationally. We show that two distinct growth laws are consistent with the growth data, and discuss how these growth laws constrain the possible biological mechanism underlying the remarkable robustness of the final adult eye size in Drosophila.
A quantitative analysis of eye disc growth
To characterize the eye disc growth process, we collected eye discs at different stages of development from early to late third instar (∼76-130 h after egg deposition), stained the apical surface using an antibody against the apical marker atypical Protein kinase C (aPKC) (Suzuki et al., 2001; Wodarz et al., 2000) and measured the anterior (A) and posterior (P) areas, which are separated by the MF, where Dpp is produced (Fig. 1A). Measurements were made both in 3D (Fig. 1B) and on 2D projections (Fig. 1C) generated from confocal z-stacks (for details, see Materials and Methods). We noted that 2D and 3D values were perfectly correlated for all considered measures, i.e. for the posterior length, LP (Fig. 1D) as well as for the total area T (Fig. 1E), the posterior area P (Fig. 1F, blue) and the anterior area A (Fig. 1F, gray). We could therefore use this correlation to extrapolate our larger 2D datasets to 3D.
During eye disc development, growth is terminated when the MF reaches the anterior side by exhausting the pool of progenitor cells. Whereas the speed of the MF is constant (after an initial lag phase; Wartlick et al., 2014), the anterior-posterior (AP) axis of the eye disc initially expands rapidly, but subsequently the expansion speed declines (Fig. 1G, Fig. S1) and the MF can reach the anterior pole and terminate growth. This is also reflected in the relative sizes of the anterior and posterior sides: initially the anterior area grows faster than the posterior area, but as the expansion of the AP axis slows down, the anterior area shrinks as the posterior area further expands (Fig. 1H). In principle, the slowdown in eye disc growth could result from a shrinkage of the proliferating anterior area due to the progression of the MF or from a decline in the area growth rate. We therefore sought to determine the area growth rate from our dataset.
The growth rate declines over developmental time
To determine the area growth rate k(LP) from the data, we determined and A(LP) at different posterior lengths LP. was obtained by plotting the total area T versus LP (Fig. 2A) and by subsequently determining the slope of the fit. The data points were not perfectly approximated by a linear fit (Fig. 2A, black line), and we therefore also considered three spline fits (Fig. 2A, gray lines). When we plotted the resulting k(LP) versus LP, we saw that the spline fits yielded slightly higher k(LP) values at the beginning and slightly lower k(LP) values at the end (Fig. 2B). The confidence in the estimated k(LP) is thus lowest for large LP values, for which we have the least data. In what follows, we will use the average of all the fits.
The eye disc shape can be approximated by an ellipse
We confirmed the applicability of the ellipse description of the eye disc (Fig. 3A) by comparing the measured maximal AP length of the eye disc LAP (Fig. 3B, green) with the one that we obtained from Eqn 8 (Fig. 3B, gray) with the measured LP (Fig. 1D), the measured T (Fig. 1E) and the measured A (Fig. 1F). The minor axis of length a then follows as a=LAP/2, and the major axis of length b as b=T/πa (Fig. 3C). Interestingly, the ratio of major to minor axis, σ=b/a, drops rapidly to about 2.3 (Fig. 3D), mainly because of an expansion of the minor axis a (Fig. 3C). After this initial phase, the eye disc grows isotropically and maintains an elliptic shape with a roughly constant ratio between the major and the minor axis (Fig. 3C,D).
Quantitative data-based evaluation of the growth models
We used the ellipse description as given by Eqns 8 and 9 to test the different growth rate models with our area measurements. To this end, we solved the dynamic model of the area increase (Eqn 1) with either a constant growth rate (CST) (Eqn 2), with a growth rate that decays according to a power law (POW) (Eqn 3), with a growth rate that declines with the increasing eye disc area (P_A) (Eqn 4) or with a growth law that declines exponentially (EXP) for k1=0 (Eqn 5) or with k1>0 (E+C) (Eqn 6). The anterior area A was determined from T and LAP according to the equation for ellipse segments (Eqn 8).
We optimized the parameter values for each growth model such that the predicted total area (Eqn 1, Fig. 4A), posterior area (Eqn 9, Fig. 4B) and anterior area (Eqn 8, Fig. 4C) would best fit the measurements (Table 1, Fig. S2). We notice that for the best-fitting parameter sets, all models with a declining growth rate also fitted the measured growth rate k(LP) (Fig. 4D, colored lines). A constant growth rate (Fig. 4, gray lines), however, fitted the data poorly and can therefore be ruled out and will not be further considered in what follows. Intriguingly, the models captured the experimentally observed decline in the ratio of major to minor axis, σ=b/a (Fig. 4E, colored lines), even though the parameter values were not optimized for these datasets.
An additional increase in the posterior area as a result of proliferation and cell differentiation, as can be incorporated by setting ɸ>1 in Eqn 1, does not affect our overall conclusion (Fig. S3). Thus, all models with a declining growth rate still fit the measured data very well and capture the measured growth rate. The growth model with an area-dependent growth rate is particularly robust to such changes, while the other growth models yield a larger final eye disc area as ɸ is increased (Fig. 4F). For all models, the match between measured and predicted ellipse shape σ=b/a worsens as ɸ is increased (Fig. 4E inset; Fig. S3, last column).
In summary, based on the analysis of the growth rate kinetics, we can reject a constant growth rate. The other growth laws fit the data similarly well. We therefore sought additional datasets to evaluate the different candidate growth laws.
Robustness of eye size
To evaluate the different growth laws, we next asked how each of the different growth laws would reproduce the naturally occurring size variations under the different conditions. Noise will introduce non-correlated changes into the parameter values, whereas changes in external conditions, e.g. temperature, can introduce correlated changes in the rate constants.
We first studied the impact of correlated changes in the kinetic parameter values on the final disc size and the time to termination. Environmental variations, e.g. temperature and nutrition, are well known to affect the final size of organs, but the impact is different in different organs. Thus, temperature changes affect the wing disc size substantially (3% per degree change), while eye disc size is considerably less affected (0.9% per degree change) (Azevedo et al., 2002). Interestingly, we found that for all growth laws, except for the power-law decay, the final disc size and the final anterior-posterior length (and thus the shape of the disc) did not change in response to parallel changes in the kinetic rate constants (v, k0, δ) (Fig. 5A,B). In case of a power-law decay in the growth rate, such robustness is observed only as long as the exponent δ (Eqn 4) is kept constant when the other parameters are changed (Fig. 5A,B, cyan symbols). Given the observed robustness of eye disc size to temperature changes, we can conclude that in case of a power-law growth law, the exponent δ must not be affected by external processes. This independence of δ would naturally be the case for the area-dependent growth rate, which is indeed robust to correlated changes (Fig. 5A,B, blue symbols).
Although the final area and AP length do not change in response to parallel changes in the parameter values (Fig. 5A,B), the time to termination changes in all cases, because a change in the MF speed, v, implies an opposite change in the time to termination (tterm∼1/v). In Drosophila, the speed of growth processes can indeed differ strongly as a result of changes in external conditions, but the final size of organs remains unaltered. Thus, when larvae develop at lower temperatures, eye discs attain about the same final size, whereas developmental time almost doubles (Azevedo et al., 2002).
To quantitatively compare the different growth models and to determine to what extent the decline in the growth rate would be adjusted when growth is slowed down, we used quantitative data for the growth process of grafted eye discs (Garcia-Bellido, 1965). Since no larval data were reported as a control in the grafting studies, we used our measurements of the eye discs as reference; we note that the final eye disc sizes were very similar for the larval (Fig. 5C,E) and grafting data (Fig. 5D,F). To compare our kinetic measurements to the grafting experiments, we had to convert the posterior lengths in our measurements to developmental time by using the measured MF speed. The average posterior length has been reported to increase by 3.1 µm/h (Wartlick et al., 2014). The average posterior length and the maximal posterior length LP are linearly related in our measurements (Fig. S4) and accordingly, the measured MF speed can be expressed in terms of LP as 3.4 µm/h.
Strikingly, the models with the declining growth laws not only maintained the same final size at a lower growth rate, but also reproduced the growth kinetics of the grafted eye discs very well (Fig. 5D,F). As before, the exponent δ had to be kept constant in the power-law growth law (Fig. 5F), and accordingly the area-dependent growth law fitted the grafting data very well (Fig. 5F, blue line).
Importantly, to achieve the slower, but size-preserving growth kinetics with the exponentially declining growth rate, the two rates, k0 and δ (Eqn 6), had to be changed in parallel in the model (to 0.15-fold the standard value). A mechanism would thus have to be in place that changes the maximal growth rate k0 and the rate at which the growth rate declines (δ) in parallel. In case of the area-dependent growth rate the decline in the growth rate depended directly on the eye disc area and was therefore self-adjusting when k0 was reduced. We conclude that an area-dependent growth law offers the most parsimonious mechanism to achieve the size-preserving slower growth kinetics in the grafts.
Variability of eye size between flies
We next wondered how sensitive our models would be to uncorrelated changes in the parameter values, as may arise because of noise. To this end, we changed the parameter values, k0, δ, σ(0), v as well as the initial sizes T(0) and LP(0) independently by 5% and recorded the impact on the final size T (Fig. 6A) and on the final length of the AP axis LAP (Fig. 6B), which is linearly related to the time to termination tterm. We noticed that the greatest robustness to uncorrelated parameter changes is obtained with an exponentially decaying growth rate (Fig. 6A,B, orange symbols) and an area-dependent growth rate (Fig. 6A,B, blue symbols), whereas a power-law growth rate results in the greatest sensitivity (Fig. 6A,B, cyan symbols).
The variation between left and right adult eyes is much smaller than the variation between eyes of different flies (Fig. 6C). In particular, even though eye size differs substantially between females and males, the sizes of left and right eyes are strongly correlated (R2=0.78; Fig. 6C), while the correlation between head width and eye size is weak (R2=0.29; Fig. 6D). This shows that a mechanism is in place that permits flies to coordinate the growth of their left and right eye discs independently of the final eye size and of the overall growth of the body. The coordination between the size of the left and right wing has recently been shown to be achieved by a Dilp8-dependent developmental delay (Colombani et al., 2012; Garelli et al., 2012). However, it has remained unclear how the discs know their final target size. Much as in case of the grafting experiments (Fig. 5), an area-dependent mechanism would naturally achieve this as long as the MF progression is delayed when the growth process becomes delayed, and the starting size and growth rate are similar for the two sibling eye discs.
Growth kinetics in differently sized eye discs
On the basis of the analysis of the growth processes in Drosophila eye discs, we found that the area-dependent growth rate and the exponentially declining growth rate are both valid candidate growth laws. We reasoned that the two growth laws could be best distinguished by comparing strains with different eye sizes because only for the area-dependent growth law does the decline of the growth rate depend on the area. With an exponential growth rate, the decline depends only on the developmental time passed.
We noticed that the eyes of the wild-type strain Oregon-R were significantly smaller than those of the GMR-GAL4 strain that we had used so far (Fig. 7A). We therefore also determined the growth kinetics for the eye discs of Oregon-R female larvae (Fig. 7B-D, red). As expected, the total area T of the Oregon-R eye discs increased more slowly than in the GMR-GAL4 strain (Fig. 7B), with smaller anterior areas at all times (Fig. 7C), whereas the eye disc shape was very similar at the same posterior length LP (Fig. 7D).
Given the differences in the growth rates and eye disc sizes, we could use the data to evaluate the plausibility of the two candidate growth rates. To this end, we determined the slopes in the diagnostic ln(k) versus LP plot (Fig. 7E) and in the ln(k) versus ln(T) plot (Fig. 7F) for the two strains. To support either growth law, the plots should be fitted well by straight lines. In case of the area-dependent growth rate, the slope of the straight line in the ln(k) versus ln(T) plot needs to be minus one, as is indeed the case, both for the GMR-GAL4 strain (black) and for the Oregon-R strain (red), even though the two strains both differed in total size over developmental time (Fig. 7B). The decline of the growth rate k(LP) is nonetheless in both cases inversely proportional to the total eye disc area T. This strongly supports an area-dependent growth rate. However, we note that the data in the diagnostic ln(k) versus LP plot (Fig. 7E) is also fitted rather well by a line. The available data thus remain insufficient to discriminate between the two growth laws.
We used Drosophila eye development as a model to study the mechanism of growth and organ size control. By combining quantitative measurements with a mathematical analysis, we found that the area growth rate declines with developmental progress in a way such that the area growth rate is inversely proportional to the total eye disc area, which was also true in a strain with smaller, more slowly growing eye discs (Figs 2 and 7). Strikingly, we found that two fundamentally different growth laws – an area-dependent growth law and an exponentially decaying growth rate – fit the eye disc growth data well.
However, the molecular mechanisms that would either lead to an area-dependent growth rate or an exponentially declining growth rate are very different. An exponentially declining growth rate could arise if a growth factor declined at a constant rate, δ, over time (Eqn 6). Such a factor could be intracellular and non-diffusive. Importantly, in the more slowly growing grafted eye discs, the δ of such a factor would have to scale directly with the slower developmental rate, potentially by being linked to the position of the MF, LP, rather than to absolute developmental time. In the case of an area-dependent growth rate, the grafting results can be explained without any such additional scaling condition. An inverse relationship between the growth rate and the total area could arise if a long-lived, diffusible, extracellular growth factor was diluted as the organ grows. The factor would need to be extracellular and diffusible because the growth rate declines uniformly, whereas cell division patterns are non-uniform in the eye disc (Wolff and Ready, 1991). Moreover, the cellular response to such a factor would need to be linear over the diluted concentration range (∼5-fold), and the dispersal of the factor would need to be limited to an area close to the apical cell membrane such that its dilution was relative to the area. Finally, the factor would need to be sufficiently long-lived to be diluted rather than degraded. This is a strong constraint, in particular in case of the grafting experiments where the eye disc grows over 2 weeks (Garcia-Bellido, 1965). However, we should note that in case of an exponentially declining growth factor the restriction would even be stronger in that the half-life of the candidate growth factor would not only have to be very long in grafted eye discs, but this half-life would have to be increased by the same factor by which development is delayed. Similarly, in case of a mechanical feedback model, it remains unclear whether the relaxation rate due to cell rearrangements (Shraiman, 2005) would be sufficiently slow to tolerate the hugely different developmental periods.
Declining growth rates are observed throughout development and aging, but, as a result of the moving MF, this would, in principle, not be necessary to terminate growth in the eye disc. However, in the eye disc, growth termination and the resulting eye size are much more robust with a declining growth rate than with a constant growth rate. In particular, an area-dependent growth rate maintains the final eye size also at vastly different developmental speeds, as has been observed previously at different culture temperatures and in grafting experiments. This independence on the developmental speed would also explain how a delay in metamorphosis, as mediated by the peptide Dilp8 in the wing disc (Colombani et al., 2012; Garelli et al., 2012) can reduce, rather than enhance, the size difference between left and right sibling discs (Fig. 6D). The final eye size can nonetheless be altered either via changes in the starting size, in the shape of the eye disc or by altering the area growth rate relative to the speed of the MF (Fig. 7). These aspects must thus be controlled globally to ensure a low variability between sibling discs.
In summary, we conclude that the area growth rate in the eye disc declines inversely proportional to the total area increase. Our failure to discriminate between fundamentally different growth laws with the help of our quantitative growth data emphasizes that the reproduction of a measured growth kinetic alone is insufficient evidence for the biological relevance of a particular growth mechanism. Further work will be required to define a mechanistic link between the area increase and the decline in the growth rate in the eye disc.
MATERIALS AND METHODS
The Drosophila strain used throughout was GMR-GAL4, a strain that is frequently used to target gene function modifications in the developing eye (Freeman, 1996) and which exhibits normal eyes. Oregon-R is a wild-type strain. All flies were raised on standard medium at 25°C. For all analysis, only female larvae were used if not indicated otherwise.
Antibody staining, fixation and imaging
Eye imaginal discs were dissected from larvae at different developmental points and fixed according to standard protocols (Casares and Mann, 2000). Rabbit anti-aPKC (Abcam, AB5813, 1:500) was used as primary antibody. Anti-rabbit Alexa Fluor 555 from Molecular Probes was used as fluorescently labeled secondary antibody (1:400). Stained discs were mounted in 80% glycerol in PBS, with DAPI (1:10,000) to counterstain nuclei, using three sticky tape layers as spacers between slide and the coverslip. Imaging of samples was carried out on a Leica TCS SPE microscope.
To measure the eye discs in 3D, we first reconstructed the 3D apical surface of the developing eye disc. To this end, the apical membrane was segmented using thresholding based on the aPKC antibody staining and subsequent manual corrections until a smooth and thin surface was achieved. Neighboring membranes (the apposing peripodial epithelium), as well as parts belonging to the antenna were removed manually. Surface reconstruction and measurements of 3D areas were done using Amira (FEI software). To measure the geometrical properties in 2D, maximum intensity projections were generated in ImageJ (Abràmoff et al., 2004) and the areas as well as the posterior length (Fig. 1B,C, yellow line) were also measured using ImageJ. The ventral eye disc folds onto itself forming a ventral ‘flap’, which can be visualized in 2D projections of the eye disc. We measured ventral flap areas separately and included these in the 2D measurements.
Eye sizes of adult flies were measured using ImageJ. The heads of adult flies were mounted in Hoyer's: lactic acid (1:1) mounting medium and cleared by heating at 70°C overnight. Images were taken focusing on the front and back planes of the eye to account for all of the eye area. Front and back areas were summed up to give the total eye area.
All models were simulated as ordinary differential equations (ODEs) of the form: , with k(LP) as summarized in Table 1, using a forward Euler scheme. In each step, the best-fitting ellipse shape σ was determined given total area T, posterior area P, and posterior length LP. Parameters for the different growth laws were estimated by fitting the models to the measured total, anterior and posterior areas. Parameter estimation was done using a trust-region-reflective algorithm as implemented in the lsqnonlin function in MATLAB R2014b (Coleman and Li, 1996). All simulations were done using MATLAB R2014b or the free software environment R version 3.1.1 (R Development Core Team). All statistical analyses were done using R.
We thank H. Sun (Academia Sinica, Taipei), E. Bach (Department of Biochemistry and Molecular Pharmacology, NYU, New York) and J. C. Hombría (CABD, Seville) for Drosophila stocks and the CABD Advanced Light Microscopy Facility and A. Iannini for technical support. We also thank Renato Paro and colleagues in our laboratories for discussions.
D.I. and F.C. conceived the study. F.C., J.V., D.I., C.S.L. and M.S.A. designed experiments. F.C., J.V., C.S.L. and M.S.-A. performed experiments and acquired data. J.V. and P.F. analyzed data. J.V., P.F. and D.I. developed and analyzed the model. D.I., F.C. and J.V. wrote the manuscript.
This work was supported by Ministerio de Economía y Competitividad [FU2012-34324 to F.C.]; and by a Swiss institute of Bioinformatics Fellowship to J.V.
The authors declare no competing or financial interests.