Identification and characterization of T reg–like cells in zebrafish

Studies of hematopoiesis in the zebrafish have defined major blood lineages. In the kidney marrow, the site of adult hematopoiesis, erythroid, lymphoid, and myeloid cells are present and can be separated based on light-scattering characteristics (Traver et al., 2003). However, subsets of cells within each of these lineages are poorly defined, primarily because of the lack of markers to distinguish a desired subset from other cells within the same lineage. In the lymphoid lineage, lck:EGFP and cd4-1:mCherry transgenes (Langenau et al., 2004; Dee et al., 2016) have been used to isolate bulk and CD4⁺ T cells, respectively, but thus far no markers have been established to further separate T cell sublineages. Indeed, it is not clear which sublineages of hematopoietic cells are present in zebrafish. Defining these sublineages in the zebrafish would further enable studies of developmental hematopoiesis, hematologic disease, and evolutionary immunology.

Among T cell subsets, regulatory T (T reg) cells are distinguished by their immunosuppressive activities. They are particularly important in maintaining peripheral tolerance by suppressing autoreactive T cells that escape negative selection in the thymus or inactivation in the periphery. A defining characteristic of T reg cells is their expression of the FOXP3 transcription factor. FOXP3 regulates expression of a suite of target genes, some of which are important for establishing the T reg lineage and others that are key mediators of T cell activation (Marson et al., 2007; Zheng et al., 2007; Birzele et al., 2011). FOXP3 function is required for T reg development, and its loss in humans and mice results in T reg deficiency, leading to severe congenital autoimmune reactivity and other immunopathologies (Bennett et al., 2001; Brunkow et al., 2001; Wildin et al., 2001). Depletion of T reg cells in adults likewise causes inappropriate immunoreactivity, which is characterized by severe lympho- and myeloproliferation accompanied by pervasive inflammation (Kim et al., 2007; Lahl et al., 2007).

Whereas in vitro studies have characterized the functional basis for T reg–mediated immune suppression, the mechanisms by which these cells exert their functions in vivo are poorly understood. Multiple mechanisms for T reg activity have been defined—production of immunosuppressive cytokines, cytolysis of target cells, cytokine deprivation or metabolic disruption of target cells, modulation of dendritic cell activity—but the relative importance of each in a variety of T reg functions has yet to be determined. Furthermore, the anatomical locations in which T reg cells mediate suppression in trans are only beginning to be explored. Liu et al. (2015) used high-resolution, spatially reconstructed histology to show that colocalization of T reg cells and target autoreactive T cells in the gut is important for T reg–mediated suppression. Although some progress has been made, these studies have been hampered by the inability to live-image cells in their native environments. In some studies, T reg cells have been adoptively transferred into mice and their interactions tracked (Mempel et al., 2006; Bauer et al., 2014), but to date there has been no way of visualizing the trafficking of these cells and their associations in a native setting. In the zebrafish, labeling bulk T cells with lck:EGFP has been instrumental in determining the locations of developmental hematopoiesis as well as providing a means to monitor acute T cell lymphoblastic leukemia (T-ALL) progression (Langenau et al., 2004; Frazer et al., 2009). More recently, cd4-1:mCherry transgenic zebrafish have been used to visualize this subset of T cells. In this study, cd4-1:mCherry-positive T cells from several anatomical sites were analyzed, and up-regulation of foxp3 expression was observed in cells from some tissues, suggesting that T reg–like cells might be present and enriched in certain tissues (Dee et al., 2016). Given the ability to directly visualize native cells in zebrafish, identification of T reg cells in zebrafish would provide a platform to further understand how these cells interact with other immune cells, traffic through tissues, and are recruited to cancers. Indeed, a preponderance of tumor-resident T reg cells is linked to cancer-driven immunosuppression and poor clinical outcomes (Roychoudhuri et al., 2015), and determining how these cells function in this milieu could further enhance the effectiveness of immune checkpoint therapies.

To address these issues, we sought to determine whether T reg cells were present in zebrafish. Using a FOXP3 orthologue as a marker, we identified CD4-enriched, mature lymphocytes with several hallmarks of T reg cells. Animals mutant for this FOXP3 orthologue display hematopoietic imbalances characterized by excess T lymphocytes. Furthermore, these mutants have a severe inflammatory defect that results from inappropriate lymphocyte activity. Together, these data indicate that T reg–like cells are present in zebrafish and are required for proper immunoregulation. The ability to isolate and modulate the function of T reg–like cells in zebrafish opens up new avenues for studying the normal regulation of immunity by T reg cells as well as the role of these cells in pathological conditions such as autoimmunity and cancer.

T reg cells are characterized by specific expression of the FOXP3 transcription factor. There are two FOXP3 orthologues in zebrafish, foxp3a and foxp3b, each of which is similarly diverged from mammalian FOXP3 (Fig. 1 A). To determine whether one or both of these could mark zebrafish T reg cells, if present, we performed quantitative RT-PCR (qRT-PCR) analysis of zebrafish whole kidney marrow (WKM), the site of hematopoiesis in adult zebrafish (Fig. 1 B). foxp3a was highly enriched in WKM, whereas foxp3b levels were only modestly higher in WKM compared with nonhematopoietic control fin tissue. For this reason, we constructed a foxp3a reporter construct that placed 2.4 kb of foxp3a promoter sequence upstream of the EGFP gene. Zebrafish embryos were injected with this construct, and a stable transgenic line was obtained.

We performed flow cytometry to determine cell subsets that express foxp3a:EGFP. Cells from WKM of Tg(foxp3a:EGFP) animals were first plotted according to side-scatter and EGFP positivity, then cells from the EGFP-positive population were plotted according to side- and forward-scatter (Fig. 1 C). Cells in the zebrafish WKM segregate into distinct blood populations based on light-scattering characteristics. Gates corresponding to erythroid, lymphoid, precursor, and myeloid-containing populations were assigned based on our own (Fig. S1 A) and previous (Traver et al., 2003; Lin et al., 2005) data. Most foxp3a:EGFP-positive cells were found in the lymphoid-containing gate. Using a different gating strategy, we examined whether foxp3a:EGFP-positive cells were present in other blood populations. We first plotted cells according to side- and forward-scatter, then from each blood population, we determined the percentage of cells that were EGFP-positive (Fig. 1 D and Fig. S1 A). The fraction of EGFP-positive cells in the lymphoid-containing gate represented a 132-fold increase in Tg(foxp3a:EGFP) animals compared with control animals. Nonlymphocytes that gated as EGFP-positive were present in wild-type AB animals; thus such cells in Tg(foxp3a:EGFP) animals were most likely background cells with intrinsic autofluorescence. To further investigate the nature of foxp3a:EGFP-positive cells from the lymphoid gate, we isolated these cells by FACS and viewed their morphology using May-Grünwald Giemsa staining. foxp3a:EGFP-positive cells had a lymphoid morphology similar to lck:EGFP-positive T lymphocytes (Fig. 1 E). foxp3a:EGFP-positive cells were also absent in rag1(lf) WKM (Fig. 1, D and F). Finally, we did not detect foxp3b expression in foxp3a:EGFP-positive cells or bulk lymphocytes, suggesting that foxp3b does not act redundantly with foxp3a, nor does it act in a separate, foxp3a-independent population of WKM lymphocytes (Fig. 1 G). Its expression in whole WKM could reflect its presence in WKM stromal cells or nonlymphocyte hematopoietic populations. The low level of foxp3b in WKM also does not exclude the possibility that it is expressed in lymphocytes present in other tissues. Collectively, based on light-scattering characteristics, cellular morphology, and requirement for rag1, foxp3a:EGFP cells are lymphocytes.

We looked for foxp3a:EGFP cells in other tissues. We assessed thymocytes from Tg(foxp3a:EGFP) animals using flow cytometry. In whole-thymus preparations, an EGFP-positive population that mapped to the lymphocyte gate was evident in Tg(foxp3a:EGFP) animals (Fig. 2 A). These foxp3a:EGFP-positive thymocytes had higher expression of the CD4 orthologue cd4-1 compared with bulk thymocytes, whereas expression of the CD8 orthologue cd8a (Fig. 2 B) was not detected. The regulator of T cell activation lck was expressed in the foxp3a:EGFP population. We also examined thymocytes from Tg(rag2:mCherry); Tg(foxp3a:EGFP) animals. Of explanted rag2:mCherry-positive thymocytes, rare cells were also foxp3a:EGFP-positive (Fig. 2 C). Because foxp3a-positive thymocytes were present in zebrafish, we speculate that development of some foxp3a-positive cells occurs in the thymus. In surveying other anatomical sites, we found that foxp3a:EGFP-positive cells were present in the skin. By flow cytometry, a population of EGFP-positive cells with characteristics of lymphocytes was evident (Fig. S1 B). Compared with bulk skin lymphocytes, these cells expressed higher levels of cd4-1, and cd8a expression was not detected (Fig. S1 C). We were unable to detect lck expression in this subpopulation. It is known that human skin-resident T reg cells express lower levels of LCK (Ali et al., 2017), although whether this lower expression reflects a functional difference in this subpopulation has not been investigated. We tracked the motility of foxp3:EGFP-positive cells in the skin with live-cell, time-lapse imaging (Fig. S1 D and Video 1). Directional motility of these cells was characterized by the formation of extensions at the leading edge of the cell followed by migration of the cell body. foxp3a:EGFP cells were highly motile in the skin, which may be evidence of their surveillance of this tissue. Finally, intermyseal connective tissue boundaries between muscle segments contained foxp3a:EGFP-positive cells (Fig. S1 E). In mammals, the thymus, skin, and muscle each contain FOXP3-expressing T reg cells (Panduro et al., 2016).

We characterized gene expression in foxp3a:EGFP-positive cells. For transcriptome analyses, EGFP-positive lymphocytes were harvested from the WKM of Tg(foxp3a:EGFP) and Tg(lck:EGFP) animals and used for low-input RNA extraction and library preparation followed by massively parallel sequencing (RNA-seq). Zebrafish genes with orthology to murine genes were used to compare expression profiles of foxp3a:EGFP-positive cells to mammalian T reg cells. Zebrafish orthologues were first compared with a narrowly defined T reg signature (Sugimoto et al., 2006). Many of these zebrafish orthologues were up-regulated in foxp3a:EGFP-positive cells compared with lck:EGFP-positive T cells (Fig. 3 A). Gene set enrichment analysis (GSEA) confirmed that there is a significant correlation between this T reg signature and zebrafish orthologues up-regulated in foxp3a:EGFP-positive cells (Fig. 3 B). GSEA also revealed that another, more broadly defined murine T reg signature (Hill et al., 2007) was enriched in foxp3a:EGFP-positive cells (Fig. 3 B). Relative expression levels of a subset of deregulated genes were confirmed using qRT-PCR (Fig. S2 A). These analyses indicate that foxp3a:EGFP cells have a gene expression signature that is consistent with expression signatures found in T reg cells. In light of these similarities, we are referring to foxp3a-positive cells in zebrafish as T reg–like cells.

We performed single-cell analyses to examine gene expression profiles and heterogeneity of foxp3a:EGFP-positive cells. foxp3a:EGFP-positive cells from WKM were isolated using FACS, and single cells were arrayed for qRT-PCR. qRT-PCR measured levels of 96 genes, which comprise mostly markers of different blood cell lineages as well as some housekeeping genes (Moore et al., 2016a). Using these expression data, unsupervised hierarchical clustering was performed to compare foxp3a:EGFP-positive cells (n = 35 cells) with unselected WKM cells from wild-type AB animals (n = 85 cells; Fig. S2 B). WKM cells fall into four distinct gene expression clusters, which correspond to T and B lymphocytes and erythroid and myeloid cells. Thirty foxp3a:EGFP-positive cells clustered with T lymphocytes. In considering individual genes (Fig. 3 C and Fig. S2 B), most cells (30/35) expressed T cell receptor genes. The largest fraction of cells (17/35) were cd4-1⁺cd8a^–. Two cells were cd4-1^–cd8a⁺. Several cells (11/35) were TCR-positive but cd4-1^–cd8a^–. Because extrathymic double-negative cells are rare and do not express FOXP3 (Hillhouse and Lesage, 2013), it is suspected that at least some of these cells express a second CD4 homologue, cd4-2, that is expressed in zebrafish lymphocytes (Yoon et al., 2015) but was not present in the single-cell analysis pipeline. cd4-2 was expressed in the foxp3a:EGFP-positive population at higher levels than in lck:EGFP-positive lymphocytes (Fig. S2 A). In addition to foxp3a:EGFP-positive cells that clustered with lymphocytes, there were five cells that segregated to the erythroid cluster. One of these cells expressed several red cell markers, and the four others were “low-expressor” cells that showed undetectable expression of most genes. Additional markers would be required to determine whether such cells represent distinct subpopulations of foxp3a:EGFP-positive cells. Overall, a comparison of gene expression to that in other T cell subsets indicates that most foxp3a:EGFP-positive cells are mature T cells with a bias toward CD4 positivity (Fig. 3 C). However, single-cell analyses reveal that there is some heterogeneity in the foxp3a:EGFP-positive population.

Mutations in mammalian FOXP3 genes cause loss of T reg cells, resulting in IPEX (immunodysregulation polyendocrinopathy enteropathy X-linked) syndrome in humans and the Scurfy phenotype in mice, both of which are characterized by an inflammatory response as a result of severe autoimmunity. To determine the consequences of foxp3a loss in zebrafish, we used site-directed transcription activator-like effector nucleases (TALENs) to induce foxp3a mutations. A region before the coding sequence for the foxp3a DNA-binding domain, in the sixth exon, was targeted. We injected TALEN mRNAs targeting this region into zebrafish embryos, allowed these animals to mature, then screened their progeny for foxp3a mutations. Animals with two independently derived mutations were recovered. The mutations, um252 and um253, created small deletions in the foxp3a gene and, by removing 11 and 14 nt, respectively, are predicted to cause premature truncations of the Foxp3a protein (Fig. 4 A). We performed qRT-PCR to measure RNA transcript levels in the WKM of foxp3a(um252) mutants and found that foxp3a transcripts were substantially reduced, presumably because of nonsense-mediated decay (Fig. 4 B).

We created a strain that combined the Tg(foxp3a:EGFP) transgene with the foxp3a(um252) mutation to determine whether T reg–like cells were affected by loss of foxp3a function. Flow cytometry of WKM from these animals showed, surprisingly, that there was an expansion of foxp3a:EGFP-positive cells (Fig. S3, A and B). To assess the nature of these cells, we isolated them using FACS and performed qRT-PCR (Fig. S3 C). As expected, the level of endogenous foxp3a transcripts in these cells was reduced. Normal foxp3a:EGFP-positive cells have low levels of cd8a; however, mutant foxp3a:EGFP-positive cells showed substantial expression of cd8a. Furthermore, expression of orthologues of mammalian T reg–enriched genes was deregulated in mutant foxp3a:EGFP-positive cells compared with their normal counterparts. foxp3b expression was not detected in mutant foxp3a:EGFP-positive cells, suggesting that it likely is not up-regulated to compensate for the loss of foxp3a. Based on these results, foxp3a:EGFP-positive cells from foxp3a mutants have an expression profile that is different from that of wild-type foxp3a:EGFP-positive cells. In some FOXP3 mutant IPEX patients, expanded populations of FOXP3-expressing cells have been observed, but these cells are similarly deranged in their misexpression of CD8 and other sublineage markers (Gavin et al., 2006). Furthermore, cells expressing a Foxp3:EGFP transgene are present in Foxp3 mutant mice, and these cells also misexpress CD8 and other T reg signature genes (Lin et al., 2007). As has been suggested in these mammalian cases, we speculate that presumptive T reg–like cells develop in foxp3a mutants but encounter a differentiation block in the absence of foxp3a function and ultimately adopt a hybrid profile with concurrent foxp3a promoter activity and cd8a expression.

The survival of foxp3a mutants was measured with cohorts of foxp3a heterozygous and homozygous mutant animals. After 20 wk, all control animals remained alive, whereas profound reductions in survival were observed in foxp3a(um252) and foxp3a(um253) homozygotes (Fig. 4 C). Reduced survival was also found in foxp3a(um252/um253) mutants, indicating that this survival defect was likely caused by loss of foxp3a function rather than background mutations present in the independently isolated foxp3a mutant strains. Aside from moderate wasting, no outwardly visible defects were present in foxp3a mutants before death.

foxp3a mutants were analyzed for hematopoietic defects. By flow cytometry, differences between control foxp3a heterozygous animals and foxp3a homozygous mutants were apparent (Fig. 5 A). Most noticeable was the increase of cells in the lymphoid-containing gate present in WKM from foxp3a(um252) and foxp3a(um253) mutants. Again, because similar hematopoietic defects were observed in foxp3a(um252/um253) mutants, this abnormality is likely caused by loss of foxp3a function. To better quantify differences in cell numbers, May-Grünwald Giemsa staining was performed on WKM, and nonerythroid cells were counted (Fig. 5 B). Lymphocytes were indeed increased, and a commensurate decrease in myeloid cells was observed. Inappropriate expansion of lymphocytes occurs in IPEX patients and Scurfy mice as a result of loss of T reg–mediated immune suppression. foxp3a mutants show a similar lymphocyte expansion.

We sacrificed and sectioned foxp3a mutants to assess whether additional internal defects could be compromising the health of these animals. Despite being smaller in overall size (Fig. S3 D), foxp3a mutants displayed marked splenomegaly (Fig. 5 C and Fig. S3 E). Visual examination of dissected spleens from foxp3(um252); Tg(lck:EGFP) animals suggested that this splenomegaly was characterized by an increased content of T lymphocytes (Fig. 5 D). Quantification of splenic lck:EGFP-positive cells from these animals confirmed this: although somewhat variable, spleens from foxp3a mutants contained a mean fourfold increase in T lymphocytes (Fig. 5 E and Fig. S3, F and G).

Pathological examination revealed profound inflammation of coelomic connective tissue surrounding the terminal colon and urogenital pore (Fig. 5 F). Inflammation extended into the ventral skin surrounding this pore (Fig. 5 G). In this and other inflamed regions, a lymphocyte-predominant mononuclear cell inflammatory infiltrate was present at much higher levels than in the same regions of control animals. In IPEX syndromic patients and Scurfy mice, persistent inflammation is present in multiple tissues and is driven by overexuberant cytokine production by and reactivity of lymphocytes. Inflammation in foxp3a mutant fish also likely results from aberrant lymphocytic reactivity, as foxp3a(um252); rag1(lf) fish have no inflammatory phenotype (Fig. 5, F and G). To assess whether an enhanced inflammatory response was associated with death of foxp3a mutants, we examined expression of inflammation marker genes in moribund foxp3a mutant and wild-type AB animals. In both colon and skin samples from foxp3a(um252) mutants, we observed an elevated level of the IFN-γ gene ifng1-2, as indicated by lower ΔCt values (Fig. 5 H). Levels of an interleukin-1β orthologue, il1b, were also elevated in colon samples of foxp3a(um252) mutants. il1b levels in skin were not significantly altered, nor were levels of a TNF-α homologue, tnfa, in skin or colon (Fig. S3 H). Based on these results, aspects of a pathological and molecular inflammatory response were associated with the death of foxp3a mutant animals.

In mammals, T reg cells play an essential role in suppressing autoreactive T cells. The adaptive immune system arose ∼500 million years ago in jawed fish (Flajnik and Kasahara, 2010), but it has not been clear whether mechanisms to curtail autoimmunity were present in the earliest adaptive immune systems. Arguing against the possibility of ancestral T reg cells, comparative analyses of FOXP3 orthologues showed that nonmammalian FOXP3 genes lack structural domains necessary for specification of the T reg sublineage (Andersen et al., 2012). However, our identification of T reg–like cells in zebrafish suggests that these cells were present at or near the emergence of adaptive immunity. In zebrafish, expression of cd4-1 is enriched in T reg–like cells, and these cells share a gene expression profile similar to that of mammalian T reg cells. Although an in vitro assessment of T reg suppression activity is currently not possible (appropriate cytokines have not yet been isolated), in vivo analyses of foxp3a mutants strongly suggest that zebrafish T reg–like cells have suppressor activity. The features of zebrafish T reg–like cells, along with the finding of foxp3-expressing blood cells in the pufferfish (Wen et al., 2011), suggest that T reg cells were an early sublineage of the adaptive immune system.

Because T reg–like cells are present in the zebrafish, the unique attributes of this model system can be harnessed to study this cell type. Ex vivo development of zebrafish embryos has facilitated studies of developmental hematopoiesis (Robertson et al., 2016), and through similar analyses, the developmental lineage and molecular specification of T reg–like cells can be examined. We have shown that direct imaging of T reg–like foxp3a:EGFP-positive cells can be performed. These cells are highly motile as they traffic through tissues. Through the combined use of the foxp3a:EGFP reporter and reporters of other cell types, it will be possible to visualize interactions of T reg–like cells with target T cells and other types of blood cells. In addition, implementation of zebrafish cancer models will enable the study of these cells in a native tumor setting. Such studies will serve as an important complement to T reg cell interaction studies that are currently performed through intravital imaging of xenografted cells in mammalian systems (Bauer et al., 2014).

Mammalian T reg cells can either develop in the thymus or differentiate without thymic education in the periphery. Development of thymic T reg cells is promoted by high-affinity interactions with self-peptide–MHC complexes, whereas peripheral T reg cells are thought to differentiate in response to non-self antigens, including allergens and commensal microbiota (Josefowicz et al., 2012). We found that a subset of thymocytes express foxp3a, suggesting that at least some T reg–like cells develop in the thymus and may be involved in suppressing autoimmunity. In human IPEX patients and Scurfy mice, autoimmunity caused by lack of T reg cells results in a gross inflammatory phenotype that is associated with excess lymphocyte activity (Lyon et al., 1990; Wildin et al., 2002). In addition, the expansion of T lymphocytes and splenomegaly are observed in FOXP3-deficient mammals (Powell et al., 1982; Kim et al., 2007). Altogether, foxp3a mutant zebrafish have a similar constellation of defects as, and can thus serve as a suitable model for, FOXP3-deficient mammals.

Although the existence of T reg cells was established nearly two decades ago, several important issues regarding their functions remain to be elucidated. In particular, the mechanisms by which they dominantly suppress immune responses in vivo, the sites at which they interact with target cells, and their involvement in pathological conditions such as tumor immune suppression are all poorly understood. Through the identification of zebrafish T reg–like cells and the development of foxp3a-deficient mutants, many of these issues can be addressed.

Zebrafish were maintained at 28.5°C with a 14:10-h day:night cycle. Tg(lck:EGFP) (Langenau et al., 2004), rag1(hu1999) (Wienholds et al., 2002) strains were used. rag1(hu1999) mutants were provided by L. Petrie-Hanson (Mississippi State University, Starkville, MS). AB was used as the wild-type strain. All zebrafish studies were performed according to Institutional Animal Care and Use Committee guidelines and University of Massachusetts Medical School Animal Care Committee protocols.

A gene structure for foxp3a was determined based on available cDNA sequences. A 2,446-bp promoter fragment upstream of the 5′-most noncoding exon was PCR-amplified (corresponding to −7,448 to −5,003 relative to the foxp3a initiator ATG) and recombined into the Gateway P4P1r vector (Thermo Fisher Scientific). Using multisite Gateway recombination, the foxp3a promoter was placed upstream of the EGFP gene followed by an SV40 polyadenylation signal in the pDestTol2CG2 vector (Kwan et al., 2007). Single-cell embryos were injected and transgenic founders selected based on heart fluorescence.

Total RNA was isolated from WKM and sorted cells using TRIzol-chloroform (Thermo Fisher Scientific) extraction followed by RNA cleanup (QIAGEN RNeasy). RNA was amplified and reverse-transcribed using the Ovation Pico WTA System v2 (NuGen). qRT-PCR was performed with SYBR green master mix (Applied Biosystems) using the following primers: foxp3a: forward, 5′-TGCTTTGTGCGTGTTGAAGG-3′, and reverse 5′-ACAGCTGAAATGGTGCCATC-3′; foxp3b: forward, 5′-TTTAAAGCAGCACGCTTCAG-3′, and reverse, 5′-AGAAACTGCATTCGCTGCTC-3′; cd4-1: forward, 5′-TCAACACCAAGACCATCAGC-3′, and reverse, 5′-GCACATGTCCATTTCACCTC-3′; cd4-2: forward, 5′-TTATTCCTGCGCACAGTCAG-3′, and reverse, 5′-TTGGGCCTGGTTTAGCAAAC-3′; cd8a: forward, 5′-TCGGAGGTTGTGGACTTTTC-3′, and reverse, 5′-TTGTAATGGTGGGGACATCG-3′; ctla4: forward, 5′-GGGAACGGCACTGTTGTTTAC-3′, and reverse, 5′-TGTCTGGCTCTTGCTTTGAC-3′; tnfrsf18: forward, 5′-AACGATGTTGAGGACAGCAC-3′, and reverse, 5′-TGCATTGTGGGTCTGGTTTC-3′; ikzf4: forward, 5′-ATTGCAATGGCCGTTCGTAC-3′ and reverse, 5′-ATGGAGTTAGCACTGAGTGAGC-3′; irf4b: forward, 5′-AACTGCAAACTGCTGGACAC-3′, and reverse, 5′-ACAGCAACACTTGGGAACAC-3′; and lck: forward, 5′-AGTGGCCCAGCATTGATATCTG-3′, and reverse, 5′-ATCAAAAGAGCCGCAGTTCC-3′. Data from three technical replicates per sample were obtained. In single-cell analyses, cells were arrayed, RNA was extracted, and qRT-PCR was performed using the Fluidigm BioMark HD as described (Moore et al., 2016a,b). The panel of genes profiled was described previously (Moore et al., 2016a). Ct values were recovered from the BioMark HD. The quality threshold was set to 0.65, and a linear derivative as baseline correction was used. Using SINGuLAR software (Fluidigm), the limit of detection was set to a Ct of 28; expression was calculated using default settings. Sample wells that failed to express housekeeping genes (eef1a1l1 and actb1) or other lineage-specific markers were deemed to lack a cell and were eliminated from further analysis. Unsupervised hierarchical clustering was performed using SINGuLAR. foxp3a:EGFP-positive single-cell samples were combined and compared with the datasets generated previously (Moore et al., 2016a), and lineage assignments were made after hierarchical clustering. The hierarchical clustering and violin plots were also generated using the SINGuLAR R package, with the limit of detection set to a Ct of 28.

EGFP-positive lymphocytes were sorted from the WKM of Tg(foxp3a:EGFP) and Tg(lck:EGFP) animals to prepare three biological replicate samples of each genotype. Sorted cells were used to prepare amplified double-stranded (ds)-cDNA using Clontech SMART-Seq v4 ultra low input RNA kit, according to the manufacturer’s recommendations. Amplified ds-cDNA was purified, fragmented, and used to prepare libraries using TruSeq RNA sample preparation kit v2 (Illumina). After validation using a DNA 1000 LabChip (Agilent), samples were submitted for RNA-seq. As described previously (Anelli et al., 2017), paired end reads were aligned using Star v2.3 (Dobin et al., 2013) to the zebrafish genome (GRCz10; Howe et al., 2013) using the Ensembl transcriptome. Analysis of differential gene expression was performed using DESeq2 (Love et al., 2014). Orthology to murine genes was determined using Ensembl.

For GSEA (Subramanian et al., 2005), a rank-ordered gene list was made with normalized FPKM values from foxp3a:EGFP-expressing cells compared with lck:EGFP-expressing cells. Default parameters of GSEA were used, and Student’s t test was used to calculate significance.

TALENs targeting the foxp3a locus were designed using the TAL Effector Nucleotide Targeter 2.0. TALEs targeting the sixth coding exon of the foxp3a gene (binding sites in caps, spacer in lowercase: 5′-CCTGGAGACAGAAGCATAggccagctgaggatgcAGAAAGACAGAGTACAGC-3′) were constructed using the GoldenGate TALEN and TAL Effector kit 2.0 (Addgene). RNA encoding foxp3a-directed TALENs was injected into single-cell AB embryos, and progeny of injected animals were screened for mutations. F1 progeny containing foxp3a mutations were outcrossed to AB animals, and strains were established. Mutant foxp3a alleles were detected by PCR-amplifying an amplicon containing the TALEN target site (forward primer, 5′-TGATTGGGTTGTGGACACC-3′; reverse primer, 5′-CAGTTCTGAAGGCAAAGGG-3′) and digesting with PvuII. Wild-type alleles are cut with PvuII, whereas mutant alleles lack the PvuII site.

WKMs were harvested and prepared as previously described (Traver et al., 2003), followed by flow cytometry using a FACSAriaII instrument (BD Biosciences). Spleens were harvested, triturated, and subjected to the same preparation and flow cytometry analysis as WKMs. FACS of WKMs was performed using a FACSVantage instrument (BD Biosciences). Flow cytometry data were analyzed using FlowJo software. Single cells were sorted into 96-well plates using a FACSAria Fusion Cell Sorter (BD Biosciences).

Cytospins were performed by applying 150 µl WKM cell suspension to funnels and centrifuging for 5 min at 800 g. Slides with cells were dried for 10–15 min. Stainings with May-Grünwald and Giemsa solutions were performed as described (Traver et al., 2003). For histopathology, humanely euthanized fish were fixed whole in 10% neutral-buffered formalin, sectioned sagittally, and routinely processed into paraffin blocks. Hematoxylin and eosin–stained sections were reviewed by A.B. Rogers, who is a board-certified veterinary pathologist. Imaging was performed using a DM5500 microscope, and LAS software was used for image acquisition (Leica).

Time-lapse videos were generated from tumor-bearing zebrafish by first anesthetizing fish in 0.17 mg/ml tricaine. Anesthetized fish were placed on nontoxic modeling clay, and Tygon tubing was placed over the fish between the caudal fin and the body to minimize movement. The anesthetized animal was fully submerged in a water bath containing 0.17 mg/ml tricaine, and a p200 pipette tip was placed at the edge of the animal’s mouth to provide a steady stream of 0.17 mg/ml tricaine in tank water through the mouth and over the gills of the fish. This was accomplished by connecting a peristaltic pump to a pipette tip, which delivered oxygenated water with 0.17 mg/ml tricaine at a rate of 750 µl/min. Cells were imaged with a 20× water-dipping lens on a DM5500 microscope (Leica). Once the imaging experiment was completed, the fish was revived by pumping fresh water through the pipette tip into the mouth until motion in the gills was evident. After intubation, the Tygon tubing and pipette tip were removed, and the animal was placed back into fresh water for recovery. Image acquisition, background subtraction, and deconvolution were performed with Leica LAS software.

Statistical analyses were performed using Prism software by GraphPad. Unless otherwise specified, comparisons between sample groups were performed using unpaired t tests. Comparisons between survival curves were performed by log-rank (Mantel–Cox) tests. Statistical analyses of single-cell profiling data were performed as previously described (Moore et al., 2016a).

Fig. S1 shows flow cytometry profiles of foxp3a:EGFP-positive cells in WKM and skin, quantification of gene expression in foxp3a:EGFP-positive cells from skin, and imaging of foxp3a:EGFP-positive cells in skin and muscle. Fig. S2 shows quantification of gene expression in and single-cell profiling of foxp3a:EGFP-positive lymphocytes from WKM. Fig. S3 shows characterization of foxp3a mutants, including expansion of Tg(foxp3a:EGFP)-expressing lymphocytes, quantification of gene expression, small stature, expansion of splenic size and splenic lymphocytes, and quantification of tnfa in colon and skin lymphocytes. Video 1 shows live imaging of motile foxp3a:EGFP-positive cells in skin.

We thank Leslie Berg for helpful discussions; Patrick White, Ed Jaskolski, and the UMMS Animal Medicine Department for fish care; and UMMS flow cytometry core for sorting experiments. We acknowledge the Epigenomics Core Facility at Weill Cornell Medical College for library preparation and next-generation sequencing services.

C.A. Painter was supported by the Cancer Research Institute Irvington Postdoctoral Fellowship Program. Research was supported by Kimmel Scholar Award SKF-13-123 (C.J. Ceol), NIH R24 OD016761 (D.M. Langenau), and NIH R01 AR063850 (C.J. Ceol). The content is solely the responsibility of the authors and does not necessarily represent the official views of the Department of Defense or National Institutes of Health.

The authors declare no competing financial interests.

Author contributions: M. Kasheta, C.A. Painter, F.E. Moore, A. Bryll, E. Freiman, and C.J. Ceol performed experiments; M. Kasheta, C.A. Painter, F.E. Moore, R. Lobbardi, D. Stachura, A.B. Rogers, Y. Houvras, D.M. Langenau, and C.J. Ceol designed experiments and analyzed data; and C.J. Ceol wrote the manuscript with help from M. Kasheta, F.E. Moore, R. Lobbardi, A.B. Rogers, and Y. Houvras.