Space colonization by branching trachea explains the morphospace of a simple respiratory organ

Depósito de investigación de la Universidad de Sevilla https://idus.us.es/ “This is an Accepted Manuscript of an article published by Elsevier in Developmental Biology on 1 June 2020, available at: https://doi.org/10.1016/j.ydbio.2020.02.005 .” 1 Space colonization by branching trachea explains the morphospace of a simple respiratory organ. A. Ruiz-Sobrino 1*, C. A. Martín-Blanco 1*, T. Navarro 1*, I. Almudi 1, G. Masiero 1, M. Jimenez-Caballero 1, D. B. Buchwalter 3, D. H. Funk 4, J. L. Gattolliat 5,6, M. C. Lemos 2, F. Jiménez 2# & F. Casares 1# 1- CABD, GEM-DMC2 Unit (CSIC-Pablo de Olavide University-Junta de Andalucía), 41013 Seville, SPAIN. 2- Department of Condensed Matter Physics, University of Sevilla, 41012 Sevilla, SPAIN. 3- North Carolina State University, Department of Biological Sciences, Raleigh, NC 27695, USA. 4- Stroud Water Research Center, Avondale, PA 19311, USA. 5- Musée Cantonal de Zoologie, CH-1014 Lausanne, SWITZERLAND. 6- University of Lausanne (UNIL), Department of Ecology and Evolution, CH-1015 Lausanne, SWITZERLAND. (*) These authors contributed equally. (#) Senior authors. fcasfe[email protected]; [email protected] Keywords: Branching morphogenesis; mayfly gill; Cloeon; space colonization algorithm; computational models; branchless; fibroblast growth factor. ABSTRACT Branching morphogenesis helps increase the efficiency of gas and liquid transport in many animal organs. Studies in several model organisms have highlighted the molecular and cellular complexity behind branching morphogenesis. To understand this complexity, computational models have been developed with the goal of identifying the “major rules” that globally explain the branching patterns. These models also guide further experimental exploration of the biological processes that execute and maintain these rules. In this paper we introduce the tracheal gills of mayfly (Ephemeroptera) larvae as a model system to study the generation of branched respiratory patterns. First, we describe the gills of the mayfly Cloeon dipterum, and quantitatively characterize the geometry of its branching trachea. We next extend this characterization to those of related species to generate the morphospace of branching Manuscript File Click here to view linked References 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 2 patterns. Then, we show how an algorithm based on the “space colonization” concept (SCA) can generate this branching morphospace via growth towards a hypothetical attractor molecule (M). SCA differs from other branch-generating algorithms in that the geometry generated depends to a great extent on its perception of the “external” space available for branching, uses few rules and, importantly, can be easily translated into a realistic “biological patterning algorithm”. We identified a gene in the C. dipterum genome (Cd-bnl) that is orthologous to the fibroblast growth factor branchless (bnl), which stimulates growth and branching of embryonic trachea in Drosophila. In C. dipterum, this gene is expressed in the gill margins and areas of finer tracheolar branching from thicker trachea. Thus, Cd-bnl may perform the function of M in our model. Finally, we discuss this general mechanism in the context of other branching pattern-generating algorithms. SIGNIFICANCE STATEMENT We introduce the gills of mayflies as model to study the generation and variation of respiratory branching patterns and show how a space colonization algorithm can explain their morphogenesis. 1. INTRODUCTION Branching patterns are pervasive in nature and essential to the efficient distribution of liquids and gasses in biological systems. Within animals, branching tubular systems are found in organs as diverse as lungs, kidneys, blood vessels, mammary glands and insect trachea, where they increase the functional surface area for gas or liquid exchange. In general, the branching patterns found in these organs show complex, self-similar tridimensional geometries [1]. When studied at the molecular and cellular detail, the mechanisms generating the respiratory trees in mouse or in Drosophila seem exceedingly complex [2]. For example, in the Drosophila embryo Ghabrial and coworkers found in a genetic screen over 200 patterning/morphogenetic genes, beyond house-keeping genes, that were necessary for shaping and developing the arborization of the tracheal system [3]. To understand this complexity, models aimed at achieving a systems-level understanding have been developed. The goal of these models is to identify the minimal set of processes that are capable of generating and maintaining the branching pattern. This global understanding, rather than an understanding of the individual parts, may be then used to focus further research on the mechanisms of those essential processes, aiding in posing the relevant experimental questions. In addition, the branching morphogenesis of an organ may vary depending on environmental conditions, or between homologous organs in related species (e.g. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 3 mouse and human). Therefore, these models should also help explain in the most parsimonious manner how the phenotypic plasticity of a branching pattern and its evolutionary variation (i.e. the pattern’s morphospace) may arise. Models of animal branching morphogenesis fall into three major categories (reviewed in [1, 4]). The realization of the self-similar or fractal nature of these patterns [5, 6] led to the development of geometric models using recursive algorithms to quantitatively reproduce the general geometry of the lung’s arborization and its interspecific variation [7-9]. However, for these types of models it is often difficult to link the recursive rules of the patterning-generating algorithm to the biological causes. To address this limitation, a second class of models, of reaction-diffusion type, has been developed that link local geometry and gene regulatory networks. Therefore, in these models, molecular details are included that allow for the investigation of the initiation of the branching process, such as bronchial budding and elongation [10, 11] and the branching pattern of the kidney ducts [12]. A third class is that of models that investigate the role of mechanical forces as agents in the branching process and, in general, have a focus on the biomechanical properties of the tube structures as emerging from cell- and tissue-level properties (see, for example, [13]). It is clear that elements of all three modeling approaches are critical to developing a robust understanding of branching morphogenesis and should be combined when possible. For example, a recent geometry-based model, that also incorporates molecular information on relevant intercellular signaling processes, has been used to explain the global morphology of the branched ducts of the mammary gland [14]. Still it is not clear yet whether the mechanisms identified by these models are universal –e.g. whether they can explain the branching patterns of both ductal systems, which are specialized in liquid transport, and those of the respiratory arborizations of lungs or insect trachea that transport gases– or to what extent they are able to generate the developmental plasticity and the evolutionary diversity of these patterns. In this paper, we capitalize on the branching patterns of trachea in the mayfly gill to investigate the global mechanism that drive their generation and morphological evolution. 2. RESULTS 2. 1. The branching pattern and morphology of the tracheal gills of Cloeon and other mayflies. The aquatic nymphs of the Baetid mayfly Cloeon dipterum bear respiratory gills on abdominal segments I to VII. These gills are movable, flat epithelial bilayers (Figure 1) traversed by ramifying trachea [15, 16]. In their structure and function, they are similar to a small bi-dimensional lung, with gas being exchanged between the thin distal 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 4 trachea (secondary trachea) and the surrounding water through the epithelium as gas is exchanged between the alveoli and the blood through the capillary endothelium. Gills on segments I to VI are bi-lamellate, while those on segment VII are uni-lamellate (and the only non-movable gill pair) (Figure 1B). Cloeon gills are oval shaped and generally considered to be of the ancestral type [16, 17](Figure 1B). Gills start developing on the second nymphal stage and soon thereafter trachea invade the pad and start branching (Figure 1C,D). The onset of gill development coincides with some other developmental processes, such as the growth of the medial cercus and the wing pads. Globally, trachea diameter thins as the tubes get close to the periphery of the gill (distal portion). Also, and as development proceeds, new thin tracheae branch out from thicker, more internal (proximal) segments (Figure 1E). Therefore, tube thickness can be considered a correlate of branching time –with thinner tubes having branched more recently. Tracheal branching is asymmetric and trachea do not anastomose. As the gills grow, their pattern scales (Figure 1F,G). Although left and right gills from the same abdominal segment are similar to each other, they are not identical, indicating that the specific branching pattern is not totally deterministic (not shown). Globally, mayflies show a large variation in overall gill number, morphology and branching pattern [16]: the shape of the gills (fibrillose, bifurcate, spine-like, operculate...) is recognized as being of phylogenetic importance. Besides respiration, gills can have additional functions such as osmoregulation, substrate attachment forming a sucking disc (in Rhithrogena) or mechanical protection (opercula in Caenidae, Leptohyphidae, Coryphoridae and others). Gills can be secondarily very reduced or even completely absent (Murphyella). Figure S1 shows a sample of gills from several related mayfly species that show a variety of branching patterns. One expectation from a biologically relevant patterning system is that it should be capable of generating this natural morphospace of branching patterns by small modifications in its parameters. 2.2. The Space Colonization Algorithm (SCA) generates branching patterns similar to those of biological gills. Since we still lack even a basic molecular understanding of gill development in mayflies, the approximation used in this work to understand the generation of this branching patterns has been macroscopic: trying to identify a simple algorithm that generates patterns similar to those of natural gills and, at the same time, which could be translated into the smallest possible set of biological mechanisms. The superficial resemblance of this pattern to a branching tree or root system set us to search for tree-generating algorithms. Most algorithms generating branched patterns are based on “recursive branching”: a series of rules that specify branching frequency, angles and relative branch length. These are all “internal” rules. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 5 However, Runions and colleagues, aiming at computer-drawing realistic tree- and shrub-like structures, used an alternative concept: “space colonization”[18, 19]. The fundamental difference with other branching-generating algorithms is that the geometry generated depends to a great extent on its perception of the “external” space available for growth. In brief, the space colonization algorithm works as follows (Figure 2): The space to be colonized by the branching structure is filled with attractive points (M); A branch will grow at a speed (CREC) towards M if it is within a certain distance from M (MAX: radius of influence) (or towards the mean position of several Ms when the branch is within the reach of these Ms). The points are then removed as the branch reaches closer to these points (within a “killing” radius, KILL). The starting point of the branching structure is a single node in the space filled with M attracting points (“root”), and the structure finishes when all these attracting points have been eliminated. Therefore, “available space” is perceived as remaining attracting points. We have implemented a 2D SCA algorithm in which the shape of the space to be colonized can be set to an ellipse, as a first approximation to the gill shape, or to any other (polygonal) shape. In addition, branching proceeds within a growing space. Although in its original form M points are placed at specific locations to determine final tree shape, in our implementation Ms are distributed randomly at a user-controlled density and, in addition, the position of each M can be affected by a stochastic motion (random walk, at speed vel, making it a dynamic space colonization algorithm, or dSCA). (Figure 2A; see also Methods for the algorithm’s pseudo-code and Supplementary Material for detailed descriptions). Figure 2B-B’’ shows three time points in the development of a “synthetic gill” generated by this SCA implementation. Branching proceeds outwards, asymmetrically, and without crossings, recapitulating the pattern of a “natural” gill. In what follows, we develop a series of metrics to compare, first, between natural and synthetic gills, to then study how some major parameters affect their geometry. 2. 3. Geometric characterization of gill branching patterns. In order to characterize the branching geometry of gills, we proceeded in three steps. First, we developed a pipeline for image segmentation (see Materials and Methods and Supplementary Material) that resulted in skeletonized images where every branch has a pixel width and branching nodes and branch tips are located in the image (Figure S2A). Second, we developed a battery of metrics to quantify the differences we are able to detect visually between branching patterns across species (Materials and Methods). Each of these metrics defines some aspect of the pattern and was developed when a previous metric was unable to discriminate between two patterns. For example, we started by defining “branch density” (branching pixels relative to total pixels (i.e. gill 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 6 surface)) as a metric able to differentiate some patterns. However, it is easy to imagine two different patterns with identical branch density (Figure S2B). Then, a new one metric, “Branch tips” (tips to surface ratio), is now able to differentiate them (Figure S2B). Following this logic, we developed a series of 8 metrics that characterize the branching geometry (see Materials and Methods for full description), including its fractal dimension. However, in a given dataset it may be that some of these metrics are partly redundant (i.e. there is correlation between them). In this case, the comparison of two or more gill sets can be done using Principal Component Analysis (PCA), so that the variability between sets leads to the selection of a specific linear combination of metrics that is used to compare the patterns. The two principal components explain a large fraction of the variability (Figure S2C and Supplementary Materials) and were renamed as “richness” (x) and “irregularity” (y) (Figure 3). “Richness” reflects the density of different elements (nodes, tips and density mainly) while “irregularity” reflects the degree of heterogeneous branch lengths and unequal distribution of nodes of the pattern. Figure 3 represents the distribution of natural gills according to these two principal components. In this morphospace the richer sampling is of Cloeon dipterum, Baetis (Rhodobaetis) tenerifensis and Siphlonurus quebecensis, for which we analyzed all abdominal gills (A1 to A7) from one individual per species, including the small lamellae (A1-A6 in Cloeon and A1 and A2 in Siphlonurus) (Figure S3), but it also includes gills from other species (see phylogenetic tree in Figure S1). The metrics can vary one order of magnitude between different patterns. For example, the density of tips per gill (“branch.tips”) is one of the metrics that contributes the most to the “richness” component and ranges from 5X10-3 to 2X10- 4 between Alainites and Siphlonurus, this latter the species with the “richest” patterns in our sample. Within the second component, “irregularity”, “average number of ancestors” measures how the nodes distribute along a branch and ranges from 0,18 in Alainites to 1,86 in Neocleon. Globally, this quantitative analysis allows the discrimination of different branching patterns and describes the gill morphospace of our sample of species. 2.4. Exploration of SCA-generated branching patterns. In order to explore the range of branching patterns generated by the SCA algorithm, we simulated synthetic gills varying one of these parameters at a time. As availability of M is key, as it defines “unoccupied space”, we first studied the effects of varying M concentrations. As expected, the total branch length saturates as [M] increases. More interesting is the fact that, at intermediate concentration the pattern varies quickly, from sparse and undulated (that we have not observed in nature, as in all species examined trachea 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 7 grow straight) to patterns resembling natural gills (Figure 4A). Therefore, within a range of M, relatively small variations in M can result in qualitatively different patterns. Another parameter is gill morphology. However, when synthetic gills differing only in their ellipticity (i.e. with all dSCA parameters remaining equal) were generated, the resulting branching patterns are qualitatively similar (Figure 4B). Finally, we explored varying the two radii, MAX (detection of M) and KILL (the radius of killing of M). While MAX above 40 generate similar realistic patterns without much qualitative change (Figure 4C), decreasing KILL increases the richness of the arborization (Figure 4D). In order to statistically assess the discrimination power of our PCA analysis, we took the 3 species with the largest sample size (Cloeon, Baetis and Siphlonurus) and compared them through different methods (see Materials and Methods for details). First, we performed an unrestricted multinomial regression model to test whether the values of the two first PCA components of each gill pattern (alone or in combination) were able to predict the species this pattern belonged to. The model using the combination of Comp.1 and Comp.2 was a perfect classifier of species, with 100% of success in predictions (see Table S1). When each PCA component was tested alone, even the Comp.1 is enough for a perfect discrimination among these groups, while Comp.2 predicted the species correctly with a 72.5% of success. Therefore both components are sensitive enough for species-specific gill pattern discrimination. Finally, just working with Comp.1, an ANOVA analysis was performed to confirm Comp.1 means are significantly different between the 3 groups (p-value∼0) and significant too in pairwise comparisons (Adjusted p-value∼0). Therefore, the morphospace described by our PCA analysis is capable of discriminating the different species-specific patterns in a statistically significant manner. 2.5. The range of SCA-generated patterns covers the natural morphospace. To have a clearer idea of the impact each parameter had on the pattern, and to compare synthetic and natural gills, we quantified and represented their metrics within the twocomponents morphospace (Figure 5). We found that variations in the SCA parameters generated patterns that cover the whole range of branching patterns in the species we have surveyed, from simple and irregular (such as some found in Baetis tenerifensis) to more complex and ramified (like those of Siphlonurus quebecencis). An especially sensitive parameter is KILL, the distance at which a growing branch erases the attracting M: reducing the KILL radius results simultaneously in richer and more irregular patterns. Therefore, a space colonization algorithm reproduces quantitatively the geometry of the whole morphospace of natural gills we have analyzed. The full set 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 8 of images (natural gills and SCA generated patterns) and extracted datasets are provided as Supplementary Material. Software (segmentation macros, raw metrics generating algorithm and spatial colonization algorithm) is available through GitHub (see Supplementary Materials). 2.6 Cloeon FGF is expressed in uncolonized areas of the gill. The FGF ligand branchless (bnl) is known to stimulate the growth and branching of embryonic trachea and the thoracic air sac primordium in Drosophila [20, 21]. Therefore, bnl could be, in principle, a good candidate for M. However, in both instances bnl is produced at specific sites to direct tracheogenesis, while our model does not require M transcription to be actively patterned. Instead, the pattern arises from the repressing action of growing trachea. Therefore, our model predicts M to be transcribed along the margin of the gill (where most of the finer trachea are growing into) and in the inter-trachea regions with lowest branching density. To test this prediction, we identified the bnl orthologue in Cloeon dipterum’s genome [22, 23](Figure s4 and Table S2) with which we generated an anti-sense probe. In situ hybridization revealed that indeed Cd-bnl is expressed most strongly along the margin, where the tips of trachea are growing, and around the fine tracheoles branching from thicker tubes, in the gill internal regions (Figure 6). Indeed, this pattern is compatible with the pattern expected for M in the SCA algorithm (see Figure 2B’’), which should be present in regions not yet colonized (or in the process of being colonized) by trachea. 3. DISCUSSION In what follows we will try to translate our implementation of the SCA algorithm into plausible biological processes. Although our model is not molecularly specific, we will draw hypotheses on potential molecular mechanisms from the model and from prior research. Finally, we will compare the space colonization model with previous models explaining tubular branching patterns to highlight the most important differences. The space colonization algorithm is based on very simple rules: attracting molecules, M, direct the growth/branching of trachea which, when reaching closer, eliminate these attracting Ms. This algorithm results in a directional growth of the pattern and yet it does not require external sources of positional information, such as signaling centers or gradients: M would be distributed uniformly and the arborization will grow from its root -the entry point of the trachea into the gill. “[M]” denotes an effective concentration, not the actual concentration of the molecular attractant. For example, higher M might reflect an increase in the affinity for M in tracheal cells. We obtain 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 15 update the growth vector of this branch; } } for (all branches){ if (there are attracting M){ branch or grow: extend a new branch of length CREC along direction v; } } while(vel) for(all M) move M randomly; if(number of M == 0){ save data; stop and finish program; } } } The full set of images (natural gills and SCA generated patterns) and extracted datasets are provided as Supplementary Material. Software (segmentation macros, raw metrics generating algorithm and spatial colonization algorithm) is available through GitHub (https://github.com/Carlosmarbla/Spatial-Colonization-Algorithm, and see Supplementary Materials). 5. FUNDING Grants: 1455906 (National Science Foundation, USA) to DHF; TIN2017-89842P (MINECO/AEI, Spain, co-financed by FEDER funds (EU)) to FJ and MCL; 657732 (Marie Sklodowska-Curie Program, EU) to IA; BFU2015-66040-P and PGC2018- 093704-B-I00 (MINECO/AEI, Spain co-financed by FEDER funds (EU)) to FC and institutional grant MDM-2016-0687 (MINECO/AEI, Spain co-financed by FEDER funds (EU)) in which FC is participant researcher. 6. AUTHOR’S CONTRIBUTIONS FC conceived the study and obtained major funding; FC and FJ supervised the study; FJ and ARS wrote software for the SCA algorithm; CMB wrote macros for ImageJ (processing and segmentation) of gill images; ARS and TN developed metrics for geometry analysis; TN carried out statistical analysis; ARS, FC, CMB, IA, DHF, J-LG generated data; FC, ARS, FJ, TN, DHF, J-LG, MCL, DBB analyzed data and results; FC wrote the manuscript with the help of DBB and inputs from all authors. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 16 7. DECLARATIONS OF INTEREST. The authors declare no competing interests. 8. REFERENCES 1. Varner, V.D., and Nelson, C.M. (2017). Computational models of airway branching morphogenesis. Semin Cell Dev Biol 67, 170-176. 2. Metzger, R.J., Klein, O.D., Martin, G.R., and Krasnow, M.A. (2008). The branching programme of mouse lung development. Nature 453, 745-750. 3. Ghabrial, A.S., Levi, B.P., and Krasnow, M.A. (2011). A systematic screen for tube morphogenesis and branching genes in the Drosophila tracheal system. PLoS Genet 7, e1002087. 4. Iber, D., and Menshykau, D. (2013). The control of branching morphogenesis. Open biology 3, 130088. 5. Mandelbrot, B.B. (1983). The Fractal Geometry of Nature, (New York: W.H. Freeman). 6. Glenny, R.W. (2011). Emergence of matched airway and vascular trees from fractal rules. J Appl Physiol (1985) 110, 1119-1129. 7. Nelson, T.R., and Manchester, D.K. (1988). Modeling of lung morphogenesis using fractal geometries. IEEE Trans Med Imaging 7, 321-327. 8. Nelson, T.R., West, B.J., and Goldberger, A.L. (1990). The fractal lung: universal and species-related scaling patterns. Experientia 46, 251-254. 9. Kitaoka, H., and Suki, B. (1997). Branching design of the bronchial tree based on a diameter-flow relationship. J Appl Physiol (1985) 82, 968-976. 10. Kondo, S., and Miura, T. (2010). Reaction-diffusion model as a framework for understanding biological pattern formation. Science 329, 1616-1620. 11. Menshykau, D., Blanc, P., Unal, E., Sapin, V., and Iber, D. (2014). An interplay of geometry and signaling enables robust lung branching morphogenesis. Development 141, 4526-4536. 12. Menshykau, D., Michos, O., Lang, C., Conrad, L., McMahon, A.P., and Iber, D. (2019). Image-based modeling of kidney branching morphogenesis reveals GDNF-RET based Turing-type mechanism and pattern-modulating WNT11 feedback. Nature communications 10, 239. 13. Kim, H.Y., Varner, V.D., and Nelson, C.M. (2013). Apical constriction initiates new bud formation during monopodial branching of the embryonic chicken lung. Development 140, 3146-3155. 14. Hannezo, E., Scheele, C., Moad, M., Drogo, N., Heer, R., Sampogna, R.V., van Rheenen, J., and Simons, B.D. (2017). A Unifying Theory of Branching Morphogenesis. Cell 171, 242-255 e227. 15. Craig, D.A. (1990). Behavioral hydrodynamics of Cloeon dipterum larvae (Ephemeroptera: Baetidae). Journal of the North American Benthological Society 9, 346-357. 16. Riek, E.F. (1973). The classification of the Ephemeroptera. In First International Conference on Ephemeroptera. (Leiden, The Netherlands), pp. 160-178. 17. Kukalova-Peck, J. (1985). Ephemeroid wing venation based upon new gigantic Carboniferous mayflies and basic morphology, phylogeny, and metamorphosis of pterygote insects (Insecta, Ephemerida). CAN. J. ZOOL. 63, 933-955. 18. Runions, A., Brendan, L., and Prusinkiewicz, P. (2007). Modeling Trees with a Space Colonization Algorithm. In Eurographics Workshop on Natural Phenomena, D. Ebert and S. Mérillou, eds. (The Eurographics Association). 19. Runions, A., Fuhrer, M., Lane, B., Federl, P., Rolland-Lagan, A.G., and Prusinkiewicz, P. (2005). Modeling and visualization of leaf venation patterns. A. ACM Transactions on Graphics 24, 702–711. 20. Metzger, R.J., and Krasnow, M.A. (1999). Genetic control of branching morphogenesis. Science 284, 1635-1639. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 17 21. Sato, M., and Kornberg, T.B. (2002). FGF is an essential mitogen and chemoattractant for the air sacs of the drosophila tracheal system. Dev Cell 3, 195-207. 22. Almudi, I., Martin-Blanco, C.A., Garcia-Fernandez, I.M., Lopez-Catalina, A., Davie, K., Aerts, S., and Casares, F. (2019). Establishment of the mayfly Cloeon dipterum as a new model system to investigate insect evolution. EvoDevo 10, 1. 23. Almudi, I., Vizueta, J., de Mendoza, A., Wyatt, C., Marletaz, F., Firbas, P., Feuda, R., Masiero, G., Medina, P., Alcaina, A., et al. (2019). Genomic adaptations to aquatic and aerial life in mayflies and the origin of wings in insects. In BioRxiv. 24. Ribeiro, C., Ebner, A., and Affolter, M. (2002). In vivo imaging reveals different cellular functions for FGF and Dpp signaling in tracheal branching morphogenesis. Dev Cell 2, 677-683. 25. Okenve-Ramos, P., and Llimargas, M. (2014). Fascin links Btl/FGFR signalling to the actin cytoskeleton during Drosophila tracheal morphogenesis. Development 141, 929-939. 26. Volckaert, T., Campbell, A., Dill, E., Li, C., Minoo, P., and De Langhe, S. (2013). Localized Fgf10 expression is not required for lung branching morphogenesis but prevents differentiation of epithelial progenitors. Development 140, 3731- 3742. 27. Nogawa, H., and Ito, T. (1995). Branching morphogenesis of embryonic mouse lung epithelium in mesenchyme-free culture. Development 121, 1015-1022. 28. Spurlin, J.W., 3rd, and Nelson, C.M. (2017). Building branched tissue structures: from single cell guidance to coordinated construction. Philos Trans R Soc Lond B Biol Sci 372. 29. Ochoa-Espinosa, A., and Affolter, M. (2012). Branching morphogenesis: from cells to organs and back. Cold Spring Harbor perspectives in biology 4. 30. Davies, J.A. (2002). Do different branching epithelia use a conserved developmental mechanism? Bioessays 24, 937-948. 31. Sweeney, B.W., Funk, D.H., Camp, A.A., Buchwalter, D.B., and Jackson, J.K. (2018). Why adult mayflies of Cloeon dipterum (Ephemeroptera : Baetidae) become smaller as temperature warms. Freshwater Science 37, 64–81. 32. Rueden, C.T., Schindelin, J., Hiner, M.C., DeZonia, B.E., Walter, A.E., Arena, E.T., and Eliceiri, K.W. (2017). ImageJ2: ImageJ for the next generation of scientific image data. BMC Bioinformatics 18, 529. 33. Arganda-Carreras, I., Kaynig, V., Rueden, C., Eliceiri, K.W., Schindelin, J., Cardona, A., and Sebastian Seung, H. (2017). Trainable Weka Segmentation: a machine learning tool for microscopy pixel classification. Bioinformatics 33, 2424-2426. 34. Napolitano, A., Ungania, S., and Cannata, V. (2012). Fractal Dimension Estimation Methods for Biomedical Images. In MATLAB – A Fundamental Tool for Scientifi c Computing and Engineering Applications, Volume 3. pp. 161-178. 35. Venables, W.N., and Ripley, B.D. (2002). Modern Applied Statistics with S. Fourth Edition. , (New York: Springer). 36. de Mendiburu, F. (2019). agricolae: Statistical Procedures for Agricultural Research. R package version 1.3-1. . 9. FIGURE LEGENDS Figure 1. The tracheal gills of the mayfly Cloeon dipterum and the development of their arborization. (A) The aquatic nymph bears gills on its first seven abdominal segments. The white ellipse marks one of them. Scale bar is 1mm. (B) Microscopic images of abdominal gills 1, 5 and 7. Except for the pair on abdominal 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 18 segment 7, gills are bi-lamellate. The scale bar in (B) is 500 μm. (C) DAPI stained gill IV of a 5 day old nymph. No signs of trachea are visible yet. (D) Bright field micrograph of a 13 day old nymph gill III where the trachea (arrows) have invaded the epithelial parenchyma of the posterior lamella (p) through its base. The anterior lamella (a) is already visible stemming from the base of the posterior one. (E) In later stages, small tracheoles (arrowheads) branch out from thicker, older trachea (asterisks). (F-G) Gills VII (from different individuals) at 34 (F) and 38 (G) days post-hatching. Scale bars are 100μm. The branching pattern scales with growth (note that the length/width of the gill approximately doubles from F to G and includes further elaborations, mostly distally). Figure 2. An algorithm based on Space Colonization generates realistic gill branching patterns. The space colonization algorithm (SCA; A-A’’): (A) A segment (blue) that falls within the radius of influence of an attracting point “M” grows in this direction. When the segment is close enough to the point, the point is eliminated. The process continues (A’) as long as there are more attracting points within reach. “MAX” is the radius of influence of M (i.e. the distance at which the segment can detect M), “KILL” is the radius of “killing” (or erasure) of M and “crec” is the increase in length of the branch per iteration. (A’’) The SCA algorithm, with a seed located in the center of a space filled with M points (red), generates a branching pattern (blue). (B-B’’) Three time-points along the development of a SCA-generated “synthetic” gill. M points in red. The younger the branches, the lighter the blue color. (C) A “natural” Cloeon VII gill is shown for comparison. Figure 3. Natural morphospace of the gill branching patterns. Distribution of the species according to the two principal component analysis. A sample skeletonized image for each species, color-coded, is included. Gills of different species fall in different regions of this morphospace. Figure 4. Effects of M and global shape on gill branching pattern. (A) Effect of increasing M concentration (M/pixels) on total trachea length. Parameters (40,11,2,1) (MAX, KILL, CREC, vel), constant shape (semiaxis (a,b)= (200,230). Red dots, simulations. Some graphic outputs (“synthetic” gills) are shown. (B) Effects of varying ellipticity (short semiaxis/long semiaxis of B1=128/230 and of B2=228/230 ) using (MAX, KILL, CREC, vel)= (40, 11, 2, 1), as in (A). (C) Variation of MAX. (D) Variation of KILL. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 19 Figure 5. Variations in parameters and their effects on branching geometry. Plot of the PCA analysis including natural gill patterns (grey circles, data in Figure 3) and simulated patterns varying the main SCA parameters. Skeletonized images shown are variations in KILL (top row, from left to right: 40,20,2,1; 40,12,2,1; 40,5,2,1; 40,2,2,1; blue symbols) and variations in CREC (left, column, top to bottom: 40,5,1,1; 40,5,2,1; 40,5,4,1; 40,5,8,1; orange symbols). Other simulations include variations of shape (ellipticity, exc, green), MAX (30,5,2,1; 40,5,2,1; 80,5,2,1; 100,5,2,1; 200,5,2,1; 250,5,2,1; pink symbols) and M (as in Figure 4). The circle symbols are Baetis (light green) and Siphlonurus (dark green) gills. Figure 6. Expression pattern of Cloeon bnl in gills. (A) White field image of a gill hybridized with an anti-sense bnl/FGF probe. The signal is stronger along the margin, where the distal tips of the growing trachea are located. In addition, signal is detected at, or around, the finer tracheoles in the inner region of the gill. (A') DAPI staining of the same gill to visualize cell nuclei. (*) indicates lateral gill regions still covered by cuticle. The cuticle precludes staining. “P” and “D” mark “proximal” and “distal”, respectively. (See Figure S4 and Supplementary Materials and Methods for further details on Cloeon dipterum bnl sequence identification). 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Figure Click here to download Figure Figure1.tiff Figure Click here to download Figure Figure2.tiff Figure 3R Click here to download Figure Figure3.tiff Figure Click here to download Figure Figure4.tiff Figure Click here to download Figure Figure5.tiff