RAD-tag sequencing
A total of 743,486,491 reads were produced from the 5 lanes of Illumina HiSeq for the 160 F2 progeny. 662,570,635 (89.1%) passed our Q20 filter. A total of 624,492,015 reads (94.3% of the filtered reads) were successfully processed for barcodes by Stacks. This corresponds to approximately 7.7 million reads per F2 progeny. After filtering and barcode processing, there were 7,535,127 reads for the F0 female and 18,850,602 reads for the F0 male used for the cross.
Map construction and anchoring
We scored 834 genetic markers in 160 male F2 progeny from the M. zebra x M. mbenjii cross. The average coverage of each genotype SNP was 49.9x (range of 6.9x-201.7x) in the F2 progeny, 254.9x in the male F0, and 105x in the female F0. The average genotype completeness was 77% and the frequencies of each genotype class were 27.7% AA, 45.5% AB, and 26.9% BB (A designated for grandfather alleles and B for grandmother alleles). 60 individuals had between 0–100 genotypes missing, 42 individuals had 101–200 genotypes missing, 25 individuals had between 201–300 genotypes missing, 22 individuals had between 301–400 genotypes missing, 7 individuals had between 401–550 genotypes missing, and 4 individuals had greater than 500 genotypes missing. The average number of missing genotypes per individual was 182. High numbers of missing genotypes can be attributed to low coverage in some of the F2 individuals, with coverage ranging from 357,000 reads to 10,500,000 reads. We obtained a linkage map that contained 22 linkage groups and spanned over 1,933 cM. This agrees with previous work indicating that there are 22 chromosomes in cichlids [25]. Marker density was approximately one marker per 2.5 cM. Additional file 1 contains the information used to create the linkage map.
This linkage map was then used to order scaffolds of the M. zebra genome assembly. To be included in the anchored map, we required that scaffolds be anchored by at least two markers in the linkage map. We found 114 scaffolds (6.5 per linkage group) that met this criterion. The average size of these scaffolds was 3,918,467 bp for a total of 564,259,264 anchored bp. This represents 66.5% of the 848,776,495 bp M. zebra genome assembly. An additional 110 scaffolds had a single hit to a particular linkage group and were not included in the anchoring. If these single scaffolds were included in the anchoring, 92.3% of the assembly would become anchored. Additional file 3 provides the placement of scaffolds on the linkage groups.
QTL scan and detection
A genome wide scan resulted in the identification of three QTL. Phenotypes with significant LOD scores included dorsal fin xanthophores, caudal fin xanthophores, and pelvic fin melanophores (Figure 2). The failure to detect QTL for the remaining pigmentation traits is most likely due to the small difference in the parental means for the other pigmentation traits. QTL were detected for the three traits for which parental means were 2.8 or greater standard deviations apart. The mapping population of 160 male F2 did not have enough power to detect QTL for traits for which the difference in parental means was smaller than 2.8 standard deviations.
Dorsal fin xanthophores
One QTL with a LOD score of 11.15 was detected for dorsal fin xanthophores on LG12 and explained 27.4% of the variance (Figure 3A). The Bayesian credible interval for this QTL spans from 62–73 cM along LG12. This matches the estimate from our previous work that a minimum of one QTL would be detected for this trait [19]. The effect plot for the marker with the highest LOD score shows incomplete dominance. Individuals with more orange on their dorsal fin possess at least one allele from the M. mbenjii grandfather (Figure 4A). The results of the effect plot are consistent with the inheritance of orange fins from the M. mbenjii grandfather [19].
Using the annotated genome assembly for M. zebra, we were able to identify candidate genes within the credible interval. Little work on the genetics of xanthophore development and carotenoid formation has been done in fishes. Csf1r, the most obvious candidate gene for xanthophores traits [8], is not present in the interval. Several other genes that might be involved in pigment cell development are present, including genes involved in vesicle formation, carotenoid synthesis, and cell aggregation (Additional file 4: Table S2). TRPM 6/1, a member of a gene family previously linked to pigment cell development in zebrafish, is located in the interval. However, this gene has only been demonstrated to be important in melanophore development. Fish possessing a mutation in this gene experience melanophore death, while the other pigment cell lineages (xanthophores and iridophores) appear to develop normally [28, 29]. We do not consider TPRM6/1 a primary candidate gene for dorsal fin color since the difference between the grandparent fins appears to be a trade-off between the presence of xanthophores versus iridophores.
Caudal fin xanthophores
One QTL with a LOD score of 5.21 was identified for caudal fin xanthophores on LG 12 in a region that overlaps with, but is broader than, the QTL region identified for dorsal fin xanthophores. The Bayesian credible interval spans from 56–81 cM along LG 12 (Figure 3B). Identification of a shared QTL is not surprising since our previous work indicated a strong correlation between these traits [19]. The QTL plot for caudal fin xanthophores is broader, possibly because this trait shows less variance in the caudal fin than it does in the dorsal fin. While the xanthophores in the caudal fin of M. mbenjji are found in stellate form and overlap to form a uniform field, M. zebra often possess small punctate xanthophores on their caudal fins.
Although the QTL region for dorsal fin and caudal fin xanthophores is shared, the percent variance explained for each trait is different, with 27.4% of the variance explained for the dorsal fin xanthophores, but only 14% of the variance explained for the caudal fin xanthophores. This difference could be due to the number of genes predicted to control each of these traits. While the Castle-Wright estimator predicted that one gene controls dorsal xanthophores, a minimum of three genes was predicted to control caudal fin xanthophores [19]. We were probably only able to detect one gene for caudal fin xanthophores due to the power limitations discussed previously.
Similar to the effect plot for dorsal xanthophores, individuals with more orange on their caudal fin possess at least one of the M. mbenjii grandfather’s alleles (Figure 4B). It should be noted that while the trait appears to show overdominance, the mean for the individuals homozygous for the M. mbenjii alleles and the mean for heterozygotes are not significantly different. Candidate genes considered for this region are the same as those considered for the dorsal fin xanthophores.
Pelvic fin melanophores
One QTL with a LOD score of 5.01 was identified for pelvic fin melanophores, with a Bayesian credible interval spanning from 68–82 cM on LG 11 (Figure 3C). This QTL explains 13.4% of the variance for this trait. We previously estimated that this trait was controlled by a minimum of three genes. The relatively low power of our mapping population may explain why only one QTL was found for this trait. The jagged QTL curve for this trait may be due to a combination of high marker density and occasional missing genotypes.
Effect plots show that individuals with more melanophores on their pelvic fin are homozygous for the M. zebra’s grandmaternal alleles (Figure 4C). This is consistent with the fact that M. zebra males possess more melanophores on their pelvic fin compared to M. mbenjii males [19]. No obvious candidate genes were identified, however various genes in the interval play a role in human skin disease, interact with other genes involved in melanophore development, or are involved in vesicle formation, packaging, and trafficking (Additional file 4: Table S2).