Integrated molecular analysis identifies a conserved pericyte gene signature in zebrafish

The vertebrate circulatory system comprises blood vessels and a heart, which are lined by a single layer of endothelial cells. The circulatory system delivers oxygen to tissues, while disposing of carbon dioxide and other waste products. Blood vessels also provide a conduit for immune cells and hormones of the endocrine system. Delivery and exchange of blood constituents to tissue is facilitated by capillary beds that are functionally tuned to their anatomical location (Holm et al., 2018). For example, in the central nervous system, capillaries of the blood–brain barrier exhibit highly selective transport activity to limit entry of circulating factors that could damage neurons (Zhao et al., 2015). By contrast, liver and kidney capillary beds are important for filtering and removing serum components from circulation (Holm et al., 2018), requiring less restricted passage through the endothelial lining. In capillaries, vascular mural cells, known as pericytes, directly contact and share a basement membrane with endothelial cells and regulate their transport function (Holm et al., 2018). In the absence of pericytes, capillary function is altered. A notable example is the blood–brain barrier, where reduced pericyte coverage causes vascular leakage due to increased endothelial transcytosis (Armulik et al., 2010). Thus, pericytes serve as gatekeepers of capillary function within the circulatory system.

The mechanisms associated with pericyte-mediated organotypic capillary function have been well studied (Armulik et al., 2011; Holm et al., 2018). However, much less is known about embryonic pericyte development. Early studies in avian embryos revealed a dual origin for pericytes: in the rostral regions of the central nervous system they derive from neural crest, whereas elsewhere they arise from mesoderm (Etchevers et al., 2001; Korn et al., 2002). Studies in mouse and zebrafish suggest that this is a conserved aspect of pericyte development (Ando et al., 2019; Simon et al., 2012). A number of signaling effectors have also been implicated in pericyte development, including Platelet derived growth factor receptor β (Pdgfrβ), which is highly expressed on pericytes, and its endothelial-expressed ligand (Pdgfβ), along with components of the Notch pathway (Ando et al., 2021; Dieguez-Hurtado et al., 2019; Hellstrom et al., 1999; Kofler et al., 2015; Levéen et al., 1994; Lindahl et al., 1997; Wang et al., 2014). However, little is known about acquisition of pericyte identity during embryogenesis. This knowledge gap stems, in part, from the lack of definitive pericyte marker genes in developmental contexts.

Numerous studies have focused on identifying pericyte-specific genes. Early efforts used microarray or RNA-seq performed on microvascular fragments from mutant backgrounds that reduce the numbers of pericytes or incorporated transgenic lines to isolate vascular mural cell populations (Bondjers et al., 2006, 2003; He et al., 2016; Jung et al., 2018). More recently, single cell RNA-sequencing (scRNA-seq) has begun to reveal new markers for pericytes (Vanlandewijck et al., 2018). However, these efforts have focused on adult mouse tissues and studies at embryonic stages, or in other vertebrate models, are lacking. In the zebrafish embryo, which has proven useful for investigating pericyte development, only a few pericyte-specific genes have been identified, limiting the utility of this model (Ando et al., 2016, 2019, 2021; Chauhan et al., 2016; Wang et al., 2014). To address this issue, we applied scRNA-seq on pdgfrb-positive cells from zebrafish larvae, with parallel integration of related datasets, to better define a developmental pericyte gene signature. We subsequently leveraged our findings to develop new reporter lines to directly visualize pericytes in vivo.

Transgenic zebrafish lines made using the pdgfrb locus drive fluorescent protein expression in pericytes and numerous non-vascular cell types (Ando et al., 2016; Vanhollebeke et al., 2015), making definitive identification of pericyte markers challenging with bulk RNA-seq. Therefore, to identify pericyte-expressed genes, we performed 3′ scRNA-seq using Egfp-positive cells isolated by fluorescent activated cell sorting (FACS) from 5 days postfertilization (dpf) larvae bearing TgBAC(pdgfrb:egfp)^ncv22 (Ando et al., 2016). After filtering and integrated analysis of replicate libraries, we identified more than 50 distinct clusters from 12,865 cells (Fig. S1A,B), underscoring the heterogeneity of pdgfrb:egfp-positive cells. We next used selected cluster-specific genes to assign cell identities (Fig. S1C; Table S1.1, S1.2; see Materials and Methods). For example, actinodin1 (and1) defined epidermal cells, whereas mitfa identified melanocytes and neurons expressed elavl3 (Fig. S1C; Table S1.1, S1.2; Kim et al., 1996; Lawson et al., 2020; Lister et al., 1999; Zhang et al., 2010). We could not assign identities to ten clusters, although cells in these cases exhibited high expression of collagens (col1a2 and col12a1a; Fig. S1C), suggesting fibroblast characteristics (Rajan et al., 2020). Surprisingly, most clusters displayed very low levels of endogenous pdgfrb transcript (Fig. S1C; Table S1.3), suggesting that these cells may have been isolated due to ectopic pdgfrb:egfp expression. Alternatively, pdgfrb:egfp expression in progenitor cells and Egfp perdurance may have contributed to their capture by FACS. Given high pdgfrb:egfp levels in pericytes (Ando et al., 2016), we focused on clusters expressing higher relative levels of endogenous pdgfrb (Fig. 1A, pdgfrb^hi cells, log₂ average expression >0.6; Table S1.3). We also included smooth muscle cell clusters independent of pdgfrb levels, as their gene expression profiles can be similar to pericytes (Vanlandewijck et al., 2018).

Among pdgfrb^hi cells, cluster 39 expressed the highest pdgfrb levels (Table S1.3). Moreover, mouse orthologs for half of the ten most enriched genes in this cluster (ranked by log₂ fold change, adjusted P<0.05, Wilcoxon Rank Sum test; Table S1.4) are expressed in pericytes, including duplicates of regulator of G protein signaling 5 family genes (rgs5a and rgs5b), endothelin receptor type Aa (ednraa), and nidogen 1b (nid1b) (Fig. 1B; Table S1.4; Bondjers et al., 2003; Kitazawa et al., 2011; Sakhneny et al., 2021). Other significantly enriched genes (log₂ fold change>0.75, adjusted P<0.05) included notch3, which is expressed in zebrafish brain pericytes (Wang et al., 2014) and potassium inwardly rectifying channel subfamily J member 8 (kcnj8), a marker of mouse brain pericytes (Bondjers et al., 2006) Fig. 1B, Table S1.4). We similarly noted cluster-enriched expression of endosialin/CD248 molecule a (cd248a), forkhead box C1b (foxc1b), melanoma cell adhesion molecule b (mcamb) and laminin beta 1b (lamb1b), orthologs of which are expressed in mouse pericytes (Fig. 1B, Table S1.4; Bagley et al., 2008; Middleton et al., 2005; Sakhneny et al., 2021; Siegenthaler et al., 2013). Consistent with pericyte identity, cells in cluster 39 did not express smooth muscle actin alpha 2 (acta2), a definitive marker of zebrafish smooth muscle cells (SMCs) (Georgijevic et al., 2007; Whitesell et al., 2014), or other SMC genes such as desmin b (desmb), myocardin (myocd), calponin 1b (cnn1b) or transgelin (tagln) (Fig. 1C; Georgijevic et al., 2007; Wang et al., 2014; Whitesell et al., 2019). Thus, cluster 39 likely comprises pericytes and we hereafter refer to genes enriched in these cells as ‘pericyte genes’.

To further characterize the expression of pericyte genes, we incorporated bulk RNA-seq analysis of TgBAC(pdgfrb:citrine)^s1010 cells (Lawson et al., 2020). We considered two differential gene sets in which pericyte-specific genes would be expected to be found. First, genes enriched in wild-type pdgfrb:citrine-positive cells compared with -negative cells (log₂ fold change>1, adjusted P<0.05; Table S2.1). Second, genes reduced in pdgfrb:citrine-positive cells from pdgfrb^um148 mutants, which lack pericytes (Ando et al., 2021), compared with wild-type siblings at 5 dpf (log₂ fold change<−1, adjusted P<0.05; Table S2.2). Intersection of these datasets revealed 23 out of 105 pericyte genes with enrichment in pdgfrb:citrine cells and reduction in pdgfrb^um148 mutants (Fig. 1D; Table S3.1, S3.2), including interleukin 11a (il11a), NDUFA4 mitochondrial complex associated like 2a (ndufa4l2a), and potassium voltage-gated channel, Isk-related family, member 4 (kcne4), which have not been previously characterized as pericyte-expressed. Most pericyte genes (90 out of 105), including many previously identified orthologs, exhibited enrichment in pdgfrb:citrine-positive compared with -negative cells, but were not reduced in pdgfrb^um148 mutants (Fig. 1D; Table S3.1). This is likely owing to their widespread expression in other pdgfrb-positive cell types not affected in pdgfrb^um148 mutant larvae. In this regard, we note that pericytes make up less than 1% of pdgfrb:egfp-positive cells based on scRNA-seq (Fig. S1A; Table S1.2). Despite its broad expression, we detected pdgfrb as downregulated, consistent with our previous observation that the um148 allele causes non-sense mediated decay of the mutant transcript (Kok et al., 2015).

An additional consideration we made for zebrafish pericyte genes was whether they were conserved among vertebrate models. Therefore, we incorporated scRNA-seq analysis of mouse pericytes (He et al., 2018; Vanlandewijck et al., 2018). Re-clustering datasets from brain and lung vascular cells identified pericyte populations characterized by high levels of Cspg4, Pdgfrb, Notch3, Mcam and Ifitm1 and low expression of SMC markers (Myocd and Cnn1) (Fig. S2A-D; Table S4). Intersection of mouse and zebrafish pericyte genes revealed a large degree of overlap (39/93 zebrafish genes found in mouse pericytes; intersection based on mouse orthologs, see Materials and Methods; Fig. 1E; Table S3.1). Notably, overlap between zebrafish and mouse pericytes from either tissue was comparable with that between mouse pericytes from brain and lung. This analysis identified a core set of 14 pericyte genes that were commonly enriched in pericytes across the three datasets, including Notch3 and Pdgfrb (Fig. 1E; Table S3.3). From this core gene set, all displayed significant enrichment in pdgfrb^pos cells by bulk RNA-seq, but only two (pdgfrb and ndufa4l2a) showed reduction in pdgfrb^um148 mutants (Table S3.3), again suggesting that most of these genes were not restricted to pericytes. This is consistent with a high proportion of non-pericyte cells expressing many of the pericyte genes from our scRNA-seq analysis (pct.2 value in Table S1.4). Nonetheless, these observations demonstrate that zebrafish larval pericytes share a conserved molecular signature with those in adult mouse tissues.

Our molecular analysis establishes a gene signature for zebrafish pericytes at the larval stage. To demonstrate the utility of our datasets for studying zebrafish pericytes, we chose two previously uncharacterized candidate pericyte genes for generation of transgenic reporter lines. The first was ndufa4l2a, which is among the most consistently enriched and specific genes from our analysis (Table S3.1). To assess ndufa4l2a expression, we constructed a recombinant bacterial artificial chromosome (BAC) with super folder green fluorescent protein (sfGFP) inserted into the first exon and used this to generate a stable transgenic line [TgBAC(ndufa4l2a:sfgfp)^um382]. Confocal imaging of the brain vasculature in TgBAC(ndufa4l2a:sfgfp)^um382 larvae at 5 dpf revealed expression in cells closely associated with blood vessels throughout the mid- and hindbrain (Fig. 2A). Cells expressing ndufa4l2a:sfgfp also co-expressed abcc9:gal4ff;uas:rfp (referred to hereafter as abcc9:rfp), a known pericyte marker (Ando et al., 2019; Vanlandewijck et al., 2018), and appeared to wrap around brain blood vessels, consistent with pericyte morphology (Fig. 2A,B). We also observed expression of ndufa4l2a:sfgfp in abcc9:rfp-positive pericytes lining retinal endothelial cells (Fig. 2C). In the trunk, ndufa4l2a:sfgfp-positive cells expressed abcc9:rfp and were closely associated with intersegmental vessels at 5 dpf (Fig. 2D). We did not detect sfGFP expression in vascular mural cells along the dorsal aorta (Fig. 2D), most of which express acta2 at this stage consistent with their identity as vascular SMCs (VSMCs) (Whitesell et al., 2014). In the trunk and cranial vasculature, we noted ndufa4l2a:sfgfp-positive cells that did not express abcc9:rfp (Fig. 2A,D, denoted by asterisks), likely owing to mosaic silencing that occurs with the Gal4/UAS system (Goll et al., 2009). Thus, the ndufa4l2a:sfgfp appeared largely restricted to abcc9-positive pericytes throughout the vasculature. We also noted non-vascular expression throughout the epidermis, as well as cells with a mesenchymal appearance on the surface of the lower jaw (Fig. S3A,B). However, this non-vascular expression was limited to these particular locations and overall the ndufa4l2a:sfgfp appeared to be more pericyte-restricted than pdgfrb-driven transgenes. Interestingly, none of the epidermis clusters from scRNA-seq analysis expressed ndufa4l2a (Fig. S1C), suggesting that the observed ndufa4l2a:sfgfp expression may be ectopic.

We also generated a reporter line using a BAC encompassing the kcne4 locus [TgBAC(kcne4:sfgfp)^um333]. Similar to TgBAC(ndufa4l2a:sfgfp)^um382, kcne4:sfgfp was expressed in brain pericytes, in this case co-expressed with a pdgfrb:gal4ff transgene driving uas:rfp (Fig. 3A). We also observed kcne4:sfgfp expression in pdgfrb-positive retinal pericytes (Fig. 3B) and along the intersegmental vessels (ISVs) and dorsal aorta in the trunk (Fig. S3C). In most cases, kcne4:sfgfp appeared to be at much lower expression levels than ndufa4l2a:sfgfp, potentially reducing the utility of the former in routine imaging of pericytes. Unlike endogenous ndufa4l2a, which was seen to be restricted to pericytes using scRNA-seq analysis, kcne4 was also enriched in SMCs (Fig. 1B; Table S1.1), consistent with previous zebrafish studies (Whitesell et al., 2019). However, we did not observe kcne4:sfgfp expression in acta2:mcherry-positive VSMCs on the dorsal aorta or the Circle of Willis in the brain (Fig. 3C,D). In both cases kcne4:sfgfp could be detected in pericytes on nearby blood vessels (Fig. 3C,D). In the intestine, which is surrounded by SMCs, we observed scattered kcne4:sfgfp-positive cells that co-expressed acta2:mcherry, along with acta2-negative cells in direct contact with blood vessels, consistent with pericyte morphology (Fig. 3E). In contrast to the dorsal aorta, kcne4:sfgfp was expressed in acta2-positive VSMCs in aortic arch vessels branching from the ventral aorta (VA) and in the cardiac ventricle, but not along the VA itself (Fig. 3F,G). sfGFP-positive cells that appear on the VA in Fig. 3F,G are located ventrally and are not directly associated with blood vessels (Fig. S3D,E). We also detected acta2-negative cells with pericyte morphology that expressed kcne4:sfgfp along the hypobranchial artery, but not prominently on other vessels in this region (Fig. 3F). Thus, kcne4:sfgfp is expressed in pericytes, as well as anatomically restricted subsets of SMCs in the gut and VSMCs in the head. As with ndufa4l2a:sfgfp, we also detected epidermal expression in TgBAC(kcne4:sfgfp)^um333 larvae (Fig. S3C). Taken together, these transgenic lines underscore the utility of our molecular analysis and provide new tools for studying pericytes and VSMCs in zebrafish.

The zebrafish has been a valuable model for studying vascular development, providing novel insights into conserved processes that contribute to blood vessel growth and function. However, efforts have focused on endothelial cells, facilitated by a wide range of technical resources for their study in the zebrafish embryo, including numerous transgenics, mutants and a broad catalog of molecular markers. By contrast, only a handful of studies have addressed the development of zebrafish pericytes, and genetic tools in this regard are limited. Our current study begins to fill these gaps by establishing a molecular signature comprising more than 100 genes expressed in zebrafish larval pericytes. We demonstrate the utility of this dataset by establishing reporter lines for visualizing pericytes in zebrafish larvae. Together, our findings will provide a foundation for identifying functional candidate genes and further developing tools to dissect the genetic requirements for pericyte development during embryogenesis.

All studies were performed under the auspices of animal protocols approved by the University of Massachusetts Medical School Institutional Animal Care and Use Committee. TgBAC(pdgfrb:citrine)^s1010, TgBAC(pdgfrb:egfp)^ncv22, pdgfrb^um148, TgBAC(pdgfrb:gal4ff)^ncv24, TgBAC(abcc9:gal4ff)^ncv34­, Tg(uas:rfp)^nkuasrfp1a and Tg(acta2:mcherry)^ca8 have been described elsewhere (Ando et al., 2016, 2019; Asakawa et al., 2008; Kok et al., 2015; Vanhollebeke et al., 2015; Whitesell et al., 2014). Lines generated in this study are described below.

For bulk RNA-seq, citrine-positive cells from approximately 25 TgBAC(pdgfrb:citrine)^s1010;pdgfrb^um148 mutant larvae at 5 dpf were obtained in the course of previous studies in parallel to cells from wild-type larvae (Lawson et al., 2020; Whitesell et al., 2014). Mutant siblings were identified and separated based on the absence of brain pericytes by confocal microscopy. Three separate groups of mutants were used to make replicate libraries. Larval dissociation, FACS isolation (University of Massachusetts Medical School Flow Cytometry Core), RNA isolation and library construction were applied as described previously (Lawson et al., 2020; Quillien et al., 2017; Whitesell et al., 2014). For scRNA-seq, ∼100 wild-type TgBAC(pdgfrb:egfp)^ncv22 larvae at 5 dpf were dissociated as described previously (Lawson et al., 2020; Quillien et al., 2017; Whitesell et al., 2014), except cells were not fixed. Live cells were subjected to FACS to isolate Egfp-positive cells and ∼10,000 cells were loaded onto a 10x Chromium (10x Genomics) for generation of single cell droplets. We constructed 3′ 10x scRNA-seq libraries according to the manufacturer's recommendations (v3, 10x Genomics) and we performed pilot depth sequencing on a MiSeq (Illumina; University of Massachusetts Medical School Deep Sequencing Core Labs) to verify quality and captured cell numbers. Subsequently, libraries were sequenced to greater depth on a HiSeq2500 (Illumina; Genewiz).

All pipelines were run remotely at the Massachusetts Green High Performance Computing Center (MGHPCC). We performed mapping of bulk RNA-seq using a DolphinNext pipeline as described previously (Lawson et al., 2020). Bulk RNA-seq datasets from wild-type pdgfrb:citrine-positive and -negative cells have been previously published and are available in GEO (GSE152759). Reads were mapped onto GRCz11 and expression levels quantified using our custom transcript annotation (V4.3.2; Lawson et al., 2020). We identified differentially expressed genes using DeSeq2 run in DEbrowser, as described previously (Lawson et al., 2020). Reads from scRNA-seq libraries were mapped onto GRCz11 using Cell Ranger (v3.1.0, 10x Genomics) in a DolphinNext environment and gene quantifications were made using the V4.3.2 annotation as described previously (Lawson et al., 2020). We used Seurat (v4; Hao et al., 2021) for clustering using the output from Cell Ranger. We used Rstudio in the Open OnDemand environment at MGHPCC to run Seurat as previously described (Lawson et al., 2020). We applied integrated analysis to combine two replicate scRNA-seq libraries from pdgfrb:egfp cells (R script for pdgfrb:egfp analysis available in Supplementary file1). Cells were filtered for feature counts and mitochondrial proportion (see Supplementary file 1), yielding a total of 12,865 cells across the two libraries. Following integration and identification of common anchors, Uniform Manifold Approximation and Projection (UMAP) was applied to reduce dimensionality based on the first 40 principal components and clusters were identified using a resolution 2. We identified all cluster-specific markers and used these to assess cluster identity by manual comparison to whole-mount in situ hybridization patterns available through the Zebrafish Information Network (ZFIN; www.zfin.org). To identify putative pericyte cell clusters, we performed text-based searching of all cluster-specific genes for selected genes previously identified as pericyte-expressed. To integrate bulk RNA-seq with the pericyte gene sets from scRNA-seq analysis, we used the top 105 pericyte-enriched genes (log₂FC>0.75 relative to all other clusters, adjusted P<0.05), along with those enriched in pdgfrb:citrine-positive versus -negative cells (log₂FC>1, adjusted P<0.05) and reduced in pdgfrb:citrine pdgfb^um148 cells compared with wild-type cells (log₂FC<−1, adjusted P<0.05). To obtain mouse pericyte gene sets, we used published data downloaded from GEO as raw count tables (GSE99235, GSE98816; He et al., 2018; Vanlandewijck et al., 2018). Each dataset was individually analyzed using Seurat run as above (see Supplementary file 2 for commands and parameters for clustering). Pericyte clusters were identified based on expression of marker genes identified previously (Vanlandewijck et al., 2018). Pericyte gene sets were generated from the top enriched genes (log₂ fold change>0.75 relative to all other clusters, adjusted P<0.05). To identify mouse orthologs for zebrafish pericyte genes, we used two data sources: a curated HomoloGene dataset from the Mouse Genome Informatics database (http://www.informatics.jax.org/downloads/reports/HOM_AllOrganism.rpt) and zebrafish:mouse orthology from ZFIN (https://zfin.org/downloads/mouse_orthos.txt). From 105 zebrafish pericyte genes, we identified 93 with mouse orthologs.

Two mouse orthologs, Lamb1 and Rbpms2, are orthologous to pericyte genes that are duplicated in zebrafish. Gene sets were intersected to generate Venn diagrams based on gene symbol using intervene (Khan and Mathelier, 2017) with the following commands: --type list --save-overlaps. Proportional Venn diagrams were generated with resulting data using eulerr (https://github.com/jolars/eulerr).

BAC plasmids containing kcne4 (DKEY-16C17) or ndufa4l2a (DKEY-11F4) loci were purchased from Source BioScience and verified for presence of the target sequence by PCR for the first coding exon (exon 2 for kcne4, exon 1 for ndufa4l2a; primers kcne4-e1-F and kcne4-e1-R, or ndufa4l2a-e1-F and ndufa4l2a-e1-R, see Table S5 for primer sequences). Recombinant BACs were constructed based on established protocols using pRedET (GeneBridges, Bussmann, 2011 #145). An Isce-Tol2-Tn5neo cassette was amplified from pGEM-Isce-Tol2-Tn5-neo (Addgene plasmid #176619) by PCR using pBelo_pGEMT-F and pBelo_pGEMT-R (Table S5), gel purified and electroporated into bacteria carrying the target BAC and pRedET. Correctly targeted clones were identified by PCR. We next introduced the sfGFP coding sequence with an optimized Kozak consensus sequence into the first coding exon of each locus. A template plasmid (p-sfGFP-FRT-amp-FRT, Addgene plasmid #176620) was constructed by inserting a gBlock encoding sfGFP followed by an SV40 polyA sequence and an FRT-flanked ampicillin cassette in place of Egfp in the pEgfp-C1 plasmid (Clontech). For kcne4 and ndufa4l2a BACs, we amplified the sfGFP/amp cassette using primers kcne4_sfGFP_f and kcne4_FRT-Amp_R, or ndufa4l2a_sfGFP_F1 and ndufa4l2a_FRT-Amp_R1 (Table S5), respectively. Templates were gel-purified and electroporated into bacteria bearing the target BAC and pRedET. Following recombination, correct clones were identified by PCR and the ampicillin resistance cassette was removed using the 707-FLPe plasmid according to recommended instructions (GeneBridges). Recombinant BACs were purified using PureLink™ HiPure Plasmid Midiprep Kit (Thermo Fisher Scientific) and injected into one-cell stage zebrafish embryos together with I-SceI (New England Biolabs). Only embryos exhibiting low mosaic sfGFP expression were grown to adulthood. We identified a single founder for each recombinant BAC. These lines are referred to as: TgBAC(kcne4:sfgfp)^um333 and TgBAC(ndufa4l2a:sfgfp)^um382. Imaging of larval stages was performed by confocal microscopy as previously described (Ando et al., 2021).

We thank members of the Lawson Lab for helpful comments on the manuscript. We thank John Polli and Patrick White for their efforts in fish care and facility maintenance. We thank Naoki Mochizuki and Koji Ando for providing transgenic lines used in this study.

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