Gene Expression Analysis of Zebrafish Melanocytes, Iridophores, and Retinal Pigmented Epithelium Reveals Indicators of Biological Function and Developmental Origin

In order to facilitate understanding of pigment cell biology, we developed a method to concomitantly purify melanocytes, iridophores, and retinal pigmented epithelium from zebrafish, and analyzed their transcriptomes. Comparing expression data from these cell types and whole embryos allowed us to reveal gene expression co-enrichment in melanocytes and retinal pigmented epithelium, as well as in melanocytes and iridophores. We found 214 genes co-enriched in melanocytes and retinal pigmented epithelium, indicating the shared functions of melanin-producing cells. We found 62 genes significantly co-enriched in melanocytes and iridophores, illustrative of their shared developmental origins from the neural crest. This is also the first analysis of the iridophore transcriptome. Gene expression analysis for iridophores revealed extensive enrichment of specific enzymes to coordinate production of their guanine-based reflective pigment. We speculate the coordinated upregulation of specific enzymes from several metabolic pathways recycles the rate-limiting substrate for purine synthesis, phosphoribosyl pyrophosphate, thus constituting a guanine cycle. The purification procedure and expression analysis described here, along with the accompanying transcriptome-wide expression data, provide the first mRNA sequencing data for multiple purified zebrafish pigment cell types, and will be a useful resource for further studies of pigment cell biology.


Introduction
Pigment cells serve as useful models for understanding many aspects of developmental and cell biology. For example, melanocytes are pigment cells studied to understand cell specification, migration, differentiation, survival, regeneration, organelle transport, secretion, and disease [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. Melanocytes produce melanin, which in humans serves as a UV protectant in skin [16]. Melanocytes also have roles in other organs such as the ear, brain, heart, and adipose tissue [17,18]. The incidence of melanomas, a disease of melanocytes and the most lethal form of skin cancer, is also increasing [19]. Methods to isolate and culture melanocytes for in vitro studies have been informative for understanding melanocyte biology [20,21]. In vivo studies of melanocyte biology and melanoma dynamics have been aided by the identification of mutants in mice and zebrafish [22][23][24][25]. Given the utility of zebrafish melanocytes to understand cell biology and disease, the transcriptome-wide characterization of genes expressed in zebrafish pigment cells would be a significant resource.
In mammalian systems, melanocytes are the only neural crestderived pigment cell type found in the dermis. In contrast, several neural crest-derived pigment cells are found in zebrafish and other poikilotherms, including reflective iridophores [5]. Several re-quirements for iridophore development from the neural crest are known [5,26,27]. However, it is unknown if markers of neural crest identity persist in iridophores following development, and whether these markers are shared by other neural crest-derived pigment cells, such as melanocytes. A further question in iridophore biology is how the guanine-based pigment is produced [28,29]. Zebrafish bearing mutations in the de novo purine synthesis enzymes gart and paics have iridophore defects, indicating purine synthesis is important for iridophore pigmentation [30]. Identifying a possible mechanism by which iridophores produce an abundance of guanine for pigment formation while maintaining adequate supplies of purines for DNA and RNA production will be informative for cell biology.
Another pigment cell type shared by mammalian and poikilothermic vertebrates is the retinal pigmented epithelium (RPE). The RPE is a group of melanin-producing cells found in the vertebrate eye. The RPE develops from the eye primordium, and is continuous with the layer of cells that forms the iris [31]. The RPE is critical for eye development and retinal health. It provides trophic support and recycles wastes from the photoreceptors of the retina [32]. The RPE forms part of the blood-retina barrier, providing the eye with an immune-privileged status [33]. Defects in the RPE contribute to diseases such as macular degeneration and retinitis pigmentosa, which result in vision problems [34,35]. Previous descriptions of gene expression in the RPE of zebrafish, chicken, and human have elucidated many of the genes playing roles in RPE biology [36][37][38]. Many of these genes are responsible for producing melanin, and defects in melanin production are often associated with reduced visual function [39]. However, it is unknown whether RPE and melanocytes use different pathways of melanin production, or if they are essentially identical. This information would be useful for understanding RPE biology, and would also inform future examinations of regulatory control for genes expressed in one or both cell types.
In order to facilitate understanding of these pigment cells, we developed a robust method to isolate these three pigment cell types from zebrafish embryos, followed by mRNA sequencing and transcriptome analysis. This purification procedure relies on the inherent densities of melanin and guanine-filled cells; hence it can be used without other complicated lineage markers. Here, we report the co-enrichment and cell-type specific gene expression profiles of melanocytes, iridophores, and RPE from embryonic zebrafish. While RPE and iridophores do not exhibit significant overlap of enriched gene expression, our analysis reveals considerable overlap among pairs of pigment cell types indicative of their common origin or function. Genes enriched in both the melanocyte and RPE lineages contain genes in the melaninproduction pathway, and suggest a more complete picture of melanin-producing machinery is contained within this set. Similarly, expression co-enrichment in melanocytes and iridophores are reflective of their neural crest origin, and suggest genes in this set may be specific to neural crest identity. Furthermore, this is the first characterization of the iridophore transcriptome. We found that iridophores specifically upregulate the guanine portion of de novo purine synthesis, as well as specific enzymes from other metabolic pathways that aid in producing the iridophore guanine-based reflective pigment. In addition to the analysis presented here, this procedure and accompanying data provide a significant resource for further biological discovery in pigment cells.

Zebrafish Strain and Sample Collection Time Points
This study was carried out in accordance with the Washington University Animal Use Committee guidelines under approved protocol #20110236. Zebrafish were reared and bred according to standard protocols [40]. The fish used in this study were homozygous for a temperature sensitive allele of micropthalmia transcription factor (mitfa vc7 ) [41]. This mutant facilitated the collection of RPE, as mitfa is not required for RPE development in zebrafish [41]. Melanocytes, iridophores, and RPE develop normally at 25uC in mitfa vc7 . When held at 32uC, the neural crest-derived melanocytes do not develop, but RPE and iridophores develop normally. All melanocyte and iridophore samples were incubated at 25uC prior to collection. We note the possibility that mitfa may be partially compromised by aberrant splicing products in this mutant at permissive temperatures, but following development at 25uC, melanocyte numbers, morphology, and pigmentation are indistinguishable from wild-type zebrafish [41]. Embryos used for RPE samples were incubated at either 25uC or 32uC until the equivalent of 3-5 days at 28.5uC, as indicated in Table S1 [42]. It should be noted that at stages later than 5dpf, choroidal melanocytes may be present adjacent to the RPE. Only one of the five RPE cDNA libraries was prepared later than 5dpf, as indicated in Table S1. We did not have separate markers to identify contamination from eye-associated neural   Table S3 for primer sequences. MMLV reverse transcriptase adds cytosines to the 39 end of the 1st strand cDNA (in black), allowing for template switching and addition of the 39 universal primer (B). PCR amplification of the library is followed by RsaI and AluI enzymatic digestion of cDNAs (C), followed by the standard Illumina library preparation steps of end-repair, a single adenine addition, Y-adapter ligation (D), PCR enrichment, and size selection (mock gel shown in E with yellow box indicating area of gel removed for DNA extraction), prior to flowcell generation and sequencing. doi:10.1371/journal.pone.0067801.g002 crest-derived melanocytes, or choroidal melanocytes in our RPE preparations. If present, we expect the contribution to total gene expression in the RPE samples to exhibit high variance, and found to be insignificant by Student's T-test upon comparison to the other cell types. However, we find the later stage RPE library to be highly correlated overall with the earlier stages, indicating that extensive contamination from choroidal melanocytes is unlikely (r = 0.92, Figure S1a). It was a further possibility that many large changes in gene expression would be present in RPE samples held at restrictive and permissive temperatures for mitfa. We found that samples held at the low and high temperatures exhibited a high degree of similarity (r = 0.89), indicating that most genes are not significantly different at the two temperatures ( Figure S1b). We also expect some genes to change during development between 3 and 5 days post fertilization. Upon inspection of melanocyte samples prepared at several time points, we find correlations increase with increasing developmental age of the embryos ( Figure  S2). We do not have the statistical power to confidently track developmental changes in gene expression within specific cell types, but we make the data available for all individual libraries to facilitate further investigation (Table S2). Figure 1 depicts the general procedure for pigment cell purification. Fish were anesthetized with Tricaine, rinsed with Ca-, Mg-DPBS (Sigma, D8537), and immersed in 100 mL TrypLE Express (Invitrogen, 12604039) per 1000 fish. Fish were incubated at 37uC and shaken at 100rpm for 15-20 minutes, followed by trituration with a Pasteur pipette to remove eyes from larva. After separation of eyes and larva, each group was placed in TrypLE Express and shaken at 100rpm at 37uC for 1-1.5 hr. Dissociated cells were filtered through a 120 uM screen into 50 mL tubes. Remaining intact tissue was triturated 10-20 times, and again filtered through a 120 uM screen into the dissociated cells. Dissociated cells were pelleted in a swinging bucket rotor (Eppendorf 5810 R) at 500 relative centrifugal force (rcf) for 5 minutes at 4uC, then resuspended in 1 mL cold isotonic Percoll (Sigma, P1644) by gentle pipetting. Isotonic Percoll was prepared by mixing 1 part 10X PBS with 9 parts Percoll. Resuspended cells were transferred to 1.6 mL Eppendorf tubes and spun at 2000rcf for 5 minutes at 4uC in a swinging bucket rotor for isopycnic separation. Pigment cells in the pellet were then resuspended in 400 mL of ice cold DPBS with 2% fetal calf serum (FCS), and placed onto preformed Percoll density gradients. Preformed gradients were prepared via centrifugation of 1 mL aliquots of isotonic Percoll in 1.6 mL tubes at 10,000rcf for 15 minutes at 4uC in a fixed angle micro-centrifuge (Eppendorf 5415 R). Tubes containing preformed Percoll gradients with overlying cell suspensions were centrifuged in a swinging bucket rotor at 2000rcf for 10 minutes at 4uC. Following centrifugation, overlying Percoll was aspirated, leaving the final 100 mL containing the pigment cell pellet. Cells were resuspended with 50 mL of cold DPBS with 2% FCS and transferred to a clean 1.6 mL tube containing 500 mL of cold DPBS with 2% FCS, and kept on ice until mRNA extraction or FACS.

FACS
We used the inherent properties of the pigmented cells to perform Fluorescence-Activated Cell Sorting (FACS). Following resuspension in 500 mL cold DPBS with 2% FCS, the enriched cell populations were screened through a 30 uM cell filter (Partec, 04-0042-2316). Cells were analyzed and sorted with a Dako MoFlo cell sorter using a 120 uM nozzle at a drop drive (DD) frequency of 22390 Hz. Cells were illuminated using a 488 nm laser. Cells were gated on two attributes to separate cells from each other and from cellular debris. Cellular debris was detected using forward and side scatter, selecting against the smallest particles (,1 mm or less). Cells were sorted based on detection using 510-530 nm and 575-595 nm filters, corresponding to FL1 and FL2 in Figure 1D, respectively. When excited by the 488 nm laser, the autofluorescence of iridophores is clearly detectable in these channels as a group of cells extending at a 45 degree line in the upper right quadrant. Melanocytes and RPE do not autofluoresce with this intensity when excited by the 488 nm laser, and cluster at the lower left of the FACS plot. Cells were collected into ice-cold DPBS with 2% FCS and kept on ice until mRNA extraction.
mRNA Extraction, cDNA Synthesis, and Illumina Library Preparation Pigment cell cDNA library construction was as follows. For mRNA extraction the DynabeadsH mRNA DIRECT Kit (In-vitrogen) was used per manufacturer's instructions. Following mRNA elution from the Dynabeads, first strand cDNA synthesis was performed using MMLV reverse transcriptase (Clontech) using an anchored polyT primer tailed with a universal primer sequence (See Table S3 for primer sequences and Figure 2 for pigment cell cDNA library construction overview.) A universal primer sequence was also added to the 39 end of the first strand by template switching, allowing for PCR-amplification of the resultant cDNA [43,44]. Following PCR amplification using the high fidelity polymerase LA Taq (TaKaRa, PCR cycle: 95C for 1 minute, followed by 20 cycles of 98C for 25 seconds, 60C for 1 minute, 68C for 20 minutes), cDNA was digested with AluI and RsaI restriction enzymes (NEB). Blunt-end enzymatic fragmentation of cDNA was used instead of sonication and gel extraction to minimize loss of sample material and eliminate the end-repair step of Illumina library preparation. Since this reduced representation strategy might miss short cDNAs that lack both restriction sites, we sought to avoid this by including enzyme recognition sites within the cDNA amplification primers. This allows for the inclusion of short cDNAs in our libraries. Standard Illumina library preparations followed, performed by the Genome Technology Access Center (GTAC) at Washington University in St. Louis (http:// gtac.wustl.edu). In brief, a single A was added to the 39 end of each strand, Y-adapters ligated, and library enrichment PCR performed, followed by gel extraction size-selection for fragments ranging from 200-400 base pairs in length. Illumina library construction of pooled 3dpf embryos was performed by GTAC from total RNA extracted with Trizol reagent as previously described [45]. No PCR amplification of whole embryo cDNA was performed prior to Illumina adapter ligation and library enrichment. Sequencing was performed on the GAIIX and HiSeq 2000 Illumina platforms. Technical sequencing replicates of the same libraries on separate lanes were essentially identical ( Figure  S3).

Sequence Analysis
We used Novoalign (www.novocraft.com), to assign resultant expression sequence tags to a customized non-redundant database of cDNA sequences consisting of 25,102 known and predicted genes (Table S4). Initial inspection of the NCBI zebrafish mRNA database of 28,286 zebrafish cDNAs (ftp://ftp.ncbi.nih.gov/ refseq/D_rerio/mRNA_Prot/) via an all-by-all BLAST search revealed multiple nearly identical sequences for cDNAs that would confound the unambiguous assignment of sequence tags. In order to generate a non-redundant cDNA database we selected single representatives for each gene as follows. In instances where a cDNA in the database resulted in a BLAST hit of greater than 94% identity to more than 70% the length of another transcript, we excluded the smaller of the two cDNAs. The resultant nonredundant zebrafish cDNA database contained 25,102 unique gene records. We further analyzed our non-redundant database by mapping all cDNAs onto the UCSC zebrafish browser (ZV9). Manual examination of 6 Mbs arbitrarily chosen along chromosome 13 revealed 130 annotated genes (combined RefSeq and Ensembl gene tracks). Of these, our non-redundant database identified 126. The 4 genes not represented in our database included 3 5 S ribosomal RNA genes and EN-SDARG00000086970, an annotated gene with a predicted ORF but no clear orthology to other species. Thus, this test shows that our method to generate a non-redundant cDNA library results in a database that identifies ,99% of annotated genes (excluding ribosomal RNA genes). Furthermore, manual examination of our non-redundant library mapped onto this 6 Mb of chromosomal sequence revealed 10 sequences that were not annotated as RefSeq or Ensembl genes. Comparison to repeat tracks on the UCSC browser revealed that most (8/10) of these unannotated matches identified ORFs from repetitive, or retrotransposon DNA. It is not clear how these sequences were initially included in the NCBI cDNA database. In summary, our efforts have generated a non-redundant zebrafish cDNA database that identifies most (,99%) of annotated genes, and includes 7% (10/126) records of dubious utility, but would tend not to confound RNA-seq analysis such as reported here. Sequencing results for each cDNA library are summarized in Table S1. For each cDNA library, the number of tags aligned to each gene was normalized by the length of the gene and the total number of uniquely aligning reads for that library using custom Perl scripts (reads per kilobase of cDNA per million mapped reads -RPKM [45]. Statistical calculations were performed in R (www.r-project.org). All sequencing data used in this study are publicly available at NCBI's GEO database under accession GSE46387 (http://www.ncbi.nlm.nih.gov/geo/query/ acc.cgi?acc = GSE46387).

Quantitative RT-PCR
Primers used for expression analysis are listed in Table S3, designed with NCBI's Primer-BLAST to be separated by at least one intron (http://www.ncbi.nlm.nih.gov/tools/primer-blast/). Expression data was normalized relative to beta actin for each sample, using cDNA produced as described above, without PCR amplification. Q-PCR was performed in a Perkin-Elmer thermocycler with the following conditions: 95C 1 min, 98C 20 sec, 60C 1 min, repeated for 40 cycles.

Results and Discussion
Purification of Melanocytes, Iridophores, and RPE from Whole Embryos to Generate Cell-specific Gene Expression Data At three days post fertilization (dpf), melanocytes, iridophores, and the RPE are readily visible in zebrafish. Melanocytes are extensively dendritic, and identifiable due to the presence of black melanin. Iridophores are round reflective cells, easily seen with a microscope using incident light. The RPE is present in the eyes as a hexagonally packed layer of melanized cells. However, the percentage of melanocytes, iridophores, and RPE compared to all other cells in the fish is less than 1%, making cell-specific gene expression analysis from the whole organism difficult. Aided by a previously reported method using density gradient centrifugation to isolate melanocytes from the caudal fins of Blackmoor goldfish [20], we developed a simple procedure to rapidly purify these pigmented cells from zebrafish embryos. This procedure relies on the inherent densities of the melanin-filled melanocytes and RPE, and the guanine-filled iridophores. In brief, cells are enzymatically dissociated from intact fish, enriched for pigment cells via density gradient centrifugation, and sorted by flow cytometry using the autoflourescent property of iridophores (Figure 1). The other pigment cells in zebrafish, the pteridine-containing xanthophores, do not pass through the density gradient and are not isolated by this method. Following cell purification, mRNA is isolated and Illumina libraries are constructed from cDNA ( Figure 2). We purified melanocytes, iridophores, and RPE as described above and performed mRNA sequencing for 11, 5, and 5 independently isolated samples, respectively. In order to confirm the quality of our cell specific libraries, we assembled a list of genes known to be expressed in melanocytes, iridophores, and RPE. When we compared expression of these genes in our cDNA libraries to this control list, we observed expression corresponding to each cell type (Table 1). Furthermore, we also found good correlation for fold changes determined by qPCR of independently prepared biological samples and those determined by our RNA-seq data, indicating our RNA-seq fold change calculations between pigment cell types are realistic values ( Figure S4). For instance, ltk, known to be a marker of iridophores, is highly enriched in iridophores compared to melanocytes (red square with x in Figure S4) [27]. Similarly, we find dct, a known requirement for melanogenesis, to be highly enriched in melanocytes and RPE relative to iridophores. Interestingly, we also find rpe65a, well-known to be expressed in RPE, to also be expressed by melanocytes. This is not entirely surprising, as RPE65 is known to be present in keratinocytes, melanocytes, and melanoma, in addition to the RPE [46,47]. Thus, using this method, cDNA from melanocytes, iridophores, and RPE can be generated concomitantly from thousands of whole zebrafish embryos in a single day. The development of a fast and robust method for purifying these pigment cells greatly simplifies the production of multiple biological replicates needed for informative analysis of highthroughput expression data.

Analysis of Gene Expression in Melanocytes, Iridophores, and Retinal Pigmented Epithelium
We first aimed to use our RNA-seq data to identify the total number of genes expressed in each of the pigment cell types. To eliminate low-level background expression, we applied a baseline expression threshold of 1 read per kilobase of transcript per million reads (1 RPKM). We chose this threshold based on the ability to detect expression of known pigmentation and neural crest genes. For example, we found mc1r, kita, and foxd3 to have RPKM values in melanocytes of 2.2 and 3.8, and 3.0, respectively (Table 1). Furthermore, we detected more than 95% of the genes expressed at 1 RPKM with a sequencing depth of approximately one million reads per library, which established our sequencing depth threshold ( Figure S5). Using this baseline threshold of 1 RPKM, we found 8,472 genes are expressed by melanocytes, 8,096 by iridophores, and 9,053 by the RPE (See Table S5 for the averaged RPKMs and T-test values). To put this into perspective, we also aligned over 700 million cDNA tags generated by the Stemple laboratory [48] from a variety of embryonic stages and adult tissues. We found this mixed data set revealed expression of at least 1 RPKM for 20,548 gene entries from our database of 25,102 unique coding sequences (Not shown). These results indicate our non-redundant database is a reasonable approximation of the protein-coding transcriptome in zebrafish, and conclude that 30-40% of all genes are appreciably expressed (.1 RPKM) by these specific cell types.

Melanocyte, Iridophore, and RPE Gene Expression is Correlated Compared to Whole Embryos
Having obtained gene expression data for these pigment cells, we set out to identify signatures of pigment cell functions from the data sets. We reasoned that genes co-enriched among pigment cell types would indicate shared pigment cell functions. As an initial Figure 3. The guanine synthesis cycle is highly enriched in iridophores. Shown is a model for guanine production based on transcriptome data as given in Table 6. Genes that are statistically enriched compared to melanocytes are shown in bold, those not statistically different are in normal font. The arrow thicknesses correspond to the fold changes in iridophores relative to whole embryos. Chemical structures are from the KEGG Compound database. doi:10.1371/journal.pone.0067801.g003 assessment of the similarity of gene expression between melanocytes, RPE, and iridophores, we applied Pearson's Product-Moment Correlation test to the datasets. Because correlation values can be artificially skewed by outlying data points [49], we calculated the average correlation from 24,102 overlapping windows of 1000 genes, sorted by increasing whole embryo expression ( Figure S6). Using this metric, we found the melaninproducing melanocytes and RPE were the most highly correlated (r = 0.90). The neural crest-derived melanocytes and iridophores were modestly correlated (r = 0.52). RPE and iridophores were also modestly correlated (r = 0.49). We expected the correlation to be dramatically lower when comparing a single purified cell type to whole embryos than when comparing two purified pigment cell types to each other. We thus determined the correlations of each pigment cell type to mRNA-seq data from whole 3dpf zebrafish. By this general assessment, melanocytes, RPE, and iridophores are distinct from whole embryos, with average correlation values less than 0.02. Presumably, this low value (0.02) reflects the baseline correlation component from expression of housekeeping genes shared by all cells, and greater values (0.4-0.9) reflect shared specific gene expression among the pigment cell types.

Gene Expression Indicative of a Pigment Cell Identity
These correlations described above suggest there are genes enriched in all three pigment cell types that are generally indicative of pigment cell identity. It is not clear a priori which genes should be shared by different pigment cell types. In order to construct an informative set of genes for this purpose, we searched for those that are expressed within a 2-fold change of each other at a minimum of 4 RPKM, and at least 100-fold greater than whole 3dpf embryos. We found 28 genes fulfilled these criteria (Table 2). Remarkably, 4 of the top 5 highest expressed genes on this list are ribosomal proteins. Pigment cell enrichment of ribosomal components is not surprising, considering several mouse coat color mutants are ribosomal proteins, although they are not the ribosomal components enriched here [50][51][52][53]. We also find the zebrafish albino gene (slc45a2) to be among the highest expressed genes in each of these pigment cell types [54,55]. The expression of slc45a2 in iridophores is interesting, considering albino fish are not reported to have an iridophore defect. One possibility is that SLC45A2 performs a common role for organelle pH homeostasis in pigment cells. However, in situ analysis reveals no enriched expression in xanthophores, suggesting that slc45a2 expression is not shared by all pigment cells [55]. This list of co-expressed genes also contains several other unexpected members, including the Jak-STAT cytokine receptor crfb5, and the BAX-inhibitor protein ghitm. It is not clear what roles these genes play in pigment cells, but this list provides a starting point for understanding their shared functions in pigment cell biology.

Gene Expression Indicative of Cellular Function or Developmental Origin
We were also interested to identify genes that were enriched in only two pigment cell types that would reveal shared function or developmental origin. For this analysis, we filtered for shared expression of genes in two cell types at least 2-fold over the third cell type, where differences exceeded a significance of p,0.05. Upon validation of expression differences via qPCR between pigment cell genes and whole embryos, we found a systematic bias of overcalling the fold changes between pigment cell values and whole embryos ( Figure S7). Based on this result, we also required an 8-fold greater than embryo expression threshold for coenrichment and cell-specific genes, as discussed below. Using these criteria, we found 214 genes were enriched in both melanocytes and RPE but not iridophores (Table S6), and 62 genes enriched in melanocytes and iridophores but not RPE (Table S7). Given that the RPE and iridophores do not share a pigment-type production or developmental origin, we expected fewer genes to be coenriched in these cells, but not expressed by neural-crest melanocytes. Only one gene, alcama (also dm-grasp or neurolin-a), the target of the zn-5/8 monoclonal antibody [56,57] was enriched in both RPE and iridophores when compared to melanocytes and whole embryos (Table S8). This is consistent with the notion there are few specific functions shared only between iridophores and RPE. It is intriguing that alcama has been found to mediate endothelin 1 signaling in cartilage development [58]. Relatedly, endothelin receptor (ednrb1) signaling is required for iridophore development in zebrafish [5]. It will be interesting to know whether alcama mutations have functional consequences in the iridophore or the RPE. Therefore, we suggest these lists of co-expressed genes are likely enriched for common metabolic functions, in the case of melanocytes and RPE, or developmental origins, in the case of melanocytes and iridophores.

Melanocyte and RPE Co-Expression
Due to their shared function of producing melanin and melanosomes, we expected many genes to be shared between melanocytes and the RPE when compared to iridophores. We found 214 genes that fit our parameters for shared enrichment (Table S6). As expected, there are many genes present involved in melanin synthesis and melanosome biogenesis, including pmela, dct, tyrp1b, tyrp1, pah, and slc24a5. There are several transcription factors in this enrichment group, including three forkhead box (foxo1b, foxp4, and foxg1b) and five homeobox-containing transcription factors (hmx1, hmx4, otx1a, otx2, and otx5). Although not required for RPE development in zebrafish, mitfa is expressed in the RPE at a relatively high level, consistent with the reported ability of mitfa to promote pigmented fate in zebrafish retinas [59]. However, as we find most of the other transcription factors coenriched in melanocytes and RPE are expressed at lower levels, it is interesting to speculate that one or more of these factors are able to compensate for mitfa specifically in the RPE. Thus, this coenrichment set may identify a more comprehensive list of genes for melanocyte and RPE identity and melanosome biogenesis.
An additional finding apparent from the 38 most highly expressed genes co-enriched in melanocytes and RPE is that iridophores also express these genes, albeit at a much lower level ( Table 3). The close developmental lineage relationship between melanocytes and iridophores suggests a possible ''leakiness'' of specificity that may result in weak expression of melanocyte genes in iridophores. This possibility is supported by the observation that the parade zebrafish mutant contains pigment cells with both melanocyte and iridophore characteristics [22].

Melanocyte and Iridophore Co-Expression
Because of the common developmental origin from the neural crest for melanocytes and iridophores, we speculated that neural crest-specific gene expression would be readily identifiable. We found 62 genes that were significantly upregulated in melanocytes and iridophores over RPE (Table S7). A list of the 15 genes expressed at least 5-fold greater than RPE and 10-fold greater than whole embryos is shown in Table 4. Included in this group are several well-known regulators and markers of neural crest and pigment cell development, including the transcription factors sox10 and the expressed repetitive element crestin [60]. The well-known pigmentation gene mc1r is also expressed in iridophores and melanocytes. Also included is the cell adhesion molecule pcdh10a, which is expressed by migrating zebrafish neural crest cells, and acts as a tumor suppressor in several cancers [61][62][63]. We also find a retinoic acid nuclear receptor subfamily member (rxrga) previously reported to be expressed in neural crest tissues [64][65][66][67]. Because we find this known set of neural crest genes in this coenrichment list, other unknown markers of neural crest identity are likely to be present. For instance, several transcription factors are in this enrichment set not previously reported in neural crest, including the forkhead box transcription factor foxo1a, as well as the cell cycle regulator cdk15. Thus, this list of genes co-enriched in melanocytes and iridophores may more broadly identify markers of neural crest origin.

Cell-Type Specific Gene Expression
We were also interested in determining gene expression specific to each cell type. Because melanocyte and RPE gene expression have been previously characterized, we do not discuss them in detail here, but present the data as the supplemental tables listed below [36][37][38]68,69]. For this analysis we required an expression level at least 2-fold over the other two cell types (p,0.05), as well as 8-fold greater expression than whole embryos. This filtering strategy resulted in a list of 108 genes specifically enriched in melanocytes (Table S9) and 24 in the RPE (Table S11). To independently validate these enrichment sets, we generated separate biological samples for each cell type and compared expression of selected genes from the enrichment sets via qRT-PCR. We observed good correlation of RNA-seq and qRT-PCR relative expression levels, indicating our RNA-seq data is a reliable indicator of the gene expression in these pigment cells (r 2 = 0.75, Figure S4). These data will be useful for further studies of melanocyte and RPE biology.

Iridophore Gene Enrichment
In contrast to melanocytes and RPE, the iridophore transcriptome has not previously been explored. Since iridophore produce a guanine-based pigment, rather than the melanin characteristic of melanocytes and RPE, we expected to find many genes to be specifically enriched in this cell type. This indeed turned out to be the case, with 346 genes passing our baseline threshold for enrichment (Table S10). Included in this enrichment list are several factors known to be important for iridophore development, including ltk, ednrb1, and pnp4a [5,26,27]. Also, in order to identify previously unreported genes that may play interesting roles in iridophores, we filtered our list to include genes expressed at least 30-fold greater than melanocytes and RPE, and 100-fold greater than whole embryos. Thirty genes met these criteria ( Table 5). The third highest expressed gene on this list, slc23l, may act as a guanine transporter. In mammals, the SLC23 gene family has roles in transporting nucleobases, such as guanine, as well as vitamin C. It is not clear how this highly expressed iridophore gene identifies a unique role for vitamin C in the iridophore, but it is tempting to speculate a role in guanine transport, either to transport guanine into the cell, or perhaps to transport newly synthesized guanine into the reflecting platelet organelles. One surprise is the finding that gpnmb is highly enriched in the iridophore. Roles for GPNMB have been described for melanocytes, melanoma, and the pigmented iris in mammals, and it has been suggested to act both as a plasma membrane protein and a component of the melanosome [70]. Our finding that gpnmb is more highly expressed in the iridophore than the melanocyte raises the possibility of a similar function in iridophore reflecting platelet organelle biogenesis. Also in this iridophore-specific enrichment list are six uncharacterized genes. Their enrichment in iridophores may aid in identifying functions for these proteins. We speculate that many iridophore-enriched genes will be indicative of novel cell-specific biological functions.
Many of the genes previously known to be expressed by iridophores are components of the purine synthesis pathway, as iridophore pigment largely consists of stacks of guanine plates [30,71]. Accordingly, we found a dramatic enrichment of enzymes comprising the pathway of guanine metabolism, from extracellular glucose import and glycolysis, through de novo synthesis and purine salvage ( Table 6). When comparing iridophores to melanocytes and RPE, we find 5 facilitated glucose transporters, 7/11 steps of glycolysis, and 9/9 enzymes for de novo purine synthesis to be enriched. Given that iridophore pigment consists largely of guanine, one might expect the guanine pathway to be specifically upregulated. Consistent with this model, we found the split in the purine synthesis pathway at IMP to favor guanine production rather than adenine. The first guanine-specific enzyme, impdh1b, is expressed at a level 71-fold greater than melanocytes. In contrast, the first adenine-specific enzyme, adssl, is not significantly different from melanocytes or RPE, at 0.7-fold the level of melanocytes. Upon inspection of the known pathway of guanine production, we observed synthesis of guanine from GMP likely results in the recycling of 5-Phosphoribosyl 1-Pyrophosphate (PRPP), the ratelimiting substrate in purine synthesis (KEGG Pathway: dre00230). The enzyme responsible for this final step of guanine synthesis, prtfdc1, is enriched 198-fold over melanocytes (p,0.01). From our expression data, we suggest a model of guanine pigment production in iridophores that illustrates a cycle of guanine synthesis utilizing PRPP as a recycled carrier molecule (Figure 3). In this cycle, specific enzymes from glycolysis, the pentose phosphate pathway, serine/glycine metabolism, and the citrate cycle, are upregulated to coordinate the extensive guanine synthesis required for the reflective iridophore pigment.
Another question in iridophore biology is how the membranous platelets containing the reflective guanine crystals are formed. Our data suggests a likely contribution from ADP Ribosylation Factors (ARFs) and Rab GTPases. ARFs are a large family of ras-related GTPases that regulate membrane trafficking and organelle structure (For review see [72]). We find two ARF-related genes to be significantly enriched in iridophores compared to melanocytes and RPE, arf6 and arfip1. Interestingly, we also found two members of the ras-related oncogene family to be enriched in iridophores; rab27b and rab38. Rab GTPases regulate many aspects of membranous vesicle formation and traffic [73]. In humans, Rab27 mutations cause hypopigmentation associated with the immunodeficiency disorder Griscelli syndrome type II [74]. In addition, the formation of COPI transport vesicles is known to be mediated by the interaction of GAPDH with RAB2 [75]. Phosphorylation of GAPDH occurs through the Src/PI3K/ AKT pathway, often downstream of receptor tyrosine kinase activation [76]. We found a single receptor tyrosine kinase enriched in iridophores, ltk, which we have previously shown to be required for iridophore development [27]. It is thus interesting to speculate that combined roles exist for GAPDH in both producing reactive carbonyl species during glycolysis for guanine production, as well as in forming the organelles within which iridophore pigment is contained.
Iridophores are a relatively less-studied pigment cell than melanocytes and RPE, and aside from the requirements for foxd3, ednrb1, and ltk, not much is known regarding the transcriptional regulation of iridophore identity. One might expect a single master regulator of iridophore identity, analogous to mitfa in melanocytes, which is highly expressed, to be readily identifiable in our data. However, we did not identify a highly expressed single candidate for a master regulator of iridophore identity. Instead, we found many moderately expressed transcription factors with known mouse or human orthologues to be significantly enriched in iridophores. Of these, there are several candidates that stand out as possible regulators of iridophore identity. One member of the basic helix-loop-helix (bHLH) family of transcription factors is specifically enriched in iridophores, tfec. Known to be expressed in iridophores, TFEC forms hetero and homodimers with other bHLH members and can function as a transcriptional activator or repressor [77][78][79]. Interestingly, there is one other bHLH in enriched in iridophores when compared to melanocytes and RPE, but it did not meet the 8-fold requirement over whole embryos. This gene, mycl1a, is a paralog of the classic oncogenic protein MYC. A key component of iridophore identity is the upregulation of glycolysis and guanine synthesis enzymes. It is established that MYC upregulates glycolysis, DNA-synthesis, and nucleotide metabolism [80][81][82]. Another oncogene shown to regulate glycolysis and associated feeder pathways is Ets-1 [83]. The vets-erythroblastosis virus E26 oncogene homolog 1a (ets1a) is also highly enriched in iridophores. In mouse, tfec transcription is activated through multiple ets-binding domains in its promoter region, suggesting a conserved regulatory mechanism for ETS1A in tfec transcriptional regulation as well [84]. Further work will be necessary to determine whether TFEC and MYCl1A coordinate the upregulation of guanine synthesis enzymes in iridophores.
We also find five homeobox-containing transcription factors are enriched in iridophores: gbx2, cart1, alx3 alx4a, and alx4b. Known to have several roles in the embryo and an early specifier of posterior neural crest in Xenopus, gastrulation brain homeobox 2 (gbx2), is highly enriched in iridophores [85]. Remarkably, the other four homeobox transcription factors we identified are all aristaless-related; cart1, alx3, alx4a, and alx4b. Members of the Cart1/Alx3/Alx4 family of homeodomain proteins are known to regulate formation of skeletal elements in organisms ranging from sea urchin to mammals [86,87]. In humans, mutation and haploinsufficiency of ALX4 are associated with skull ossification defects [88,89]. It is not clear how these aristaless-related transcription factors might regulate iridophore identity, but recognizing their expression here will aid in understanding their function. Together, these iridophore-enriched transcription factors likely play key roles in regulating iridophore identity.

Conclusion
To contribute to the understanding of pigment cell biology, we have developed a method for rapidly and reliably purifying melanocytes, iridophores, and retinal pigmented epithelium from zebrafish embryos, followed by global gene expression analysis by mRNA sequencing. This work represents the first concomitant comparison of three pigment cell types and whole zebrafish embryos, which uniquely allowed us to identify co-expressed genes indicative of shared function or developmental origin, as well as those that are specifically enriched in single cell types. We thus have identified many genes previously not reported to be enriched in these pigment cell types. In particular, we discovered a dramatic upregulation of specific enzymes from several metabolic pathways that coordinate guanine synthesis in iridophores, along with many membrane-trafficking components and transcription factors that are likely critical for iridophore identity. This characterization of global gene expression data from multiple purified zebrafish pigment cell types will provide a resource for further biological analysis of these cells. Table S1 Sequencing Summary. Gross sequencing results for cDNA libraries used in this analysis, with their respective developmental time points at 28.5C. Unless indicated, fish used for each library were held at 25uC. The type of sequencing runs are indicated by the suffixes of technical replicate names: 351: single end 36, 433, 436, and 440: single end 42, 399s61 and 399s62: paired end 101. *Uniquely mapped to the Washington University non-redundant cDNA database. Tags that were not uniquely assigned to a gene typically mapped to polyA, genomic, and Illumina-adapter sequences. **Total RNA was used for Illumina library preparation of 20 pooled 3dpf embryos, resulting in a lower fraction of uniquely mapping sequences due to a large proportion of ribosomal RNA in library. Abbreviations: Mel -Melanocyte, Irid -Iridophore, RPE -Retinal Pigmented Epithelium, hpfhours post fertilization, dpf -days post fertilization. (XLSX)