Biological tubes are essential for animal survival, and their functions are dependent on tube shape. Analyzing the contributions of cell shape and organization to the morphogenesis of small tubes has been hampered by the limitations of existing programs in quantifying cell geometry on highly curved tubular surfaces and calculating tube-specific parameters. We therefore developed QuBiT (Quantitative Tool for Biological Tubes) and used it to analyze morphogenesis of the embryonic Drosophila trachea (airway). In the main tube, we find previously unknown anterior-to-posterior (A-P) gradients of cell apical orientation and aspect ratio, and periodicity in the organization of apical cell surfaces. Inferred cell intercalation during development dampens an A-P gradient of the number of cells per cross-section of the tube, but does not change the patterns of cell connectivity. Computationally ‘unrolling’ the apical surface of wild-type trachea and the hindgut reveals previously unrecognized spatial patterns of the apical marker Uninflatable and a non-redundant role for the Na+/K+ ATPase in apical marker organization. These unexpected findings demonstrate the importance of a computational tool for analyzing small diameter biological tubes.
Epithelial and endothelial tubes are essential to the survival of almost all multicellular organisms, and the physical properties of a tube, including size and shape, are crucial to its function. Mis-regulation of tube development or maintenance can have severe implications, including polycystic kidney disease in humans (Harris and Torres, 2009; Song et al., 2017). However, the mechanisms by which small tubes control their lengths and diameters is largely unknown (reviewed by Affolter and Caussinus, 2008; Iruela-Arispe and Beitel, 2013; Wang et al., 2017). Even basic descriptions of how cell shape and organization change during development, which is a fundamental prerequisite for understanding size control mechanisms, has not been determined for most tubes.
A major contributor to the lack of characterization of cell shape and organization in tubular epithelia is that existing computational tools for quantifying epithelial cell shape and organization (Barbier de Reuille et al., 2015; Förster and Luschnig, 2012; Khan et al., 2014) are not optimized for tubular epithelia, particularly for small diameter tubes with high curvature. Moreover, currently available programs do not calculate essential tube parameters such a centerline, branch points, cross-sectional areas and cell orientation, and do not unroll or otherwise project tube lumenal surfaces. To address this problem, we have created QuBiT (Quantitation Tool for Biological Tubes), a set of command line-based computational tools for measuring cell and tube morphology.
As a testbed for developing QuBiT, and to identify novel biological processes, we used QuBiT to characterize the early development of the Drosophila trachea system, which is one of the best-studied systems of tubular epithelia (reviewed by Manning and Krasnow, 1993; Samakovlis et al., 1996). The tracheal system is the gas exchange organ of the fly and thus functions as a lung, but its branch structure more resembles a vascular system because it is a ramifying network that directly delivers oxygen to specific tissues. Tracheal tubes are epithelial monolayers that are approximately the size of small capillaries or kidney tubules in mammals, but there are no associated muscle cells, pericytes or other accessory cells that are known to contribute to tracheal tube size control. Thus, tracheal tube size directly results from interactions of the tracheal cells with each other and with a secreted apical extracellular matrix (aECM) that transiently fills the tube lumens as they expand during their initial development.
Using QuBiT, we obtained several unexpected results, including: (1) anterior-to-posterior (A-P) gradients are present in many cell characteristics, including apical orientation and aspect ratio; (2) there exists a periodicity at the tube segment level to these characteristics within the A-P gradient; (3) inferred cell intercalation during development dampens an A-P gradient of the number of cells per cross-section of the tube, but does not change the connectivity distributions of tracheal cells; (4) cell connectivity distributions in the main tracheal tube are not influenced by the complex shapes of, or possible stresses on, cells that interface the side branches with the dorsal trunk (DT); (5) the apical marker Uninflatable (Uif) has supracellular A-P stripes of higher expression in the trachea and hindgut; (6) the long isoform of Na+/K+ ATPase α subunit, ATPα, has a non-redundant role in levels and subcellular localization of the apical marker Uif. These results demonstrate both the utility of QuBiT for analyzing tubular epithelia and the importance of quantitative analysis in understanding the cell biology of tubular epithelia.
Overview of analysis using QuBiT
To maximize maintainability, accessibility and extensibility of a tool for epithelial tube analysis, we developed QuBiT using commonly available and well-supported software platforms, rather than develop entirely new programs. At present, QuBiT uses a predominantly command line interface. It uses the signals from markers of the tube lumenal surface and cell junctions to define the lumenal surface and demarcate individual cell surfaces. QuBiT then calculates tube- and cell-specific parameters. Although this approach does not yield a full 3D reconstruction of the entire cell bodies that comprise a tube, it focuses on the apico-lateral junctions and apical regions that control tube size in tubes such as the Drosophila tracheal system (Beitel and Krasnow, 2000; Laprise et al., 2010; Sollier et al., 2015; Wodarz et al., 1995) and greatly simplifies the reconstruction problem. Moreover, as many tubes, including endothelial blood vessels and larval tracheal tubes, have thin cell bodies, the apical surface can closely approximate the location and shape of the entire cell.
Fig. 1A shows a schematic of the workflow. Image stacks are generated by confocal microscopy using settings that produce cuboidal voxels (Fig. 1Bi). Image segmentation is performed on the entire stack using Ilastik, a general-purpose image segmentation program (Kreshuk et al., 2011). We then analyze segmented images using custom-written code in Matlab (MathWorks) (open source available at http://github.com/gjbeitel/QuBiT). Tube analysis proceeds by segmenting the boundary of the tube lumen and creating a skeleton, which enables robust calculations of parameters of interest, including length, surface area and cross-sectional area (Fig. 1Bii, gray tube). Separately, cell junctions are masked onto the tube surface, resulting in apical cell surfaces that can directly be analyzed for parameters such as size and orientation (Fig. 1Biii).
QuBiT – data collection
QuBiT is designed to analyze images collected with cuboidal voxels containing either tube surface or junctional information with enough resolution for the desired analysis. For the measurements of tracheal length described below, basic length measurements were obtained using image stacks of entire whole-mount embryos collected using a 40× oil objective with a 0.38 µm voxel size (Fig. 1Bi). Between 75 and 250 optical sections were collected per embryo. For cell parameter analyses, images were collected using a 100× oil objective with 0.15 µm voxels. Two image stacks were tiled to recreate an entire tracheal tube.
QuBiT – surface mapping and defining tube centerlines
QuBiT defines a tube surface using markers that either visualize the cell lumenal/apical surface or the lumenal contents themselves (Fig. 1Bi,Bii). For the Drosophila trachea, we used either Uif (Zhang and Ward, 2009) or the aECM marker Vermiform (Verm) (Wang et al., 2006). An image stack is then segmented using the carving module in Ilastik, which allows isolation and extraction of a tube as a 3D binary image file in HDF5 format. This file is then directly imported to Matlab, and we applied a marching cubes algorithm (Cline et al., 1987; Hassouna and Farag, 2007) to ‘skeletonize’ or create the centerlines of the primary DT (Fig. 1Bii). If side branches are present, they are automatically included in the skeleton, and each branch intersection is used to identify a specific node or segment of the tube. The tube centerlines, which are not calculated by other epithelia analysis software, are of crucial importance because centerlines define the lengths of tube segments, mark branch intersection points that define tube segments, are the local reference for determining cell orientation, and are necessary for calculating orthogonal planes that are used to quantify tube parameters including diameter, cross-sectional area and surface area (e.g. Fig. 2C,D). To allow rapid investigation of different parameters without recalculation, computed results from the segmented data are stored with the image stack. For determining tracheal tube parameters, we used ten cross-sections per segment, which samples the tube approximately once every 15 voxels or 2.3 µm.
QuBiT – cell mapping
To map the apical cell faces onto the tube surface, QuBiT uses information from apico-lateral junctions (Fig. 1Biii). For the tracheal system, we used the claudin Kune-Kune (Kune) as a marker of septate junctions (SJ), as we have previously done (Nelson et al., 2010), because Kune staining gives a more continuous signal along the lengths of junctions than does E-cadherin staining, which improves segmentation. Kune signal also overlaps sufficiently with the Uif signal to allow demarcation of the apical cell surfaces.
Junctional segmentation was performed using the pixel classification tool in Ilastik. We imported this data into Matlab, and masked and inverted the junctions on the tube surface, as defined by either the apical marker staining or the edge of a lumenal stain, to obtain apical cell surfaces (Fig. 1Biv). We then calculated cell parameters that included area, cell orientation and aspect ratio.
Because the small cell size, narrow diameter and junctional organization of the tracheal DT leads to unresolvable segmentation errors, we used several parameters to filter our tracheal cell data. We applied a radial filter to exclude cells on dorsal branches and transverse connectives, an axial position filter to remove cells outside of regions of interest on the DT, and a size filter to exclude improbably small or large cells (Fig. S1A). After these steps, we detected on average 19 to 23 cells per segment. Previous work has shown that there are ∼21-28 cells per tracheal segment, depending on A-P position and genotype (Beitel and Krasnow, 2000; Samakovlis et al., 1996), resulting in our segmentation efficiency of ∼75% in the trachea.
QuBiT – re-imaging: unrolling the tube to make a 2D surface model
To provide a planar representation of the cells that allows simultaneous visualization of the entire tube surface and direct application of a broad array of existing 2D analysis tools, QuBiT can computationally unroll and flatten tubes. To do this, we calculate the exact cross-section using the orthogonal plane to the local centerline at regular voxel-length intervals along the DT. We then extract the cell junction data from projections radiating from the centerpoint in each orthogonal plane and write the data on a 2D plane (Fig. 1Bv, detailed in Fig. 4A).
As an example of the power of simultaneous visualization of the entire tube surface by unrolling, a previously unreported A-P striped organization of Uif is strikingly visible in the unrolled tracheal tube but is not apparent in conventional confocal slices or projections of the trachea (Fig. 1Biv, detailed in Fig. 5).
As an example of an operation that is much easier to perform in 2D, we used watershed segmentation to reduce cell junction widths to single pixels (Fig. 1Bv). This creates a more accurate representation of cell interfaces and enables analysis of parameters such as cell connectivity (the number of cells any given cell touches) and the number of cells in a tube cross-section, which we used in the analysis of the tracheal system below. In future work, 2D projections could be used to track tracheal cells over time more easily than in 3D, which is important in determining changes in cell arrangements.
QuBiT – re-imaging: recalculating and rerolling a 2D projection into a 3D model
The process used in QuBiT to unroll a tube can also be used to recreate a 3D tube from a 2D projection (Fig. 1Biv). As an example of how this can be employed, we used tube unrolling and rerolling to test for evidence of active cell rearrangements during tracheal expansion. We computationally expanded embryonic stage 14 tracheal tubes to the size of stage 16 tubes, and then compared the 3D parameters of the computed and actual stage 16 apical cell surfaces.
Tubes were computationally expanded by projecting the raw unrolled 2D stage 14 cell data back onto 3D tube surfaces with corresponding lengths and radii of stage 16 tubes. To do this, we increased the number of pixels along the long (x) axis of each 2D projection proportionally to our measured changes in tube length using nearest neighbor interpolation. Then, at each position along the model tube axis, we sampled the 2D image for cell surfaces and wrote the data on the orthogonal (yz) circumference with the corresponding radius based on stage 16 tube radii, which was normalized with respect to axial position. This resulted in a re-creation of stage 14 cell surfaces after tube inflation to stage 16 with no other changes to cell apical geometry or topology. One limitation of the resulting model is that the projected tube is perfectly straight, whereas the original tubes had small irregularities, which introduces a small amount of error. However, because of the asymmetric expansion along the length of the tube, the alternative approach of expanding along a curved path would result in over- and/or underexpansion of some regions, which would also result in slight deviations. To determine whether systematic errors were introduced, we compared the effect of unrolling and rerolling stage 14 and 16 tubes (Fig. S2A,B). For both stages, there was excellent agreement between the original and rerolled tubes, with most parameters differing by no more than 6%, which is sufficient for the analyses described below.
QuBiT – data analysis
For statistical analysis of tube parameters, QuBiT takes advantage of the extensive tools of the Matlab framework to allow users to perform statistical analyses on individual tubes, and on multiple independent tubes simultaneously, without the need to export the data into other analysis packages. Notably, for tubular systems with reproducible features, QuBiT has the ability to track, align and compare these features. In the case of the trachea, we aligned individual DT segments rather than just normalizing to length along the tube. As detailed below, this enables calculations that reveal tracheal cell parameters differ not only between segments, but also within segments.
Tracheal DT dimensions change in A-P gradients, but posterior is not always bigger
We tested the utility of QuBiT by investigating the growth of wild-type (WT) tracheal DT that occurs during ‘tube inflation’, a ∼2.5 h event in which a burst of lumenal secretion (Tsarouhas et al., 2007) doubles tube diameter and increases tube length by ∼15% between stages 14 and 16 (Fig. 2A). We used QuBiT to analyze six stage 14 and six stage 16 DT, yielding a total of 692 and 852 identified cells, respectively, in segments 2-9. By comparison, Nelson et al. (2012) and Förster and Luschnig (2012) each analyzed ∼55 cells from segments 6 or 7 in stage 16 embryos. The order-of-magnitude increase in number of cells analyzed, combined with analysis of most DT segments at both early and late stages in tracheal development, allowed us to identify statistically significant patterns to cell apical orientation along the length of the trachea and developmental changes that were not possible with the limited data obtained by previous groups.
Our measurements using QuBiT show that the DT length increased by 13.7±1.3% (Fig. 2B). These data were not statistically different from our manual length measurements of the same embryos [250±18 µm by hand in Volocity (PerkinElmer) versus 263±6.9 µm using QuBiT at stage 16, P=0.17, two-sample t-test], but had the advantage of automatically generating measurements on each individual segment that allowed a more detailed analysis of the size growth of tracheal tubes than had previously been performed. For example, whereas one might have expected the individual DT segments to all expand equally, most DT growth resulted from increases in the lengths of segments DT4 through DT9 (+22±8%), and the remaining segments lengthened only slightly (Fig. 2B). This variability in segment length ratio is in contrast to the steady tube diameter ratio (Fig. 2D), which is consistent with previous observations that tube length and diameter are regulated independently (Beitel and Krasnow, 2000; Forster et al., 2010; Luschnig et al., 2006; Wang et al., 2006), and further indicates that length is individually regulated by segment. Interestingly, anterior segments are generally longer than posterior segments, despite anterior segments having fewer cells (e.g. ∼21 in DT3 versus ∼28 in DT9; Robbins et al., 2014). These results predict that anterior and posterior cells have different arrangements and/or apical shapes in their respective segments, which, as described below, is indeed the case.
Many DT cell properties have an A-P gradient
To investigate how the shapes of the tracheal cells might contribute to control of tube size, we measured shape parameters of tracheal cell apical surfaces at stages 14 and 16. We focused on the apical surface because the available evidence indicates that tracheal size is regulated at the apical cell surfaces (Laprise et al., 2010; Nelson et al., 2012; Olivares-Castiñeira and Llimargas, 2017; Robbins et al., 2014) and that the shape of the basolateral surface is not regulated (Beitel and Krasnow, 2000).
We first examined how the apical area of individual tracheal cells changes as a function of segmental location and stage. Before tracheal tube expansion at stage 14, tracheal cells in all segments have similar apical areas (Fig. 2E, blue line), and though there is some variability within each segment, there is no apparent A-P gradient. By the end of stage 16, the area of all tracheal cells increases (P<10−47, comparing stage 14 with stage 16, paired t-test), but posterior cell area increases more, resulting in a shallow A-P cell apical area gradient (Fig. 2E). This shallow gradient is in marked contrast to the much more pronounced and uniform A-P gradient in anterior versus posterior segment surface area where the A-P difference is almost 2-fold (Fig. 2C). The discrepancy between the fairly uniform apical cell area and graded segment surface area predominantly results from posterior segments being both shorter (Fig. 2B) and having more cells than anterior cells (∼21 cells in DT3 versus ∼28 cells in DT9) (Samakovlis et al., 1996).
We next investigated apical cell shapes by calculating their aspect ratios (Fig. 2F). In contrast to the mostly uniform apical cell area of tracheal cells, the apical shape of tracheal cells is not uniform. Anterior cells have higher aspect ratios than posterior cells (P<10−3 at both stages 14 and 16, unpaired t-test of DT2 versus DT10) but, overall, aspect ratio does not change significantly during development (P=0.21, paired t-test).
As with cell apical size and aspect ratio, apical cell orientation also shows an unexpected A-P gradient (Fig. 2G). It has previously been shown from an analysis of a small number of cells that the long axes of tracheal cell apical surfaces in tracheal segment 8 (DT8) are not aligned with the long axis of the tracheal tube: in agreement with these previous measurements, we found that the average cell orientation in DT8 was 47±28° (Nelson et al., 2012: 42±39°; Förster and Luschnig, 2012: 62±18°). However, when we measured cell apical orientation in all tracheal segments, we found that anterior cell apical surfaces tend to be oriented along the DT axis with an average angle of 38±26°, whereas posterior cell apical surfaces tend to be oriented along the circumference of the tube with an average angle of 63±23° (P<10−11, unpaired t-test using DT2 against DT9 at stage 16). The gradient also exists at stage 14 (P=0.005) and cell orientation differs significantly between stages 14 and 16 (P<10−14, paired t-test).
Many DT cell apical properties have segmentally repeating variations
Despite tube surface area increasing fairly linearly from anterior to posterior (Fig. 2D), apical cell areas and aspect ratios at stage 16 show a strongly periodic pattern, with cell area and aspect ratio having local minima close to the points at which the transverse connectives (TC) join the DT, and local maxima in the middle of a segment (Fig. 2E-F). Cell apical orientation also shows a sinusoidal pattern, but the phase is opposite to those of surface area and aspect ratio, with orientation having local maxima where the TC branches connect to the DT (Fig. 2G). This sinusoidal pattern is robust with respect to the window size used to calculate average cell apical surface area along the tube, with the periodicity being clearly visible with window sizes ranging from at least 0.2 to 0.5 of a segment (Fig. S1B-C). Notably, the sinusoidal pattern corresponds to the physical organization of the trachea. The DT is comprised of nine distinct segments that derive from nine clusters of epidermal cells during development. Rather than the tube being a continuous cobblestone of cells, the cobblestone pattern is punctuated by pairs of thin, washer-like cells called ‘fusion cells’ that join adjacent segments immediately posterior to where the TCs branch from the DT (Fig. 2H, arrowheads). DT cells close to the fusion cells have reduced apical sizes relative to cells more distant from fusion cells.
One explanation for this segmental sinusoidal periodicity could be that the fusion cells constrain tube growth such that the tube is wider between fusion points. However, WT stage 16 DT tubes do not show constrictions at the fusion cells (Fig. 2A,H), and the plotted measurements of tube cross-section area and radius do not show segmental periodicity (Fig. 2D). Thus, the differential areas of cells that are closer to or further away from fusion cells appear to arise from differences between the cells and/or difference in the local lumenal environment. As no periodicity has been observed in lumenal proteins or organization, it appears likely that the DT cells closer to fusion cells expand their apical surface area less in response to secreted lumenal contents than do DT cells in the middle of the segment. Consistent with this hypothesis, uniform overexpression of the glycoprotein Tenectin in the DT using the btl-Gal4 driver results in tracheal segments in which the diameter of the tube increases more between fusion cells than at fusion cells, with the change in diameter along the segment being a smooth curve rather than a step at the fusion cell (figure 6F in Syed et al., 2012).
A simple expansion model explains changes in cell parameters during tracheal inflation
Given that expansion of the tracheal length and diameter is highly asymmetric, with length and diameter increasing by 13% and 100% respectively, we asked whether the changes in cell size, aspect ratio and orientation could result simply from an inflation of the tube, as if one were to inflate a cylindrical balloon with cell apical surfaces drawn on the surface. This is an appropriate model for the tracheal system because there is no change in cell number during tracheal development (Samakovlis et al., 1996); therefore, changes in tube size must result from changes in cell shape and/or cell arrangement. As described above, we used QuBiT to computationally expand stage 14 tubes, and consequently the cells that make up the tube, to their respective stage 16 sizes (Fig. 3Aiii). We then analyzed the cell shape and orientation parameters of the resulting ‘computationally expanded’ cells and compared them with cells from stage 14 and with 16 tubes that had been similarly unrolled and rerolled, but without computational expansion (Fig. 3Ai and Aii, respectively, termed ‘model’).
The cell parameters for stage 16 computationally expanded tubes are shown as purple lines in Fig. 3B-D, with stage 14 and stage 16 models shown as blue and orange lines, respectively. As an internal control, we first examined cell apical area (Fig. 3B). Apical area in the computationally expanded tube and in the actual stage 16 data was comparable (P=0.15, paired t-test), with the computational model still showing the strong segmental periodicity that was observed in the actual data. Interestingly, cell aspect ratio (Fig. 3C) and orientation (Fig. 3E) were also comparable between in vivo stage 16 data and computationally expanded tubes (P=0.09 and 0.31 respectively, paired t-test), indicating that changes in cell apical shape during tracheal growth are consistent with tracheal cells passively expanding their apical surfaces in response to inflation of the lumen.
The distribution of DT cell connectivity does not change during tube expansion
If tracheal cells simply expanded their apical surface areas in response to lumen inflation, in the simplest case, there would be no changes in the cell organization and the connections between cells would remain static. However, previous work has demonstrated that cell intercalation occurs during tracheal DT growth (Förster and Luschnig, 2012). To test whether QuBiT can find evidence of such intercalations, and to investigate the possible effects of such intercalation, we used QuBiT to measure two quantifiable parameters of cell organization at an early and a late stage of tracheal tube expansion: the number of connections each tracheal cell has with neighboring cells (cell connectivity) and the number of cells per tube cross-section.
Cell connectivity can refer to both the number of cells that a given cell shares contacts with, as well as the specific connections a cell has with neighboring cells. Because we analyze data from discrete time points rather than time-lapse images, we focus on the number of contacts that tracheal cells have, and how that parameter changes along the length of the DT and during development. Changes in the distribution in the number of contacts cells have is indicative of migration or other rearrangement events. However, as detailed below, an unchanging connectivity distribution does not demonstrate an absence of cell rearrangement.
We used QuBiT's capabilities to unroll the tube to create 2D representations of tracheal cells on the apical surface. A simplified diagram depicting the unrolling algorithm is shown in Fig. 4A. Starting with the segmented apical surface of the DT (blue), we systematically dilated or eroded this image one pixel at a time, resulting in ‘shells’ or ‘masks’ (dotted lines) that were a fixed relative depth (R) with respect to the baseline surface. At the same time, at one-pixel intervals along the length of the DT (x), we calculated each orthogonal, or cross-sectional, plane. Within each plane, we calculated vectors that correspond to a radial sweep at interval angles (θ, purple lines) along the ‘circumference’. We sampled the circumference every 2° or 4° depending on the expected radius of the tube. For most tracheae, 4° intervals are sufficient; lower intervals do not meaningfully change the overall result. Masking the radial vectors at a particular depth allowed us to reconstruct the contents of the tube at that depth as a 2D image; doing this at multiple depths resulted in ‘image stacks’ of an unrolled tube, in which the mask depth R takes the place of the usual z-axis; circumferential position θ takes place of the usual y-axis; and the A-P position x occupies the x-axis (right half of Fig. 4A).
Using this method, we unrolled the SJ channel at a fixed depth (R=+2) for all tracheae, and then applied a watershedding algorithm to allow quantification of cell-cell connectivity (Fig. 4B, degrees of connectivity coded by color). One limitation of watershedding stage 16 trachea is that, because the lateral faces of the washer-like fusion cells in many places are so closely apposed that the junctional signal for the opposing sides merge with each other and the neighboring cells, those regions of the fusion cells are lost from the watershedded diagram. Thus, although unrolled junction images clearly reveal the washer-like fusion cells that join tracheal segments (e.g. Figs 1Biv, 4A right, 5C,G), fusion cells are not accurately represented in connectivity analyses. However, as fusion cells account for only two out of ∼24 cells per tracheal segment, and are consistently misrepresented in all analyses, this limitation does not significantly affect conclusions drawn from connectivity analyses.
Stage 14 and 16 tracheae are predominantly composed of a mix of pentagons and hexagons (Fig. 4C, teal and green bars). There was no obvious clustering of any particular type of polygon at either stage, and the distribution of polygon types is relatively uniform along the length of the tube.
The distribution of DT cell connectivity is not impacted by side branches
We also took advantage of QuBiT's ability to accurately map tube branch points to ask whether there were any particular patterns of apical connectivity associated with the points at which the side branches of the DB or TC exit the DT (black and white circles in Fig. 4B, respectively). No obvious patterns of cell apical surface connectivity or organization of surrounding cell connectivity were observed in 2D maps of tracheal tubes, indicating that neither the complex shapes of the individual cells that contribute to both the DT and the side branches, nor any forces resulting from tension on the side branches, significantly impact the connectivity surrounding DT cells.
Importantly, there was no significant change in the distribution of polygons between stages 14 and 16 (Fig. 4C, P=0.23, χ2 test). Thus, it appears that there is little or no change in cell connectivity during tube expansion. This result is consistent with intercalation events occurring infrequently enough that they cause little change to the overall connectivity, but is also consistent with the occurrence of intercalation events conserving their cell connectivity, such as occurs in the T1-T2-T3 transitions and rosette formation in epithelial sheets that drive cell intercalation during convergent extension in zebrafish (Sepich et al., 2000; Tada and Heisenberg, 2012), flies (Blankenship et al., 2006; Irvine and Wieschaus, 1994; Munjal et al., 2015; Par et al., 2014; Simões Sde et al., 2014; Zallen and Wieschaus, 2004) and vertebrates (Lienkamp et al., 2012). We further investigated these two possibilities below.
The number of cells per cross-section changes during tube expansion
To distinguish between a paucity of cell rearrangements versus active changes in cell organization that conserve cell connectivity, we determined the number of cells per cross-section of the tube. Infrequent cell intercalations should not significantly change the number of cells per cross-section but, in a tube, even intercalation events that conserve cell connectivity would lead to changes in the number of cells per tube cross-section if they occurred with significant frequency. We used our 2D diagrams to examine the number of cells in orthogonal cross-sections of the DT at stages 14 and 16 (Fig. 4D-E), which yields two notable observations.
First, similar to the patterns observed for cell aspect ratio, orientation and area, a sinusoidal pattern of the number of DT cells per cross-section emerges (Fig. 4F). Second, in contrast to the lack of change in connectivity distributions during development, the number of cells per cross-section in the DT changes from stage 14 to 16, most notably in the anterior and posterior (P<10−3 from TC1 to TC4 and P<10−10 from TC7 to TC10, paired t-test). Crucially, as our simple tube expansion model preserves the organization of tracheal cells, the number of cells per cross-section in the computationally expanded stage 16 tube is the same as the original stage 14 tube and is also statistically different from measured stage 16 tubes. Thus, although a simple apical expansion model was sufficient to account for the changes in cell parameters during DT expansion, our analytical software indicates cell intercalation events are actually common during tracheal maturation. QuBiT's ability to predict the existence of cell intercalation events demonstrates that QuBiT can be used to reveal developmental changes in cell organization, even under circumstances in which live imaging of tube morphogenesis is not feasible.
Tube unrolling reveals unexpected patterns of Uif localization
In addition to using tube unrolling to examine cell connectivity, we used unrolling to simultaneously visualize the entire DT tube surface and visually search for patterns in the organization of the apical marker Uif in WT stage 16 trachea (Fig. 5A). Surprisingly, whereas Uif was thought to be a uniform marker of the tracheal apical surface, in the unrolled projections it is immediately apparent that the tube surface has previously unrecognized A-P ‘bands’ of Uif signal, twice per circumference, that are elevated compared with the rest of the tube and that are not interrupted by the washer-like fusion cells that disrupt cuticle pattern at this stage (Fig. 5B, yellow arrowheads). Notably, these bands are in register with the TCs and DB branch exit points, with one band being aligned with the DB exits and the other with TC exits (Fig. 5B, white and black circles indicate TC and DB exits, respectively).
To verify these observations, we examined unrolled images at various depths (ranging from R=−6 to +10, the equivalent of −1 µm to +1.5 µm) from the segmented apical surface and found the same phenomenon at any surface level (data not shown). The banding does not result from a technical artifact such as imaging depth because no such banding is visible in the junctional staining channel (Fig. 5C), but banding is visible in tube cross-sections generated by the Volocity program (Fig. 5B, yellow arrowheads). Moreover, the dimmest regions do not correspond with the deepest sections. Closer examination reveals that the Uif bands have a supracellular organization that is contiguous across cells. Although transitions from high to low Uif levels can happen at cell-cell boundaries (Fig. 5D insets, white dots), many cells show discrete regions of high and low Uif levels on a single apical surface (Fig. 5D insets, white arrowheads). The biological significance of these Uif bands will be investigated in future studies, but as the existence of Uif bands is not obvious in conventional confocal projections (Fig. 5A) or slices (Fig. 5B) and has not previously been recognized despite the common use of Uif as a marker, these findings highlight the benefits of using QuBiT to analyze biological tubes.
Src42 acts in conjunction with other systems to regulate cell orientation
To test the ability of QuBiT to analyze mutants with altered tracheal development, we examined Src4226-1 (Fig. 6A-G), the only mutant for which apical cell orientation has been analyzed by independent software/workflows (Förster and Luschnig, 2012; Nelson et al., 2012). Analysis using QuBiT confirmed that, compared with WT, Src42 DTs are 24% shorter (P<10−3, Fig. 6B) and 13% larger in diameter (P=0.01, Fig. 6D). DT cell apical surface area is reduced by ∼28% (P<10−30, Fig. 6E), but cell apical aspect ratio is largely unchanged (Fig. 6F). Importantly, the cell apical orientations determined by QuBiT for Src4226-1 in segment 8 are in good agreement with previous studies by Förster and Luschnig (2012) and Nelson et al. (2012) (this report: 60±23°; Nelson et al., 2012: 66±28°; Förster and Luschnig, 2012: 73±12°) and are oriented more circumferentially than in WT (+13±36°; +24±36°; and +11±21°, respectively). However, the order-of-magnitude higher throughput of QuBiT extends the previous analyses by showing that, although the intrasegment orientation gradient appears to be dampened, the overall A-P cell apical orientation gradient is still present in Src4226-1 mutants (Fig. 6G). As loss of Src42 does not ‘lock’ all tracheal cell apical surfaces in a strongly orthogonal orientation, Src42 is unlikely to be a ‘master regulator’ of tracheal cell apical orientation. Rather, normal activity of Src42 biases the orientation of cell apical surfaces in which orientation is predominantly determined by a Src42-independent mechanism.
The Na+/K+ ATPase mutant disrupts localization of the apical marker Uif
To test the ability of QuBiT to follow tubes that have significant bends, and to use QuBiT to analyze an uncharacterized mutation, we examined a new CRISPR mutation in the Na+/K+ ATPase α-subunit, ATPα (Figs 5E-H, 6H-K). In Drosophila, the Na+/K+ ATPase has an ion transport-independent function in SJ formation (Paul et al., 2007) and a redundant role in apical/basal polarity (Laprise et al., 2010). This new allele, ATPαLNFS, results from a frame-shift mutation in the exon that encodes the long isoform of ATPα, which is thought to contain the SJ activity (Paul et al., 2007). Confirming this hypothesis, SJs are disrupted in the ATPαLNFS homozygotes, as evidenced by reduced Kune levels and failure of Kune to localize to the apical region of the lateral membrane (Fig. 5G-G′). As expected for a SJ mutation, ATPαLNFS increases tracheal tube length up to 30% (P<0.01, Fig. 6I), resulting in the tubes buckling and following a convoluted path between transverse connectives (Fig. 6H). Despite the tortuous paths, QuBiT is able to calculate tube centerlines and parameters. Tube cross-sectional area is generally increased compared with WT (P=0.01), but the diameter at fusion cells is not increased, conferring a saw-tooth-like appearance to the cross-sectional area graph (Fig. 6K) that would be expected from the sausage-link appearance of many SJ mutant trachea, including ATPα and wurst (Paul and Beitel, 2003).
Unrolling of ATPαLNFS mutant tubes reveals effects of Na+/K+ ATPase mutations that had not been previously observed. Uif staining in ATPαLNFS mutants was markedly different from that in WT: whereas in WT there are occasional areas where Uif is locally reduced, Uif-deficient regions in ATPαLNFS mutants were larger and more frequent (Fig. 5F, white arrowheads). Interestingly, these regions correspond to the surfaces of the trachea that were furthest away from the outer epidermis of the embryo, were confined to regions between the brighter Uif bands (yellow arrowheads) and did not coincide with any similarly irregular patterning of Kune (Fig. 5H). The disruption of the apical surface marker Uif in an ATPα single mutant extends our previous results regarding the role of the Na+/K+ ATPase in polarity (Laprise et al., 2009) and suggests that the polarity function of the Na+/K+ ATPase is not entirely redundant with that of Yurt. These finding also further illustrate that the tube unrolling function of QuBiT, a function not available in general epithelial analysis packages, can help identify patterns that are not readily evident from conventional confocal slices or 3D projections (e.g. Fig. 5E,F′).
QuBiT can analyze larger diameter tubes
To demonstrate that QuBiT can analyze tubes of different sizes, shapes and architectures, we examined the Drosophila embryonic stage 16 hindgut (partial max projection in Fig. 7A; 3D reconstruction in Fig. 7A′), which consists of a monolayer of cells like the trachea but can have a larger diameter (up to 14 µm and 18 cells per cross-section versus 5 µm and four cells, respectively). As with the trachea, we observed banding of Uif (Fig. 7B, yellow arrowheads), which is also visible in the tube cross-section (Fig. 7B′, yellow arrowheads). However, unlike the trachea, the intensity of the banding varied strongly along the tube length. In addition to the banding, the unrolled projection revealed two narrow stripes of missing Uif (white arrowheads in Fig. 7B-D) that coincide with boundary cells (red boxes in Fig. 7B-D, inset in D). The absence of Uif in boundary cells is unexpected given that boundary cells show elevated levels of the apical marker Crumbs (Kumichel and Knust, 2014). Thus, QuBiT can analyze both small and larger diameter tubes, and the use of QuBiT again revealed unexpected expression patterns of an apical marker.
QuBiT can analyze tubes with suboptimal apical staining
To test the ability of QuBiT to process low-quality data such as might be present in mutant or diseased tissue, we analyzed the Drosophila salivary gland (Fig. 7E-H′; partial max projection in Fig. 7E; 3D reconstruction in Fig. 7E′). Uif staining of the salivary gland apical surface is more sparse and patchy (Fig. 7F) than in the trachea, with significant sections of the circumference nearing background levels of Uif (Fig. 7F′, white arrowheads), but QuBiT is nonetheless able to render the apical surface and junctions (Fig. 7G) and overlay them (Fig. 7H, representative cross-section in Fig. 7H). Similarly, given that either physical tissue sections or optical confocal sections may fail to encompass an entire tube, we asked whether QuBiT could still segment and analyze the rest of the tube. As seen in Fig. 5B, areas of trachea tube that go beyond the bounds of a stack are projected as the absence of signal (Fig. 5B, red arrow), but the analysis of the rest of the tube proceeds. Thus, QuBiT can process sub-optimal image data that has either poor or incomplete signal.
We created QuBiT to quantify epithelial tube and cell apical properties in 3D, and to unroll tubes into 2D projections to enable analyses that are either difficult or impossible to do in 3D. The power of this approach is demonstrated by the fact that, although the Drosophila tracheal system is one of the most studied epithelial tube systems, QuBiT revealed novel findings. The ability of QuBiT to analyze tubes of different sizes and its robustness to suboptimal image data is demonstrated by our use of QuBiT to analyze the Drosophila hindgut and salivary gland. As these analyses again revealed previously unidentified biological phenomena, it is likely that there is much to be found by using QuBiT to analyze tubes in other systems.
MATERIALS AND METHODS
Flies used in this study were: w1118 (M. A. Krasnow, Stanford University, CA, USA), Src4226-1 (Takahashi et al., 2005), spider-GFP (Morin et al., 2001), and ATPαLNFS (see below).
Immunofluorescence and imaging
Drosophila embryos were fixed using the protocol in Samakovlis et al. (1996). Primary antibodies used were: guinea pig anti-Uif (1:750; Zhang and Ward, 2009); rat anti-Uif (1:500; this work, made per guinea pig anti-Uif, signal matches guinea pig anti-Uif), guinea pig anti-Verm (1:750; Luschnig et al., 2006) and rabbit anti-Kune (1:750; Nelson et al., 2010). Secondary antibodies, all from Thermo Fisher Scientific, were used at 1:1000 and were: Alexafluors 488 (A11006, A11073), 568 (A11011, A11075) and 647 (A21244). Confocal images were obtained using a Leica TCS SP5 confocal microscope with 40× and 100× objectives. Images were processed using Volocity, Fiji and Adobe Illustrator.
We used CRISPR to mutate the first exon of ATPα and create a frameshift eliminating the long ATPα isoform, while leaving intact the short isoform that has a separate transcriptional start site before the third exon. We assembled the repair template with ∼1 kb homology arms based on the method used in Bruckner et al. (2017) using Gibson assembly (Wang et al., 2015) of GBlocks fragments (Integrated DNA Technologies). The relevant sequence features and primer sequences are described below.
ATPα long N-terminus frame shift (LNFS) amino acid and corresponding DNA sequence
Amino acid sequence: start-MALRSDYEHGRADSFRVATVIATDDDKSHSRWPVQVKTQDAGQS* (underline indicates a frameshift; bold indicates a point mutation; * indicates the stop codon). Corresponding DNA sequence: ATGGCGTTAAGGTCGGATTACGAGCATGGGAGGGCGGACTCTTTCCGAGTCGCCACCGTTATAGCCACCGATGACGACAATCGCACAGCCGATGGCCA GTACAAGTCAAGACGCAAGATGCCGGCCAAAGTTAA (underline indicates the mutated NGG PAM; bold indicates a mutation). Primer sequences for left homology arm: forward, 5′-ggctatcgcttgctggaaag-3′; reverse, 5′-aatcgaaaggggtcggtatt-3′. Primer sequences for right homology arm: forward, 5′-gcgctttactaacccacttc-3′; reverse, 5′-acgtaaatgaattgaggtgc-3′.
QuBiT – code and sample datasets
The MATLAB based code for QuBIT, a user manual, and sample datasets are available online at http://github.com/gjbeitel/QuBiT.
We thank Rob Ward for anti-Uif anti-sera, Richard Carthew for helpful discussions, and J. Hornick and the Northwestern University Biological Imaging Facility for imaging support.
Conceptualization: R.Y., M.M., G.J.B.; Methodology: R.Y., E.L., Y.-J.K., M.M., G.J.B.; Software: R.Y., E.L., M.M.; Validation: R.Y., Y.-J.K., M.M., G.J.B.; Formal analysis: R.Y., M.M., G.J.B.; Investigation: R.Y., Y.-J.K., M.M., G.J.B.; Resources: G.J.B.; Data curation: R.Y.; Writing - original draft: R.Y., G.J.B.; Writing - review & editing: R.Y., E.L., M.M., G.J.B.; Visualization: R.Y., G.J.B.; Supervision: G.J.B.; Project administration: G.J.B.; Funding acquisition: M.M., G.J.B.
This work was supported by National Institutes of Health (RO1GM108964 to G.J.B.) and the Simons Foundation (597491 to M.M.). M.M. is a Simons Foundation Investigator. Deposited in PMC for release after 12 months.
The authors declare no competing or financial interests.