Comprehensive RNA-seq Analysis to Evaluate the Pigmentation-Related Genes Involved in Albinism of Cichlid Fish, Aulonocara baenschi

The Nkhomo-benga peacock, Aulonocara baenschi is a species inhabiting Lake Malawi and its population has decreased by over 70% over the past 10 years. Since 2018, the species has been included in the IUCN Red list and its status has been set to the critically endangered (CR) level. And the fish is known to have rare but naturally occurring albino strain. Despite that, little has been done to develop its genetic resources such as sequence-based data and expression patterns of genes. Accordingly, this study performed high-throughput sequencing on the caudal fin of the two different-colored A. baenschi in order to identify their pigment-related genes and analyze the expression patterns of the genes on the transcriptome level. Through the Illumina sequencing platform, we obtained 77,035,702 and 88,689,756 reads from the caudal fin of the wild-type and the albino-type, respectively. In addition to that, using de novo assembly and various annotation databases, a total of 141,116 contigs, including 129,369 unigenes were obtained. In the comparative analysis of differentially expressed genes between the wild-type and the albino-type in terms of melanogenesis, the tyrosinase gene family involved in melanin synthesis tended to be less expressed in the albino-type than the wild-type. However, the expression of the critical proteins involved in the pigment synthesis of xanthophore was higher in the albino-type than in the wild-type. In addition, the yellow-orange pigment was more apparently expressed in the caudal fin of the albino-type than of the wild-type. This indicates that the external characteristics of the fish were in accordance with the analysis results on the transcriptome level. The results of this study will be the foundation of research on expression characteristics and gene networking analysis of the pigmentation-related genes of fish.


INTRODUCTION
Albinism, characterized by a lack of melanin pigment in the skin, hair, and eye [oculocutaneous albinism (OCA) or ocular albinism (OA)], is known to be associated with mutations in genes involved in melanin synthesis (Oetting and King, 1999). A human with albinism mostly has whiteor bright-colored skin and hair, and the color of the eyes varies from dull gray-blue to brown. However, a certain type of albinism affects only the color of the eyes of the affected individual (Fertl and Rosel, 2009). A lack of or abnormal distribution of pigment in the iris can cause the eyes to be seen as red or violet (White and Rabago-Smith, 2011). Similarly, albinism is a phenomenon which can be observed not only in a human body but also in the body of a fish. It has long been reported and until recently, studies have been conducted on various species with albinism (Johnson, 1968;Gong et al., 2019). Most cases of albinism are attributed to a decrease or inhibition of gene expression caused by missense, nonsense, frame-shift and deletion mutation of the genes involved in the melanin synthesis pathway (Imes et al., 2006;Nakajima et al., 2012). In addition, tyrosinase (TYR), tyrosinase-related protein-1 (TYRP1; 5,6dihydroxyindole-2-carboxylic acid), tyrosinase-related protein-2 (TYRP2; dopachrome tautomerase; DCT), solute carrier family 24 member 5 (SLC24A5), and microphthalmia-associated transcription factor (MITF) have been reported as genes involved in albinism (Yasumoto et al., 1997;Amiel et al., 1998;Lamason et al., 2005). Melanin is synthesized in melanosome which is a highly specialized membrane-bound intracellular organelle in melanocyte, and melanin induces the production of hydrogen peroxide and quinone intermediates (Seiji et al., 1963;Sulaimon and Kitchell, 2003;Solano, 2014). There are different kinds of pigment cells in various species of vertebrate animals. In the case of fish, more varied kinds of pigment cell types have been found than in other species. To date, several kinds of chromatophores have been reported including melanophore, cyanophore, erythrophore, iridophore, leucophore, and xanthophore in teleosts (Schartl et al., 2015). Due to the diversity of pigment cells, fish have a wide range of body color and make use of them. For example, some fish camouflage themselves and use nuptial color to attract mates (Kodric-Brown, 1998;Fujii, 2000;Leclercq et al., 2010). It has been reported that inbreeding in the same species is one of the causes of body-color diversity (Shikano et al., 2007;Ohshima et al., 2013;Hattori et al., 2020).
The Aulonocara baenschi belongs to African cichlids and inhabit Lake Malawi, which is home to a wide variety of vertebrate species and lacustrine fishes (Danley et al., 2012;Oliver, 2018). While males have a blue face and a very bright yellow body, showing their breeding color all year, females have a creamcolored body with brown stripes and some show yellow color. Because of their appearance, the A. baenschi is also called sunshine peacock or yellow peacock (Meyer and Riehl, 1985;Alderton, 2012). Lake Malawi is where a total of 500 to 800 species of cichlid with various traits live, following millions of multiple repeated colonization and hybridization events (Joyce et al., 2011). The population of the A. baenschi has decreased by over 70% over the past 10 years due to over-collection by the ornamental fish trade market (The IUCN Red List of Threatened Species 2018: e.T61057A148658365). The species has been included in the IUCN Red list and its status has been set to the critically endangered (CR). In addition, genetic resources of the A. baenschi available are mostly limited to opsin genes which are classified according to the wavelength of light, and little genetic information is currently known about the fish other than mitochondrial DNA partial fragments. Furthermore, the albino mutant A. baenschi has hardly been studied genetically. Hence, this study attempted to compare on a molecular level the wild and albino-type of the A. baenschi, a species of African cichlid fish, which is characterized by its colorful body. To that end, we performed high-throughput sequencing to identify the types of pigment-related genes of A. baenschi with albinism and analyzed the expression patterns of the genes. The findings were compared with those of other species with albinism. Also, for the comparative analysis of the albino and the wild-type, the transcription factors interacting with the pigment-related genes were examined and the expression patterns were analyzed. Furthermore, we revealed which expression of genes decreased or increased in the albino-type compared with those in the wild-type, and the functions of the genes were examined.

Ethics Statement
All the experimental procedures with the fish were performed according to the guidelines of the Institutional Animal Care and Use Committee (IACUC) of Gangneung-Wonju National University (GWNU-2019-24). And all protocols were approved by GWNU ethical committee. Furthermore, all the authors of this study have completed Animal Welfare and Ethics Course certification under the CITI program, research ethics and compliance training program.

Fish Sample Collection
The fish used in the experiment were in the juvenile stage and their average weight and length were 3.56 ± 0.12 g, 6.8 ± 0.65 cm in the wild-type (n = 5), and 4.39 ± 0.25 g, 7.37 ± 0.26 cm in the albino-type (n = 5), respectively. They were purchased at a local aquarium (Seoul, South Korea). To ensure morphological clarity of pigment melanin in melanocyte cells, the fish were kept in a recirculation system with no light at 23 ± 0.5 • C for 1 week (Rodgers et al., 2010). The fish were anesthetized with 200 ppm of clove oil (Thermo Fisher Scientific, Waltham, MA, United States) prior to their caudal fin being dissected. The images of the caudal fin tissues were taken using a stereo microscope (Motic, SMZ-171) with 1-2× magnification (Figure 1). Subsequently, the caudal fin was maintained at −80 • C for total RNA extraction and ethanol fixation for DNA extraction. The species used in the experiment was confirmed based on control, partial control region of mitochondrial DNA (Supplementary Figure S1; Hashem et al., 2020) and the 5flanking regions of red-sensitive opsin and ultraviolet-sensitive opsin genes (JF262699.1, JF262735.1).

Library Construction and Illumina Sequencing
Total RNA isolation was conducted using RNAiso Reagent (Takara Bio, Shiga, Japan) according to the instructions given by the manufacturer. Additionally, TaKaRa MiniBEST Universal RNA Extraction Kit (Takara Bio, Japan) was used to extract the total RNA. The quality of the total RNA was confirmed by using wavelength of the ND-1000 nanometer (Thermo Fisher Scientific, Waltham, MA, United States) and Bioanalyzer (Agilent, Santa Clara, CA, United States). The entire transcriptome next generation sequencing (NGS) libraries were constructed using TruSeq stranded mRNA prep kit (Illumina, San Diego, CA, United States). RNA-seq libraries were sequenced using Illumina NovaSeq 6000 on 101 × 2 stranded (pairedend) read module, and raw sequencing files were produced using Illumina base-calling software and CASAVA v1.8.2 with ASCII Q-score (offset 33). Read-through library adaptor sequences, low-quality sequences (limit = 0.05), ambiguous nucleotides (maximal 2 nucleotides) were removed using CLC Genomics Workbench version 11.0 (CLC bio, QIAGEN, Denmark).

De novo Transcriptome Assembly and Annotation
Since the genome assembly of the A. baenschi had never been constructed to which high through-put sequencing reads were to be mapped, a de novo assembly was constructed for RNA-seq reference. The following default parameters for CLC Genomics Workbench 11.0 were used: kmer size 45, bubble size 50 and minimum contig length 150 bp. Transcriptome completeness was evaluated using the Benchmarking Universal Single-Copy Orthologs (BUSCO) dataset for Vertebrata_odb9, with the default parameters (number of species: 65, blast e-value: 1.0E-3, mode: transcriptome) (Simão et al., 2015). An annotation based on BLAST× was conducted with the unigenes obtained as a result of the de novo assembly. The conditions of the BLAST× annotation were as follows: e-value cut-off was 1.0E-3, the word size of BLAST parameters was 6, HSP length cut-off was 33 and HSP-Hit coverage was 0. The annotation was performed using the NCBI non-redundant protein sequences database (NR_v5). Additionally, functional annotations were performed using the InterProScan v 5.34-73.0 (Finn et al., 2017), Gene ontology (GO) 1 and GO-Slim and EggNOG v 4.5.1 (Huerta-Cepas et al., 2016). The results were merged with those of the BLAST× annotation to construct a RNA-seq reference.

Differentially Expressed Genes (DEGs) Analysis and RNA-seq Data Validation by Quantitative Real-Time RT-PCR (qPCR)
To investigate the differences of the genes expressed in the caudal fin between the wild-type and the albino-type, the reads of each group were mapped to the RNA-seq reference assembly. The differentially expressed genes (DEGs) analysis was based on the count of the reads mapped and the selection of the DEGs was performed with the false discovery rate (FDR) p-value below 0.05. To functional analysis in differential expressed genes, the GO analysis was conducted at annotation mean level of 7, which was determined with GO distribution. We analyzed the functions of the DEGs in GO terms of three aspects: biological process (BP), molecular function (MF), and cellular component (CC). The expression patterns of DEGs were analyzed based on the fold change (FC) value using up-regulation (FC > 1.0) and down-regulation (FC < −1). The various expression patterns identified through DEGs analysis were divided into six categories provided by the website color genes 2 : components of melanosomes and precursors, melanocyte development, eumelanin and pheomelanin, melanosome construction/protein routing, melanosome transport and systemic effects. Furthermore, qRT-PCR was performed on 13 kinds of pigment-related genes/transcription factors in order to validate the values estimated by RNA-seq. After normalizing the concentration difference of the total RNA between the two groups, cDNA was synthesized using PrimScript RT Kit (Takara) which included random primers and oligo-dTs. Ribosomal RNA protein7 (RPL7) was used as an internal control gene to correct the concentration of each template in qRT-PCR (Øvergård et al., 2010). The pair of primer sequences used in qRT-PCR of RPL7 and pigment-related genes was constructed based on the sequence of the reads resulting from RNA-seq ( Table 1). The primers with qRT-PCR amplification efficiency between 1.91 and 2.02 determined using a cDNA serial dilution series acquired from a single sample, were used for RNA-seq validation. qRT-PCR reactions were performed using Thermal Cycler Dice real-time PCR system (Takara) and SYBR premix Ex TaqII Kit (Takara) with the following cycling conditions: 1 cycle of initial denaturation at 95 • C for 30 s, followed by 45 cycles denaturation of 95 • C for 5 s, and annealing at 60 • C for 30 s. All the data in qRT-PCR were represented as Mean ± standard error of measurement (SEM). A one-way ANOVA was performed using SPSS 25.0 software, and the level of statistical significance was set at p < 0.05.

Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway and Investigation and Analysis of Protein-Protein Interaction (PPI) Networks Analysis
The assembled contigs of the A. baenschi were annotated with Kyoto Encyclopedia of Genes and Genomes (KEGG) Automatic Annotation Server (KAAS) and a pathway analysis of pigmentrelated genes was conducted (Moriya et al., 2007). To investigate the expression of the genes involved in the pigment of epithelium and learn about the genes interacting with the expression, we used STRING v11.0 3 based on the zebrafish (Danio rerio) database (minimum required interaction score > 0.5). Furthermore, key enzymes of melanogenesis including tyrosinase (TYR), TYRP1, and DCT were analyzed, and various protein genes known to interact with the key enzymes were investigated in the aspect of the gene expression patterns using the RNA-seq data. Finally, an interaction networks map was constructed using multiple proteins used in the analysis, and the networks were visualized with Cytoscape v3.7.1 software (Shannon et al., 2003).

Sequencing, Assembly, Annotation
We conducted whole transcriptome sequencing with the A. baenschi caudal fin of both wild-type and albino-type, and the results were 88,689,756 reads (48% GC contents) and 77,035,702 reads (48% GC contents), respectively (  Figure S2). We carried out GO analysis based on a node score/value distribution on the genes expressed in the caudal fin of the A. baenschi.
The genes related to epithelium development were found the most in BP, the genes related to hydrolase activity in MF and the genes related to protein-containing complex in CC (Supplementary Figure S3). Furthermore, regarding the enzyme code distribution within the GO annotation, hydrolases, transferases and oxidoreductases were the ones hit the most in order of priority (Supplementary Figure S4). As for the KEGG annotation, the genes concerned with pathways in cancer, PI3K-Akt signaling pathway, human papillomavirus infection and MAPK signaling pathway were hit the most (data not shown).

Clustering Analysis of Differentially Expressed Genes (DEGs)
Differentially expressed gene analysis identified gene expression level differences while GO and KEGG analysis addressed the function of these genes. In terms of BP were both upregulated and down-regulated genes concerned with nucleic acid templated transcription (GO:0097659) and regulation of nucleic acid templated transcription (GO:1903506) accounted for the largest portion of the BP category ( Figure 2). As for MF, nucleoside-triphosphatase activity (GO:0017111) and zinc ion binding (GO:0008270) took were the most common. With regard to CC, microtubule (GO:0005874) and lysosomal membrane (GO:0005765) were the most common terms overall. KEGG analysis identified pathways in which DEGs are involved (Figure 3). We focused on the 20 functional KEGG pathways most represented by our DEGs. For DEGs whose expression either up-or down-regulated the greatest number of genes were involved in the metabolic pathways (ko01100). However DEGs whose expression up-regulated in the albino-type were also involved in the pathways in cancer (ko05200) and PI3K-Akt signaling pathway (ko04151), while, the DEGs whose expression down-regulated were found to be in biosynthesis of secondary metabolites (ko01110) as well as pathways in cancer (ko05200). Interestingly there were several pathways that were almost exclusively represented by DEGs whose expression upregulated, including focal adhesion (ko04510), Ras signaling pathway (ko04014), and axon guidance (ko04360). In the DEGs between the wild-type and the albino-type, the genes expressed the most were T-cell surface glycoprotein CD3 gamma chain followed by myosin light chain 1 and endogenous retrovirus group PABLB member 1 Env polyprotein-like ( Table 3). The genes expressed the least were interferon-induced very large GTPase 1-like followed by barrier-to-autointegration factorlike protein and TPA: putative transposase. The KEGG analysis revealed that DEGs with high expression between the wildtype and albino-type in the caudal fin were associated with viral infection and cancer pathway or found to take control of tight junction and calcium signaling pathway. Also noted, the DEGs whose expression decreased were found to be concerned with mineral absorption, melanogenesis, tyrosine metabolism and longevity regulation pathway. The findings indicate that most of the genes involved in the components of melanosomes and their precursors sharply decreased their expression rate in the albino-type (Figure 4). Except MITFA and SLC45A2, all the other genes decreased their expression more than 20-fold in the albino-type in comparison to the wild-type. Of the genes involved in melanocyte development, except for PAX6A, KRT1-19D, KRT1-C5, and RECQL4, some showed no significant differences in the expression of DEGs between the wild-type and the albino-type and others showed the increase of expression in the albino-type rather than in the wild-type. In addition, the genes involved in eumelanin and pheomelanin tended to be expressed more abundantly in the albino-type than in the wild-type, which is contrary to the expression patterns of the genes involved in components of melanosomes and precursors. Furthermore, the expression of transcription factors involved in xanthophore including CSF1RA, AOX3, and GCH was higher in the albino-type than in the wild-type. Of note, most of the genes involved in melanosome transport and systemic effects tended to be expressed more remarkably in the albino-type than in the wild-type.

Gene Networking and KEGG Analysis
The variation in pigmentation in epithelium is determined by various pigment cell types, and the degree of the expression of the genes related to melanogenesis determines albinism. Hence, it is important to analyze the interactional network of the genes. Accordingly, we performed an analysis of protein-protein interaction (PPI) networks based on studies with albinism related protein genes or transcription regulators including TYR, TYRP1, and DCT. The PPI network map consisted of 28 nodes and 159 edges. The average node degree was 11.4. and the average local clustering coefficient was 0.794. This suggested that the protein genes related to melanogenesis and the transcription regulators interacted with each other (Figure 5). In the PPI network map, DDC, GOT1, GOT2A, IL4I1, and TPO, proteins involved in melanogenesis but also in tyrosinase metabolism,   metabolism of amino acids and derivatives formed a different group from tyrosinase family genes. In addition, the expression of these genes tended to be slightly higher in the albino-type than in the wild-type, as opposed to the expression of tyrosinase family genes. At the center of the PPI network map was TYR, a connection protein, and TYR functioned as a bridge between melanogenesis and tyrosinase metabolism. That was clearly proven when the analysis was performed with the interaction score set at high confidence (0.7). TYR was involved in regulating the expression of other genes except MAPK1, MAPK3, and GCH2, whereas AOX5, EDNRB1A, FOXd3, SOX10, and TFE3B were involved in regulating the expression of TYR. The analysis of the functional enrichments in the PPI network indicated that in reactome pathways DCT, DDC, GOT1, GOT2A, IL4L1, OCA2, SLC45A2, TH, TYR, and TYRP1 which were involved in metabolism of amino acids and derivatives (DRE-71291) interacted with one another when they were expressed. Among those, DCT, OCA2, SLC45A2, TYR, and TYRP1 turned out to be involved in melanin biosynthesis (DRE-5662702) as well. As for the analysis of PFAM protein domains, it was revealed that DCT, TYR, and TYRP1 were involved in the common central domain of tyrosinase and interacted with one another. Regarding the analysis of INTERPRO protein domains and features, it was discovered that transcription factors MITFA, MITFB, and TFE3 were involved in the control of the expression of the genes related to melanogenesis. The analysis of SMART protein domains indicated that helix loop helix domain was present in transcription factors including MITFA, MITFB, and TFE3. In the melanogenesis pathway (map04916), the expression of α-melanocyte stimulating hormone (α-MSH), a pigmentationinducing hormone, was higher in the albino-type than in the wild-type (Figure 6). α-MSH secreted by keratinocyte was combined with melanocortin-1 receptor (MC1R), a receptor present in the membrane of melanocyte, to activate adenylate cyclase, a signaling protein. Also, the expression of protein kinase (PKA), a signaling substance, and cAMP responsive FIGURE 6 | The KEGG signaling map for melanogenesis and tyrosinase metabolism. The genes displayed with different colors according to FC value.
Frontiers in Marine Science | www.frontiersin.org element binding protein (CREB) which controls the expression of MITF was higher in the albino-type than in the wild-type. In the A. baenschi, MITF is expected to exist in two isoforms: MITFA and MITFB, and the expression patterns of the genes are as follows. The expression of MITFA was lower and that of MITFB was higher in the albino-type than in the wildtype. The genes involved in the formation of melanin including TYR, TYRP1, and DCT showed the same expression pattern as MITFA. In the tyrosine metabolism pathway (map00350), the expression of TYR, a key enzyme in the conversion of tyrosine into L-DOPA (3,4-dihydroxyphenylalanine) and DOPA quinone, was minimally observed in the albino-type. Analysis of KEGG pathway indicates that the caudal fin of albino-type fish produced little pheomelanin and eumelanin in the melanoma inside the melanocyte.

Validation qRT-PCR
To validate the expression patterns of DEGs, pigmentationrelated genes and transcription factors were selected for qRT-PCR validation. qRT-PCR validation was performed on the genes involved in components of melanosomes and their precursors, melanocyte development, melanosome transport, eumelanin and pheomelanin. Most of the genes which were validated with qRT-PCR aligned with the result of the analysis using the high-throughput sequencing data generated with a statistical significance. As expected, the expression of the tyrosinase gene family of the albino-type tended to decrease when compared with the wild-type, the control group (Figure 7). When compared with the wild-type, the expression of the genes including TYRP1, TYR, TYRL, PMELB, DCT, and SLC24A5 which were involved in components of melanosomes and their precursors decreased 348. 1-, 68. 8-, 28. 0-, 8. 7-, 18. 9-, 9.5-and 20.6-fold, respectively. The expression of the genes including EDNRB-L and MITFA which were involved in melanocyte development decreased 2.6and 2.3-fold, respectively. Meanwhile, the expression of MITFB increased 2.2-fold. The expression of the genes including MREG-L and MLPH which were involved in melanosome transport decreased 15.2-and 6.0-fold, respectively. And SLC7A11, a gene involved in eumelanin and pheomelanin, decreased 2.3-fold, whereas POMCA increased 1.4-fold.

DISCUSSION
Whole transcriptome NGS was performed to identify the pigment-related genes and analyze their expression patterns in the caudal fin of both the wild-type and the albino-type of the A. baenschi. A phenotypic difference of the caudal fin of the two types of the A. baenschi is related to the presence of the dark pigment of the melanosome in the melanocyte. The results of the analysis on the transcriptome level revealed that the expression of the tyrosinase gene family (TYR, TYRP1, TYRP2/DCT) which are key factors involved in the melanin synthesis in the melanosome was much lower in the albino-type than in the wildtype. The aforementioned expression patterns have been reported in research not only on Osteichthyes such as the catfish, tilapia, carp and sturgeon but also on vertebrate animals including the chicken and goat (Jiang et al., 2014;Zou et al., 2015;Zhu et al., 2016;Yu F. et al., 2018;Bhat et al., 2019;Gong et al., 2019;Xiong et al., 2020). TYR is known to be involved in melanin synthesis and a key enzyme in the hydroxylation of tyrosine to DOPA (3,4-dihydroxyphenylalanine) and the oxidation of DOPA to dopaquinone (del Marmol and Beermann, 1996). The expression of TYR plays a major role in the melanogenesis pathway involved in the induction of albinism. Also, the decreased gene expression caused by mutation in the gene exon/intron structure or mutation in the transcription binding sites in the promoter hinders the synthesis of pheomelanin or eumelanin (Beermann et al., 2004;Boonanuntanasarn et al., 2004). The expression of not only TYR found in the caudal fin but also TYRL belonging to a different clade from TYR remarkably decreased in the A. baenschi albino group compared with the wild group (Yu S. et al., 2018). In order to investigate the cause of the differential expression, we mapped NGS reads of both the albino and the wild-type to TYR and TYRL obtained through de novo assembly. The results indicate that no read was mapped to TYR and only two reads were mapped to TYRL in the albino-type. The low expression of TYR was not due to the mutation in the exon/intron of the TYR gene but because of the binding affinity of the transcription factor in the promoter which regulated transcription in the A. baenschi albino. In addition, a gene sequence analysis of NGS reads of the A. baenschi wild-type which were mapped to TYR genes confirmed the existence of transcript variants of the species. This finding agrees with the results of other relevant research on the chicken, macaque and medaka, whose TYR genes produce multiple transcripts with alternative splicing (Hidehito et al., 1994;Kim et al., 2017;Yu et al., 2019). MITF, LEF-1, and transcription factor binding to IGHM Enhancer 3 (TFE3) are known to be representative transcription factors involved in the activation of the TYR gene promoter (Yasumoto et al., 1997;Verastegui et al., 2000). MITF is a key transcription factor and the master gene for melanocytic survival which is involved in the activation of TYR, TRP-1 and DCT, members of the TYR gene family (Fang et al., 2002;Vachtenheim and Borovanský, 2010). MITF was present in two different isoforms of the A. baenschi, and the two types of MITF showed differential expression patterns from each other, which was the same tendency as that observed in the skin of the crucian carp . It is interesting to note that the expression of DCT, MITFA, and TYR genes involved in melanin synthesis is lower in the red-colored crucian carp than in the back-gray colored crucian carp, which tends to be correlated with DNA methylation levels (Zhang et al., 2017a). In the color pattern variability of skin pigment of the trout, the expression of MITF decreased to a greater extent in the light region than in the black spot, which is a similar pattern MITFA showed in the A. baenschi. However, MITFB showed a different expression pattern in the same species (Djurdjevič et al., 2019). MITF found in the chicken did not show any significant differential expression patterns between the white type and the black type . The regulation of MITF expression is known to be affected by SOX10, PAX3 and LEF-1, transcription factors, and the expression patterns of MITF and the three transcription factors are known to be similar to each other. In this study, however, MITFA which FIGURE 7 | Validation of the RNA-seq FC data by qRT-PCR. Target genes were normalized to the reference gene, RPL7.
had the same expression pattern with the TYR gene family had a different expression pattern from SOX10, PAX3, and LEF-1. In contrast, MITFB, a different clade from MITFA, had the same expression pattern with SOX10, PAX3, and LEF-1, which was reported in previous studies. PAX3A and SOX10 found in the crucian carp and SOX10 in the tilapia showed different expression patterns from those observed in this study. In contrast, SOX10 found in the trout showed a similar expression pattern with that observed in this study (Zhu et al., 2016;Zhang et al., 2017b). The expression patterns of LEF-1, DCT, MITF, and TYR found in human embryonic kidney 293T (HEK293T) cells are similar to one another, which was not aligned with the findings of this study . It is interesting to note that the expression patterns of the involved genes varied in accordance with the pigment cell type. For example, the main genes involved in melanocyte including DCT, MITFA, TYR and TYRP1 in the A. baenschi were less expressed in the albino-type than in the wild-type. Note that most of the main genes involved in xanthophore including CSF1RA and GCG were more expressed in the albino-type than in the wild-type. SOX10 is known to begin to be expressed in the neural crest cell before the onset of melanin synthesis, and then take control of the expression of xanthophore, melanocyte and iridophore. In this study, it is believed that SOX10 takes control of the expression of CSF1RA and is involved in the expression of MITFB since the expression of the genes related to xanthophore increased in the albino-type (Singh et al., 2016). In the case of PAX3, it is reported to be a major transcription factor involved in the formation of not only melanocyte but also neural crest and xanthophores. In this study, too, the expression of the gene increased by being involved in the control of the xanthophores-related genes (Minchin and Hughes, 2008). Even though the 5 -flanking region of MITFA and MITFB was not obtained and transcription factor binding sites were not examined, it is assumed that the increase of the expression of LEF-1 and PAX3 was because they were involved in the expression of MITFB just as SOX10 was. SLC24A5, found in the zebrafish and known to be responsible for lighter pigmentation caused by hypopigmentation of melanocytes in the skin and the retinal pigment epithelium, was also expressed very little in the A. baenschi albino-type compared with in the wild-type (Lamason et al., 2005). Also, MATP (SLC45A2) and MLPH involved in the transport of melanocytes of inorganic ions and charged and uncharged organic molecules were less expressed in the albino group, showing the same expression pattern as the one observed in the Japanese flounder albinism and the crimson snapper (Zhang Y.-P. et al., 2015;Wang et al., 2017).

CONCLUSION
This study identified pigment-related genes found in the caudal fin of the two types of the A. baenschi with two different body colors which belongs to the African cichlid and has been listed on the IUCN Red list and explored the expression patterns of the genes. When we observed both the wild-type and the albinotype, a phenotypic difference in melanin pigment between the two types was observed. On the transcriptome level, remarkable differential expression among the genes involved in melanin synthesis was observed, too. In addition, the genes in various cell types involved in xanthophore were less expressed than in the wild-type. In regards to the body color, the albino-type had a yellow-orange colored body affected by xanthophore. The expression pattern analyses of certain enzymes and transport membrane proteins involved in albinism examined in this study may contribute to facilitating better understanding of albinism in the aspects of its cause and phenomena. Also, the study results may be the foundation of research on dermatopathy in vertebrates including mammals and teleostei and could be used as genetic resources for research on African cichlids with various genetic traits.

DATA AVAILABILITY STATEMENT
The raw data and the sample information used in this study were submitted to NCBI. The Sequence Read Archive (SRA) rawdata can be accessed in NCBI with the following information: BioProject (PRJNA578449) and BioSample (SAMN13061232 and SAMN13061233).

ETHICS STATEMENT
The animal study was reviewed and approved by Institutional Animal Care and Use Committee (IACUC) of Gangneung-Wonju National University. Written informed consent was obtained from the owners for the participation of their animals in this study.

AUTHOR CONTRIBUTIONS
SL prepared the fish samples (including critical care of fishes), analyzed the gene-expression and data, and conducted the manuscript draft. HL helped the quantitative real-time PCR experiment and improved the manuscript. Both authors developed and supervised the experiment, and modified the manuscript.