A gene expression study of dorso-ventrally restricted pigment pattern in adult fins of Neolamprologus meeli, an African cichlid species

Fin sampling and melanophore counting

Four captive bred, adult individuals of N. meeli, two males and two females, were used in this study. The size of the fish was between 7–10 cm and their fins looked intact. The fishes did not show sex dependent differences or individual variation in their stripe patterns. Prior to the experiments, the fish were transferred to separate aquaria and fed on identical diets for two weeks. To obtain fin tissue biopsies, the fish were anesthetized in water containing 0.04 gram per litre of MS-222 and parts of their dorsal fin (D), anal fin (A), and dorsal and ventral parts of their caudal fin (Cd and Cv) were cut under a stereomicroscope (the clipped fin areas are specified by blue dashed squares in Fig. 1A). From each of the cut fin tissues, a piece was dissected and immersed in 0.5 mg/ml epinephrine solution (Sigma No. E4375) for 60 min in order to aggregate the melanosomes. The remaining fin tissue was dissected carefully based on its distinct colours (specified with red circles in Fig. 1A), and each piece was immersed in RNAlater (Qiagen) and stored frozen until RNA isolation. Throughout the paper, the individual fin tissues investigated in the study are addressed as “fin regions” and identified by fin type and a number referring to location, such that, for instance, A-1 stands for the distal region of the anal fin, A-2 identifies the middle region and A-3 the proximal region (Fig. 1A). Anatomically equivalent fin regions are grouped into distal, middle and proximal “classes of fin regions” (e.g., the distal class consists of the regions D-1, Cd-1, Cv-1 and A-1; Fig. 1A).

Photographs of the epinephrine treated fin tissues were taken with a camera mounted on a stereomicroscope (eyepiece micrometre x20, Leica). For each fin region, melanophores were counted in four inter-ray fin tissue sections and the area of the inter-ray tissue was measured based on pixel counts. Melanophore density in each fin region was calculated as number of melanophores/mm². Anaesthesia and fin biopsies were performed under permit number BMWFW-66.007/0013-WF/V/3b/2016 issued by the Federal Ministry of Science, Research and Economy of Austria.

RNA isolation and cDNA synthesis

Tissue samples of the individual fin regions (see above) were transferred from RNAlater to tubes containing TRI Reagent (Sigma) and 1.4 mm ceramic spheres, and homogenized by FastPrep-24 Instrument (MP Biomedicals, Santa Ana, CA, USA). RNA was extracted according to manufacturer’s Trizol protocol and dissolved in 30 µl RNase-free water. RNA samples were treated with DNase (New England Biolabs) to remove contaminating DNA. RNA concentration was measured by spectrophotometry using a Nanophotometer (IMPLEN GmbH, Munich, Germany). The quality of the RNA samples was evaluated in a R6K ScreenTape System on an Agilent 2200 TapeStation (Agilent Technologies) to ensure that the integrity number (RIN) of all samples was higher than 7. cDNA was prepared from 1,000 ng of RNA using the High Capacity cDNA Reverse Transcription kit (Applied Biosystems), according to the manufacturer’s protocol. Negative controls, i.e., reactions without addition of reverse transcriptase (-RT samples), were prepared to confirm the absence of genomic DNA. cDNA was diluted 1:3 times in nuclease-free water for further use in quantitative real-time PCR.

Gene selection, primer design and real-time qPCR

To validate suitable reference genes for accurate expression analysis, we selected 7 genes expressed in a variety of tissues and frequently used as reference gene candidates in qPCR studies of teleost fishes (Table S1) (Tang et al., 2007; Fernandes et al., 2008; Olsvik, Søfteland & Lie, 2008; Small et al., 2008; Zheng & Sun, 2011; Yang et al., 2013; Ahi et al., 2013). In addition, we selected 15 target genes that are known to be involved in adult pigmentation and/or stripe formation in zebrafish (Table 1 and Table S1). Later, we extended our expression analyses to six additional candidate genes which were predicted as potential upstream regulators of two differentially expressed genes in our study (Table S1).

The qPCR primers were designed within sequences conserved across African cichlids, based on recently released transcriptome data from a distantly related species, Oreochromis niloticus, and a closely related species from Lake Tanganyika, Neolamprologus brichardi (Brawand et al., 2014). The sequence alignment was conducted using CLC Genomic Workbench, version 7.5 (CLC Bio, Aarhus, Denmark) and locations overlapping the exon boundaries of the genes were determined based on the Nile Tilapia annotated genome sequences in the Ensembl database (http://www.ensembl.org/Oreochromis_niloticus). The qPCR Primers were designed on exon boundaries of the conserved regions using Primer Express 3.0 software (Applied Biosystems, Foster City, CA, USA) and checked for self-annealing, hetero-dimers and hairpin structures with OligoAnalyzer 3.1 (Integrated DNA Technology) (Table S1).

Real-time PCR was performed in 96 well-PCR plates on an ABI 7500 real-time PCR System (Applied Biosystems) using Maxima SYBR Green/ROX qPCR Master Mix (2X) as recommended by the manufacturer (Thermo Fisher Scientific, St Leon-Rot, Germany). Each biological replicate was run in duplicate together with no-template control (NTC) in each run for each gene and the experimental set-up per run followed the preferred sample maximization method (Hellemans et al., 2007). The qPCR was run with a 2 min hold at 50 °C and a 10 min hot start at 95 °C followed by the amplification step for 40 cycles of 15 sec denaturation at 95 °C and 1 min annealing/extension at 60 °C. A dissociation step (60 °C–95 °C) was performed at the end of the amplification phase to identify a single, specific product for each primer set (Table S1). Primer efficiency values (E) were calculated with the LinRegPCR v11.0 programme (http://LinRegPCR.nl) (Ramakers et al., 2003) analysing the background-corrected fluorescence data from the exponential phase of PCR amplification for each primer-pair and those with E less than 0.9 were discarded and new primers designed (Table S1).

Data analysis

Three different ranking algorithms, BestKeeper (Pfaffl et al., 2004), NormFinder (Andersen, Jensen & Ørntoft, 2004) and geNorm (Vandesompele et al., 2002), were employed to identify the most stably expressed reference genes The standard deviation (SD) based on Cq values of the fin regions was calculated by BestKeeper to determine the expression variation for each reference gene. In addition, BestKeeper determines the stability of reference genes based on correlation to other candidates through calculation of BestKeeper index (r). GeNorm measures mean pairwise variation between each gene and other candidates, the expression stability or M value, and it excludes the gene with the highest M value (least stability) from subsequent analysis in a stepwise manner. NormFinder identifies the most stable genes (lowest expression stability values) based on analysis of the sample subgroups and estimation of inter- and intra-group variation in expression levels (Ahi et al., 2013; Pashay Ahi et al., 2016).

The Cq values of the best ranked reference genes was used as Cq _reference in the ΔCq calculations. For the analysis of the qPCR data, the difference between Cq values (ΔCq) of the target genes and the selected reference gene was calculated for each target gene; ΔCq _target = Cq _target − Cq _reference. All samples were then normalized to the ΔCq value of a calibrator sample to obtain a ΔΔCq value (ΔCq_target − ΔCq _calibrator). For comparisons of gene expression involving all three classes of fin regions, one biological replicate of D-3 was arbitrarily chosen as calibrator sample. For comparisons restricted to anatomically equivalent fin regions, one biological replicate of D-1, D-2 and D-3 served as calibrator for proximal, middle and distal fin region samples, respectively. Relative expression quantities (RQ) were calculated based on the expression level of the calibrator sample (E^−ΔΔCq) (Pfaffl, 2001). A two-way ANOVA followed by post hoc Tukey’s HSD test was implemented for each target gene to compare RQ values among fins (averaged across regions) and fin regions. To assess similarities in the expression patterns of the target genes, Pearson correlation coefficients (r) were calculated in R (http://www.r-project.org).

To identify the potential upstream regulators, we performed motif enrichment analysis on 1 kb promoter sequences of two differentially expressed genes, based on the annotated genome of the Nile tilapia (Flicek et al., 2013) using three programs: MEME (Bailey et al., 2009), SCOPE (Carlson et al., 2007) and XXmotif (Luehr, Hartmann & Söding, 2012). We retained the enriched motifs that were present in both promoters and screened for potential transcription factor (TF) binding sites using STAMP (Mahony & Benos, 2007), with the motif position weight matrices (PWMs) retrieved from the TRANSFAC database (Matys et al., 2003).

Characterization of melanophore distribution in the fin regions

Comparisons of melanophore density and gene expression patterns were conducted between anatomically comparable regions of the dorsal, caudal and anal fins; i.e., either among the total areas cut from each fin (e.g. D-1 + D-2 + D-3 for the dorsal fin), or separately within each of the corresponding fin region classes; the distal fin regions (D-1, Cd-1, Cv-1 and A-1); the middle fin regions (D-2, Cd-2, Cv-2 and A-2); and the proximal fin regions (D-3, Cd-3, Cv-3 and A-3). Expression comparisons between different classes of fin regions were avoided because given the different histological properties of distinct anatomical fin regions along proximal-distal axis, gene expression differences between non-analogous regions could arise for various reasons not associated with colour patterning.

Melanophore density was significantly higher in the black distal regions of the dorsal fin and dorsal part of the caudal fin (D-1 and Cd-1) than in their dark grey ventral counterparts, A-1 and Cv-1 (Fig. 2C). The white dorsal middle regions (D-2 and Cd-2) contained almost no melanophores, whereas melanophore densities were intermediate in their grey ventral counterparts, A-2 and Cv-2, and finally the most proximal regions in all the fins had almost similar numbers of melanophores (Fig. 2C). Melanophore density clearly corresponded with the impression of darkness/lightness of the investigated fin regions. Additionally, the white colour of D-2 and Cd-2 regions appeared to be the result of both melanophore absence and an accumulation of white reflecting iridophores (Fig. 2B). Interestingly, total melanophore numbers summed across fin regions did not differ between fins (Fig. 2D). This indicates that the different fin patterns result from variation in the distribution and perhaps also pigmentation of a constant number of melanophores. In other words, aggregation of melanophores in the black regions and their absence in white regions, determines the dorsal stripes, whereas in their ventral counterpart regions the same number of melanophores is distributed more evenly.

Validation of reference genes

qPCR-based gene expression analyses depend on comparisons with stably expressed reference genes (Kubista et al., 2006), which have to be validated for the species and the specific experimental conditions in each study (Ahi et al., 2013). To this aim, we tested the expression of seven reference gene candidates on the cDNA generated from each of the 12 fin regions. The expression levels of the reference gene candidates varied from actb1, with the highest expression (lowest Cq) (Fig. S1), to ef1a and hprt1 with the lowest expressions (highest Cq). Next, the genes were ranked based on three algorithms, i.e., BestKeeper, geNorm and NormFinder, and standard deviation (SD) as described in Ahi et al. (2013) (Table 2). Among the reference genes, tbp ranked first by geNorm and NormFiner and second by BestKeeper analyses (Table 2). Hence, the data indicated high expression stability of tbp and suggested it as a suitable normalization factor to accurately quantify differences in gene expression between the fin regions.

Expression analyses of candidate genes

To assist the cellular characterization of the fin patterns, we investigated the expression levels of an iridophore lineage specific marker ltk (Lopes et al., 2008) and of the melanosome marker slc24a5 (Lamason et al., 2005) in the individual fin regions. The expression level of ltk was significantly higher in the white-coloured regions of the dorsal and caudal fins (D-2 and Cd-2) than in the corresponding ventral regions (A-2 and Cv-2), confirming that iridophores are accumulated in the white stripes. In contrast, expression of ltk in the distal and in the proximal fin regions was homogeneous across fins (Fig. 2E) (Table S3). The expression level of slc24a5 varied significantly among fins for each of the three fin region classes (Fig. 2E). Moreover, slc24a5 expression was significantly correlated with melanophore density across all fin regions (r = 0.89, p < 0.0001), although the increase of slc24a5 expression levels with melanophore number was stronger in the dorsal (D, Cd) than in the ventral (A, Cv) fin tissues (Fig. S2). Variation in the association between melanophore numbers and slc24a5 expression levels was also observed among fin regions, as the correlation was significant in the distal and the middle regions (r = 0.73, p = 0.001; r = 0.76, p = 0.0007, respectively) but not in the proximal regions of the fins (r = 0.22, p = 0.42). This suggests that variation in melanosome densities within the melanophores may also contribute to fin color contrasts.

Since the total number of melanophores, summed across regions, was constant in all fins, we next determined the average expression levels of each of the candidate target genes across the three regions of each fin and compared them among fins in order to identify spatial (particularly dorso-ventral) differentiation in gene expression. The selected candidate genes are known to be involved in stripe formation and/or adult pigmentation in zebrafish (see the details in Table 1). One of the genes, kir7.1 (kcnj13), had very low expression levels (Cq >  35) and therefore was discarded from rest of the analysis. This, however, suggests that kir7.1 expression is not required for chromatophores in fin regions of adult N. meeli, which is in agreement with findings in zebrafish,where the function of kir7.1 was not necessary for pigment cell survival in adults (Iwashita et al., 2006). Seven genes, igsf11, ltk, mitfa, mpv17, pmel, slc24a5 and slc45a2, showed differential expression between the fins (Fig. 3A). The expression level of the melanosome formation gene slc24a5 was significantly lower in the anal fin than in the other fins. Given that melanophore counts did not differ across fins (Fig. 2D), this indicates a reduced melanosome number in the melanophores of the anal fin. Furthermore, owing to the dorsal white stripes, expression levels of the iridophore marker ltk were significantly higher in the dorsal than in the ventral fin tissues. Dorso-ventral differences were also observed in the expression levels of igsf11, mitfa, pmel and slc45a2 (D, Cd > A, Cv), whereas mpv17 showed differential expression along the posterior-anterior axis (D, A > Cd, Cv).

Next, we were interested in those genes which showed overall dorso-ventral expression differences across the fins, and compared their expression within the distal, middle and proximal regions separately. Most interestingly, two genes, igsf11 and pmel, displayed higher expression in the black and in the white regions (D-1, Cd-1, and D-2, Cd-2) than in the corresponding ventral regions (Fig. 3B) (Table S3). Importantly, the elevated expression of these genes in the dark stripe was not simply due to the higher density of melanophores in these regions, as the differences persisted after correcting for the number of melanophores (last row in Fig. 3B). Elevated expression in both black and white stripes suggests that igsf11 and pmel have similar expression and potential function(s) in both iridophores and melanophores, which might emanate from their shared neural crest cell origin (Curran et al., 2010).

In contrast, the elevated expression levels of two further genes, mitfa and slc45a2, in the black stripe regions compared to the corresponding ventral regions (first row in Fig. 3B) levelled out after correction for differences in melanophore density among fin regions (last row in Fig. 3B). As discussed before, ltk expression was elevated only in the white stripes (Fig. 2E), consistent with its expression in  iridophores.

Finally, we extended our study by predicting potential TF biding sites in the upstream promoter sequences of igsf11 and pmel in Nile tilapia, an African cichlid with high quality annotated genome. Using different motif enrichment tools, we identified tens of motifs enriched in the promoter sequences of both genes. After parsing the motifs against the known TF binding sites in vertebrates, we compiled lists of top potential TFs binding to each motif (Table S2). We then analysed the expression of six of the TFs, which were predicted by all of the three motif enrichment tools and implicated in pigmentation processes (Barrallo-Gimeno et al., 2004; Cook & Sturm, 2008; Li et al., 2009; Agarwal et al., 2011; Besch & Berking, 2014; Natarajan et al., 2014) (Fig. 4). Five of these TFs showed slight but significant expression differences between the fins (Fig. 4A), however only one TF, irf1, showed dorso-ventral differences (D, Cd <  A, Cv). This pattern was confirmed in comparisons among anatomically equivalent regions, as the reduced expression of irf1 in the black stripe regions (first row in Fig. 4B) was robust against the correction for melanophore density (last row in Fig. 4B).

We also tested for expression correlations among the seven differentially expressed target genes (identified in Fig. 3A) and the six predicted TFs (Fig. 4C). The results showed positive expression correlations (blue shadings in Fig. 4C) between several pairs of target genes, including igsf11 and pmel, as well as negative expression correlations (red shadings) between target genes and three TFs, irf1, irf2 and ap2a. Notably, strong expression correlations with several target genes were observed for irf1, which was negatively correlated with igsf11, mitfa, pmel and slc45a2.