Whole-body clonal mapping identifies giant dominant clones in zebrafish skin epidermis

The growth of vertebrate animals during post-embryonic stages requires massive proliferation of stem cell populations that were specified during early development. As each organ or tissue-specific stem cell has the ability to produce progeny cells, intricate control mechanisms are needed to coordinate and synchronize their individual contributions or so-called ‘clonal growth’. Despite extensive application of lineage tracing and clonal analyses to study stem cell dynamics during early development and adult hemostatic conditions (Byrd et al., 2019; Chatzeli and Simons, 2020; Gonzales and Fuchs, 2017), little is known about the clonal behavior that enables extensive post-embryonic growth of tissues and organs. In particular, it is almost completely unknown how individual stem cells or progenitor cells determine their extent and mode of contribution to larger populations.

Skin epidermis is a relatively simple model with which to study the mechanisms of clonal growth. In the stratified structure of skin, basal epithelial cells (BECs) proliferate to generate more BECs or differentiate to generate suprabasal and superficial epithelial cells (Blanpain and Simons, 2013). Using mouse tail and paw epidermis as a model, a pioneering study from Dekoninck et al. demonstrated that skin stem cells are more uniform, in terms of transcriptional profiles, than their adult counterparts (Dekoninck et al., 2020). Intriguingly, the growth of each clone is determined by a fixed imbalance of self-renewing divisions and an ever-decreasing proliferation rate, while the shapes of individual clones appear to be locally regulated by the underlying collagen fiber orientation (Dekoninck et al., 2020). Yet, owing to technical challenges with simultaneous monitoring of clonal behavior over an extensive area of the mouse body surface, it has not yet been determined whether dominant clones (defined as clones with significantly more progeny cells than competent neighbors) may arise to support vertebrate skin growth and development.

The zebrafish skin epidermis offers a well-characterized system for studying clonal dynamism in the BEC population. Similar to those in the mouse model, BECs in zebrafish skin function as a stem cell pool, generating a continuous supply of epithelial cells during development and in adult homeostatic tissues (Guzman et al., 2013; Lee et al., 2014). Also like their mouse counterparts, epithelial cells in all strata of the zebrafish skin are derived from BECs during post-embryonic growth (Lee et al., 2014). One major difference between mammalian and zebrafish models is that fish skin has no keratinized, non-transparent dead-cell layer on the outermost body surface. This key feature and a relatively flat body shape make zebrafish uniquely amenable to high-resolution intravital imaging of large numbers of skin cells over an extensive body area (Chen et al., 2016).

In this study, we sought to delineate clonal growth behavior over the entire zebrafish body surface during the massive expansion of skin epidermis that occurs in post-embryonic growth periods. We created a multicolor cell tagging tool and a centimeter-scale imaging platform to capture the growth and location of each stochastically tagged BEC over the course of a 28-fold body surface expansion. Through same-animal whole-body monitoring, we found that ∼86% of BECs survive and grow, regardless of their anatomical location. Unexpectedly, we also found that regional body surface areas proportionally define the size and shape of individual BEC clones, thus explaining variations in clonal heterogeneity that may seem random when monitoring small body surface areas. Systematic mapping of 742 BEC-derived clones further revealed that clonal dominance is a key mechanism supporting zebrafish skin expansion, a recurring theme in vertebrate organ development (Gupta and Poss, 2012; Nguyen et al., 2017; Rios et al., 2014). Finally, we showed that transient perturbation of a growth-induced ECM component in the BEC population is sufficient to decouple clonal growth from surface growth. Overall, this study identifies dominant clone growth as an adaptive stem cell behavior that supports massive expansion of vertebrate skin epidermis during post-embryonic development.

During post-embryonic growth, the zebrafish skin epidermis undergoes a gradual transition from a simple bi-layered to a multi-layered structure (Guzman et al., 2013). We took advantage of the 3.9 kb krt19 promoter, which can ectopically drive transgene expression in the entire skin epidermis, to create a panel of transgenic lines based on the Brainbow 1.0L cassette (Livet et al., 2007). After an extensive screen, we identified one line, Tg(krt19:Brainbow)^as65, that had labeling throughout the whole skin epidermis, including both the BEC and the superficial epithelial cell (SEC) populations over different developmental stages (Fig. 1A-C). To enable controlled, multicolor tagging of individual skin cells, we then conducted another screen for inducible Cre lines based on the ERT2CreERT2 design (Matsuda and Cepko, 2007) (Fig. 1A). By pairing the Cre lines with the Floxed reporter line, we identified one stable Cre line, Tg(krt19:ERT2CreERT2)^as66, that can be used for stringent control of Cre recombinase activity. In the double-transgenic progeny, we found that the duration of tamoxifen (Tam) treatment quantitatively determined the number of tagged BECs (Fig. S1). In contrast, we detected no leaky background recombination in the absence of Tam in any of the examined animals (10/10; Fig. S1B), a crucial prerequisite for long-term cell tracking studies. As we intended to use this Brainbow technology-based, cell-tagging system to visualize clonal growth behavior in the BEC population, we hereafter refer to the system as basebow.

To determine whether basebow may enable long-term monitoring of the BEC growth, we first induced Cre activity for 4 or 8 h at 4 days post-fertilization (dpf), and then we imaged animal body surfaces at 8, 21, 56 and 90 dpf (Fig. 1D). Of note, basebow labeled only a small proportion of SECs at 8 dpf upon a prolonged activation of Cre (i.e. 4 to 8 h of Tam treatment; Fig. S1C,D). At this developmental stage, zebrafish skin epidermis is bi-layered, and the outermost SEC population is entirely shed by 32 dpf (Lee et al., 2014). Thus, each remaining skin clone at later stages can be used as a proxy for clones derived from the 8 dpf BEC pool. By imaging singly labeled BECs from the larval to adult stages, we found that individual BECs readily undergo clonal expansion (Fig. 1E), consistent with their role as the epidermal stem cells in the skin (Lee et al., 2014). Notably, the overall clone size dramatically increased, while the average cell size within each clone was moderately reduced (16% decrease at 21 dpf; Fig. 1F). We further determined that most, if not all, BEC-derived clones remained rooted in the basal layer upon expansion, as we failed to detect any skin clones that were exclusively superficial from orthogonal views of confocal images (n=0/78; Fig. 1E). Thus, we concluded that direct measurements of clone size in 2D space should provide conservative, yet reasonable, estimates of the extent of lateral expansion and the cell numbers in respective clones. Taken together, these findings led us to conclude that basebow is suitable for visualizing clonal growth behavior in all strata of the skin epidermis, and it might allow large-scale, long-term mapping of each BEC-derived clone in intact and live zebrafish.

To unambiguously define the growth of each tagged BEC, we transiently incubated basebow animals in Tam for only 2 h at 4 dpf. This treatment scheme not only restricted basebow-mediated cell tagging to the BEC population (Fig. S1D), but it also allowed us to achieve sparse cell tagging over the entire skin epidermis (24.7 BECs per animal; Fig. 2A-D). Thus, each potential clone would be sufficiently isolated from the others at later stages, allowing for more accurate identification and quantification. To monitor the fate and growth of each tagged BEC, we performed same-animal whole-body imaging at two time points: 7-8 dpf and 35-36 dpf. Of note, during this 4-week period, the zebrafish body length and body surface undergo drastic expansion, as measured by Standard length (Parichy et al., 2009) and Body surface area (364% length increase and 28-fold expansion, respectively; Fig. 2E,F). To enable whole-body monitoring of all clones on a large animal, we set up a macroscopic imaging platform with an epifluorescence microscope and a highly sensitive sCMOS camera. We then conducted z-scans of the entire volume of the zebrafish body and performed a post hoc parallax correction, 2D-deconvolution and maximal projection to compare BEC growth on different body regions (Fig. 2G,H; see Materials and Methods). Of note, although the platform enables fluorescence imaging on a centimeter scale, the system falls short of being able to differentiate fluorescence signals at varying levels, which is a requirement for tracking secondary and/or tertiary colors generated with Brainbow technology. Remarkably, we found that ∼86% of the 8 dpf BECs are able to survive and grow (n=191/222 from 9 animals; Fig. 2I,J). On average, each BEC-derived clone expanded to occupy about 0.013 mm² of the body surface area (Fig. 2K). To exclude the possibility that clones may be occasionally derived from unlabeled cells, we performed whole-body imaging on 8 and 53 dpf basebow zebrafish with no prior exposure to Tam. We detected no trace of spontaneously labeled cells or clones throughout the entire skin epidermis of all examined animals (n=0/18 and 0/24, respectively; Fig. S2). Thus, based on same-animal whole-body monitoring of 222 stochastically labeled BECs from larval to early juvenile stages, we concluded that nearly all BECs survive. Furthermore, maximizing the survival rate of each competent stem cell and/or progenitor cell in a population may serve as a key cellular strategy to support massive expansion of skin surfaces during post-embryonic growth.

To determine whether clonal growth features, such as size and shape, may be spatially coordinated on a large scale, we conducted a whole-body, head to tail scan of 742 BEC clones across different body regions (n=37 animals; Fig. 3A,B). Of note, we captured clonal growth features from both sides of the animal to double the number of examined clones (Fig. S3A,B). We detected no notable differences between the left-side BEC clones and the right-side clones on the same animals (Fig. S3C-E). Intriguingly, although all BECs were singly labeled at 8 dpf along the anterior-posterior axis of the fish body, we found that the resulting BEC clones on the anterior part of the fish body surface were on average 1.7-fold larger and less rounded than those on the posterior part, as determined by the two-dimensional size of the clone, major axis length and ratio of each clone area to its squared convex perimeter (Fig. 3C-E). Of note, a vertical line drawn at the anterior of the anal fin was taken as a midline dividing the fish body surface into anterior and posterior halves (AS versus PS; Fig. 3B). Intriguingly, despite local modulation of BEC clone morphologies by scale structures (Fig. 3F), the orientations of anterior clones on a whole-animal scale were markedly biased toward the dorsal-ventral axis of the fish body; posterior clones were more isotropic, with no specific shape orientation (Fig. 3G). To determine whether regional differences in body surface areas may affect clonal growth, we further compared the anterior half of the body surface area with the posterior half. Intriguingly, we found that the AS is about 1.8-fold larger than the PS (Fig. 3H), a fold difference that is rather close to the average fold difference in size of clones found in each half (1.7-fold; Fig. 3C). This similarity might suggest that the average growth of each BEC clone on a regional scale could be proportionally coupled to the extent of skin surface expansion on distinct body regions. Of note, we could not have detected this coupling if the clonal growth behavior were examined only at the scale of local body regions.

Because some BEC clones appeared to have undergone extensive growth (Fig. 3C), we wanted to determine whether ‘elite’ dominant clones may exist, making outsized contributions to skin surface expansion during growth. To define a priori cutoff criteria for dominant clones, we first examined total number of BECs in the zebrafish skin epidermis at 8 dpf, the time point at which each labeled BEC remained as a single cell. Using a transgenic line that labels only the BEC population [Tg(krt19:H2A-mCherry)^as54; Fig. S4], we determined that an 8 dpf zebrafish larvae hosts about 3000 BECs on one side of its body (3109±117; n=8 animals; Fig. 4A,B). Notably, the BEC density remains more or less constant on different body surface areas during post-embryonic growth (Fig. 4C,D). Thus, when a clone grows proportionally with the fish body surface area expansion, an ordinary BEC clone is expected to occupy 0.03% of the body surface area in 2D space (i.e. 1/3109=0.03%; assuming all 8 dpf basal cells can grow and make the same contributions). We then surveyed the above-mentioned BEC clones for their respective ratio of body surface area coverage. Remarkably, we found that during post-embryonic skin growth, a few dominant clones readily emerge. Some dominant clones occupy up to 0.6% of the animal body surface, which is about 20-fold larger than expected (0.03%) for ordinary clones that grow proportionally (Fig. 4E). Among a total of 742 BEC clones, we identified 3.6% giant dominant clones (cutoff ratio: 0.2%), 7.7% dominant clones (ratios of 0.1%-0.2%), 20% near-dominant clones (ratios of 0.05%-0.1%) and 67% ordinary clones (ratios of 0.007%-0.05%; Fig. 4F). Of note, we set a specific size cutoff at 0.1% for dominant clones, which is about threefold higher than the theoretical ratio of an ordinary clone (0.03%). BEC clones that occupied less than 0.007% of the body surface or were smaller than 0.002 mm² were not included in our semi-automated identification pipeline (Fig. S5A,B; see Materials and Methods). Intriguingly, we found that ‘giant dominant clones’ grow exclusively on the anterior region of the fish body surface (27/27; AS in Fig. 4G,H). We failed to detect any clones that could occupy more than 0.2% of the body surface area on the posterior region (0/102; PS in Fig. 4G,H). Similarly, as many as 96% of the ‘dominant clones’ were found on the AS region (54/56; Fig. 4I). Of note, we determined clone size based on 2D maximal-projection images, which are likely to underestimate actual clone sizes due to large variations in z-position over the wide body surface area (Fig. S5C). Therefore, there may be more dominant BEC clones to support skin growth than we estimate here. To determine the expansion dynamics that may lead to the appearance of a dominant clone, we performed same-animal same-clone tracking at 21, 25 and 28 dpf (Fig. 5A,B). Intriguingly, although individual BEC clones appeared to grow at distinct rates, we determined that dominant clones expand at a faster pace on average than non-dominant clones (Fig. 5C,D). Of note, we failed to detect any clearly split or shrunken clones in the examined cases (n=17). Altogether, the whole-body monitoring of single cell growth at a centimeter scale allowed us to conclude that the extent to which a BEC can grow is varied yet spatially constrained. Furthermore, the dominant clone populations are not anatomically widespread but confined to specific body regions.

As the extent of clonal growth is proportionally coupled to the body surface area of distinct regions (Fig. 3), we speculated that there might exist specific molecular regulators that enable this coupling to occur. From a RT-qPCR analysis of six candidate genes (i.e. lef1, fgf20a, pola2, ddx52, lamb1a and lamb1b) that might affect clonal growth of the BEC population by various mechanisms (Armstrong et al., 2017; Chen et al., 2015; Poss et al., 2000; Tseng et al., 2021; Wang et al., 2019; Whitehead et al., 2005), we found that lamb1a displays differential expression between AS and PS skin tissues (Fig. S6). The expression in AS is about 2.2-fold higher than it is in PS of adult tissues and 1.3-fold higher in 35 dpf juvenile tissues (Fig. S6F). As it encodes an extracellular matrix component, lamb1a is transiently induced in the BEC population of adult regenerating tailfin tissues; the protein is ectopically deposited to the basement membrane to mediate epithelial-mesenchymal interactions that are essential for jumpstarting tailfin regeneration (Chen et al., 2015). Intriguingly, we determined that lamb1a expression is also elevated in juvenile developing skin tissues (an increase of 2.7-fold in AS skin tissues when comparing 35 dpf juveniles to adults; Fig. S6F), suggestive of a role in supporting post-embryonic skin growth. Thus, we further hypothesized that Lamb1a activity may determine the extent and mode of clonal growth in the BEC population. To perturb Lamb1a activity in a spatiotemporally controlled manner, we took advantage of a temperature-sensitive allele, lamb1a^sde1 (Chen et al., 2015), to manipulate Lamb1a localization and function in intact live basebow zebrafish (Fig. 6A). Upon shifting the animals from 28°C to 34°C, homozygous lamb1a^sde1 allele carriers had Laminin protein sequestered in skin BECs (sde1, Fig. 6B,C), preventing its proper deposition at the basement membrane (Chen et al., 2015). We then examined the BEC clonal growth behavior after a 2-week heat treatment. We found that clone number, clone size and major axis length were significantly reduced in the homozygous sde1 mutants compared with their heterozygous siblings (sde1 versus sde1/+; deceases of 48%, 33% and 22%, respectively; Fig. 6D-G). Correspondingly, the homozygous mutants had fewer dominant clones within the population (Fig. 6H). Despite the profound effect on the clonal growth profile, direct inactivation of Lamb1a activity appeared to have no immediate impact on cell proliferation, as determined by the number of EdU and P63 double-positive cells (Fig. S7). Taken together, we speculated that the Lamb1a effect on clonal growth behavior is likely to be specific, as the long-term perturbation of Lamb1a activity has no notable influence on fish body length or surface area (Fig. 6I,J). Of note, we detected the fin atrophy phenotype in homozygous sde1 mutants upon Lamb1a inactivation (Fig. 6K), consistent with the reported role of Lamb1a in maintaining the morphology of tailfin tissues (Chen et al., 2015). Altogether, these findings led us to conclude that growth-induced, differential expression of Lamb1a on distinct body compartments facilitates local clonal growth in a cell-autonomous manner. Thus, the formation of dominant skin clones is promoted on specific body regions that require more surface coverage.

Here, we developed a cell tagging and imaging system to monitor clonal growth behavior over the entire skin epidermis of a zebrafish during post-embryonic growth. By tracking the clones of single skin stem cells during the 28-fold expansion of body surface area, we determined their overall survival rate, the extent and mode of their growth, and molecular, local and global features of their expansion. We propose that clonal dominance is an integral, adaptive stem cell mechanism for supporting post-embryonic skin growth in vertebrate animals. To the best of our knowledge, the macroscopic imaging platform described here may represent the first successful attempt to capture clonal growth behavior of single stem cells on a whole-body, centimeter scale.

Through same-animal whole-body monitoring of clonal growth, we determined that: (1) majority of skin stem cells survive and clonally expand during post-embryonic growth; (2) within each expanding clone, the cell size is only moderately changed, regardless of clone size or developmental stage; and (3) local environmental cues (i.e. fish scales) shape the orientations of expanding clones. Intriguingly, these key findings are consistent with clonal growth behaviors discovered in the mouse tail skin epidermis. By imaging a skin area of ∼1 mm² on a mouse tail from P4 to P60, Dekoninck et al. (2020) identified that: (1) skin clone numbers are stable during post-embryonic growth, with the gradually increasing number of cells within each clone trailing the expansion of the tail surface; (2) average skin cell sizes are relatively constant, regardless of developmental stage; and (3) local orientation of collagen fibers evidently affects the shape of the clones. Thus, based on the overall similarities of clonal growth behaviors between the zebrafish and mouse skin epidermis, we speculate that dominant clones and/or giant dominant clones may also exist in mammalian tissues to support the expansion of the skin epidermis. These elite skin clones are not widespread but are region specific, and they may only be found in distinct body regions where surface areas expand most during post-embryonic development.

Our findings highlight clonal dominance as a recurring mechanism for supporting the development of vertebrate organs, including the heart, skeletal muscles and the mammary gland (Gupta and Poss, 2012; Nguyen et al., 2017; Rios et al., 2014). Why some stem cells or progenitor cells become dominant is an intriguing question. Cell competition-based models, either active or passive, have been used to explain why dominant clones might emerge in a population (Amoyel and Bach, 2014; Claveria et al., 2013). Notably, most cell competition models include the inherent assumption that many less-competent neighbors will be eliminated by the dominant clone. In contrast to this assumption, we provide evidence that the majority of skin stem cells survive and thrive during post-embryonic growth, regardless of their dominance. Thus, instead of adapting a competitive elimination strategy, individual stem cells may work collaboratively to ensure robust growth of tissues and organs on a large scale, which may give rise to large-sized clones in either a stochastic or an opportunistic manner (Krieger and Simons, 2015).

Here, we provide evidence that an ECM component, Lamb1a, is transcriptionally upregulated in skin tissues during post-embryonic growth, and its expression level appears to be spatially correlated with the extent and mode of clonal growth on distinct body regions. Although the specific cellular mechanisms that enable Lamb1a to regulate clonal growth remain to be determined (Diaz de la Loza et al., 2017; Urbano et al., 2009), our findings support the idea that there is active crosstalk between the BECs and the underlying dermis; our study results are also well aligned with the emerging theme that ECM components are dynamically regulated and have diverse functions during development (Derrick and Noel, 2021; Keeley et al., 2020; Matsubayashi et al., 2020). Further studies will be needed to determine how local BECs are able to perceive their location and how proper laminin levels are maintained in different body regions to support appropriate clonal growth. It is known that a network of growth factors, such as Fgfs and Wnts, can regulate laminin expression during tissue growth (Chen et al., 2015; Nagendran et al., 2015), and an intriguing study from Priya et al. demonstrated that local tension heterogeneity can modulate organ-scale patterning (and even cell-fate choices) during cardiac trabeculation in zebrafish (Priya et al., 2020). In particular, it will be interesting to determine whether dynamic changes in local skin tension or intra-tissue tension (Ning et al., 2021) may similarly modulate local expression of lamb1a. If so, this mechanism could instruct context-specific production of laminins to regulate clonal growth at a large scale.

In summary, this study reveals how clonal heterogeneity may arise from a seemingly homogeneous stem cell population in a vertebrate tissue. We identify bona fide dominant skin stem cell clones and where to find them, and we determine how their growth may be coordinated on a population scale. Our findings reveal an intricate, global control of local clonal growth behavior during large-scale tissue morphogenesis.

Tg(krt19:Brainbow 1.0L)^as65 was generated with a transgenic construct consisting of a 3.9 kb krt19 promoter (Lee et al., 2014) and the Brainbow 1.0L sequence (Livet et al., 2007). Tg(krt19:ERT2CreERT2)^as66 was created with a transgenic construct consisting of a 4 kb krt19 promoter and the ERT2CreERT2 sequence (Matsuda and Cepko, 2007). In Tg(krt19:ERT2CreERT2)^as66, a cmlc2:mCherry fragment flanked with I-SceI sites was co-injected as a selection marker. Tg(krt19:H2A-mCherry)^as54 was created with the krt19 promoter and the histone 2A-tagged mCherry sequence for determining the in toto BEC number. All constructs were flanked with I-SceI sites to facilitate transgenesis. The lamb1a (sde1) and pola2 (mem) ts mutants were as described previously (Chen et al., 2015; Wang et al., 2019). Zebrafish larvae were fed with paramecia starting at 8 dpf. Animal density was maintained at one fish per 3 l in most experiments, unless otherwise specified. For Lamb1a inactivation experiments, sde1 and sde1/+ animals were transferred from 25°C to 34°C for the indicated periods. All animal experiments were approved by the Institutional Animal Care and Utilization Committee (IACUC) at Academia Sinica.

Cre activity was transiently induced in 4 dpf embryos by tamoxifen (2 μM; Sigma, T5648) in all experiments. The induction time was 2 h, unless otherwise specified.

Larvae were sedated with tricaine (0.4 mg/ml) and mounted in 1.5% low melting point agarose prior to imaging. For same-animal tracking experiments, individuals were revived by transferring to aquarium water with an air pump. All images of 8 dpf larvae were captured using a Leica SP8 upright confocal microscope with a 25× water dipping lens (25×/0.95 HCXIRAPO). For imaging basebow, only the 448 nm laser line was used for excitation. Emission filters were set up as follows: (1) 455-480 nm bandpass (blue channel); (2) 510-540 nm bandpass (green channel); and (3) 600-700 nm bandpass (red channel). For imaging Tg(krt19:H2A-mCherry), the 552 nm laser was used. The emission filter was 590-700 nm bandpass (red channel). For whole-mount imaging of juvenile and adult basebow, animals were anesthetized in fish water containing tricaine (0.2 mg/ml) before imaging with an epifluorescence microscope (Leica M205). Clone images were captured using a monochrome camera with a highly sensitive sCMOS sensor (Leica DFC9000T). For histology, images were captured using SP8 confocal with either a 63× oil lens (HC PL APO 63×/1.40 Oil CS2, WD 0.14 mm) or a 20× lens (Leica SP8; HC PL APO 20×/0.75 CS2, WD 0.62 mm).

Fish body surface and standard length were measured manually in FIJI (ImageJ) (Schindelin et al., 2012). For whole-mount epifluorescence images, signals from a z-stack (FOV: 6.65 mm×6.65 mm; system optimized mode) were collected, parallax corrected and processed using the Extended Depth of Focus (EDOF) function of LAS X software (Leica). Maximum intensity projection images were then analyzed for clone size (i.e. area and major axis length) and shape factors (i.e. circularity and orientation) using FIJI. For identifying individual BEC clones in an unambiguous manner, clone images captured by green and blue channels were analyzed separately. In brief, images of known dimensions (6.65 mm; 2048×2048 pixels) were converted to binary images. Individual clones larger than 0.002 mm² were identified using LoG detector (estimated blob diameter: 0.4 mm; 0.01 as threshold with sub-pixel localization) in TrackMate (Tinevez et al., 2017). Each TrackMate-identified clone was visually inspected for incorrect labeling prior to further analyses. Wand tool in FIJI was used to outline individual clones for computing size and shape factors. For counting total BECs in the Tg (krt19:H2A-mCherry), TrackMate was used (estimated blob diameter: 13 µm; 1 as threshold with sub-pixel localization). Only BEC clones covering the trunk body surface area were included for quantification (Fig. 3). Whole-animal images were stitched using Image Composite Editor (Microsoft Research). For quantitative analysis, non-stitched, maximum intensity projection images were processed and analyzed. For optical section views of BECs and SECs, Bitplane Imaris 9.6 software (Oxford Instruments) was used.

Tissues were fixed with 4% paraformaldehyde (PFA) in phosphate-buffered saline (PBS) at 4°C overnight. Fixed samples were washed in fish fixation buffer (0.4% sucrose in PBS) three times, for 5 min each, before being embedded in 1.5% agarose with 5% sucrose. Samples were then incubated in 30% sucrose for overnight at 4°C. Frozen blocks were sectioned at 16 µm with a cryostat for antibody staining. The primary antibodies used were: mouse anti-p63 (1:250, Biocare Medical, CM 163A), rabbit anti-mCherry (1:250, GeneTex, GTX128508) and rabbit anti-laminins (1:250, Sigma, L9393). The secondary antibodies were: Alexa Fluor 488 anti-mouse (1:400, Invitrogen, A-11001), Alexa Fluor 594 anti-mouse (1:400, Invitrogen, A-11005) and Alexa Fluor 594 anti-rabbit (1:400, Invitrogen, A-11037).

Larvae at 10 dpf were immersed in EdU solution (5 mM in 5% DMSO; ThermoFisher, A10044) at 26°C for 1 h before fixation with 4% PFA for 2 h. Samples were then washed with PBST three times for 5 min each, before dehydration with methanol through a 25%, 50% and 75% series (5 min each), and stored in 100% methanol overnight at −20°C. The next day, samples were rehydrated through the reverse methanol series and then washed three times in PBST for 5 min each. After re-fixation with 4% PFA at room temperature for 20 min, larvae were washed three times in PBST for 5 min each. Next, larvae were incubated in EdU staining solution [1 mM CuSO₄, 50 mM ascorbic acid, 100 mM Tris buffer (pH 8.0)] and 10 µM Alexa Fluor 488 azide for 30 min at room temperature in dark. After three washes in PBST for 5 min each, larvae were incubated in blocking buffer [10% heat-inactivated newborn calf serum (NCS, Gibco, 26010066), 4% goat serum (Gibco, 16210072) and 0.1% DMSO in PBST] for 30 min at 37°C. Blocking buffer was drained before adding primary antibodies in NCS-PBST: mouse anti-p63 (1:250, Biocare Medical, CM 163A) against BECs for 3 h at 37°C. Samples were washed three times in PBST for 5 min each, before incubation with secondary antibody: Alexa Fluor 488 anti-mouse (1:400, Invitrogen, A-11001) for 1 h at 37°C. Nuclei were stained with DAPI in PBST, then washed three times in PBST for 5 min each. Samples were mounted in 1.5% low melting point agarose prior to imaging (Leica SP8; 25x/0.95 HCXIRAPO).

Skin tissues were collected from four individuals by scale-plucking under a dissecting microscope and homogenized in 1 ml Trizol (Sigma, Cat. No: T9424-200ML), using the TissueLyser II (Qiagen). cDNA was synthesized from 0.6 μg RNA using the SuperScript III First-Strand Synthesis System (Invitrogen, 18080051). qPCR analysis was performed with a Roche LightCycler 480 following the manufacturer's instructions. Primer sequences are listed in Table S1. Each sample was analyzed in biological quadruplicate and technical triplicate. Results were analyzed by the ΔΔCT method, using the level of actb cDNA as an internal control.

Statistical significance was analyzed by Prism 9.0 (GraphPad Software). All statistical values are displayed as mean±s.d. unless otherwise specified. A two-tailed Student's t-test was used for data with parametric distributions when distributions passed the D'Agostino-Pearson normality test. Otherwise, the two-tailed Mann–Whitney test was used for data with non-parametric distributions. Sample sizes and statistically significant differences are reported in the figure or the corresponding figure legends.

We thank Taiwan Zebrafish Core Facility (TZCAS; MOST 108-2319-B-400-002) for maintenance of zebrafish lines, the Chen laboratory members for comments on the manuscript, and Marcus J. Calkins for English editing and comments.

The peer review history is available online at https://journals.biologists.com/dev/article-lookup/doi/10.1242/dev.199669.