Identification of Gene Variants Associated with Melanocyte Stem Cell Differentiation in Mice Predisposed for Hair Graying

Age-related hair graying is caused by malfunction in the regenerative potential of the adult pigmentation system. The retention of hair color over the life of an organism is dependent on the ability of the melanocyte stem cells and their progeny to produce pigment each time a new hair grows. Age-related hair graying is variable in association with genetic background suggesting that quantitative trait loci influencing this trait can be identified. Identification of these quantitative trait loci may lead to the discovery of novel and interesting genes involved in stem cell biology and/or melanogenesis. With this in mind we developed previously a sensitized, mouse modifier screen and discovered that the DBA/1J background is particularly resistant to melanocyte stem cell differentiation in comparison to the C57BL/6J background. Melanocyte stem cell differentiation generally precedes hair graying and is observed in melanocyte stem cells with age. Using quantitative trait loci analysis, we have now identified three quantitative trait loci on mouse chromosomes 7, 13, and X that are associated with DBA/1J-mediated variability in melanocyte stem cell differentiation. Taking advantage of publicly-available mouse sequence and variant data, in silico protein prediction programs, and whole genome gene expression results we describe a short list of potential candidate genes that we anticipate to be involved in melanocyte stem cell biology in mice.

study that has successfully identified a genetic locus, IRF4, that is involved in age-related hair graying in humans (Adhikari et al. 2016). As an alternate to studying pigmentation in humans, the mouse model system has long been used to investigate genetic and molecular mechanisms related to melanocyte biology (Jackson 1994(Jackson , 1997Steingrímsson et al. 2006). Thus, we sought to assess whether mice, as an alternate to humans, could help in the identification of genetic variants that contribute the phenotypic diversity of hair graying. Previously we reported the development of a sensitized screen to evaluate genetically diverse inbred mouse strains for their ability to influence hair graying (Harris et al. 2015). In this screen we employed the Tg(DctSox10) transgenic mouse line to predispose mice to hair graying. The cellular mechanism responsible for hair graying in these mice is melanocyte stem cell (McSC) differentiation. This phenomenon precedes hair graying and is positively associated with hair graying severity (Harris et al. 2013). Differentiated McSCs have also been observed in human hairs and their number increases with age (Nishimura et al. 2005), which makes McSC differentation a relevant cellular phenotype to evaluate for genetic loci that may modify the extent of age-related hair graying in both mouse and humans. Mechanistically we predict that modifier genes that effect McSC differentiation may directly regulate the process of melanogenesis, but could also be involved in initial McSC establishment, proliferation or migration of McSC progenitors.
McSCs that are undergoing differentiation produce visible ectopic pigmentation when viewed by light microscopy and the number of hairs that contain these ectopically pigmented McSCs (EPMs) varies in animals of different genetic backgrounds (Harris et al. 2015). Mice that are hemizygous for a conditional, Sox10-expressing transgene, designated C57BL/6J-Tg(DctSox10)/0, are extremely susceptible to McSC differentiation. In contrast, progeny derived from mating C57BL/6J-Tg(DctSox10)/0 mice to other inbred genetic backgrounds (C3H/HeJ, 129SvEvTac, FVB/NTac, DBA/1J, BALB/CJ) produced F1 hybrids that all exhibit reduced numbers of hairs with EPMs in response to the transgene (Harris et al. 2015). Tg(DctSox10)/0 hybrid animals produced by mating C57BL/ 6J-Tg(DctSox10)/0 mice to DBA/1J mice exhibit a particularly low level of transgene-mediated EPMs. Reduction of EPMs in these F1 hybrids suggests a dominant mode of inheritance for EPM resistance and we sought to identify these DBA/1J-associated resistance loci using QTL linkage mapping.

Phenotype Analysis
Assessment of hairs with EPMs was performed as described previously (Harris et al. 2015). Briefly, between 9-11 weeks of age, hair along a 2x2 cm region of the lower back was plucked by hand to synchronize and initiate hair regrowth. Hairs in the plucked region were allowed to regrow for seven days (equivalent to hair cycle stage anagen III/IV). Skin from these animals was dissected and processed for cryosectioning. Using light microscopy, approximately four 10 mm sections were analyzed in total skipping at least three sections between those analyzed to prevent counting the same hair twice. EPMs occur in the hair bulge at the insertion point of the arrector pili muscle, thus only sectioned hair follicles that spanned the entire region from the sebaceous gland and past the junction of the dermal/subcutis were counted. Between 50-100 hairs were examined for the presence or absence of EPMs. The percentage of hairs with EPMs was calculated by dividing the number of hairs with EPMs by the total number of hairs analyzed. 122 N2 mice were phenotyped. For QTL analysis, selective genotyping was performed on 79 of the N2 mice exhibiting the highest (. 50%, n = 39) and lowest (,20%, n = 40) percentage of hairs exhibiting EPMs. The EPM phenotype was converted to a binary trait with the high EPM phenotype scored as 1 and referred to as 'affected', while the low EPM phenotype was scored as 0 and referred to as 'unaffected'. Graphing was performed using Graphpad Prism (Graphpad Software). Brightfield microscopy was performed on a Zeiss Observer.D1 compound microscope. Images were obtained with an Axiocam Hrc camera (Zeiss) using the ZEN software (Zeiss) and processed with Adobe Photoshop (Adobe).

Genotyping
Presence of the Tg(Dct-Sox10) transgene was determined by PCR using primers that generate an amplicon spanning the Dct promoter and the Sox10 cDNA: 59-AGCAGTATGGCTGGAGCACT-39; 59-TCCAGTCGTAGCCGCTGAGCA-39. PCR cycling was performed as published previously (Harris et al. 2013). SNP genotyping was performed on a custom panel of 1449 SNPs (equivalent to the Mouse Medium Density Linkage Panel, Illumina) using the GoldenGate Genotyping Universal-32 Assay Kit with UDG (Illumina). Complete SNP genotyping data are available in Supplemental File 1 (sheet name-Original Sample Genotypes). In preparation for QTL analysis using R/qtl, 559 non-informative SNPs and SNPs with a high number of no call (NC) values or were omitted. Sample genotypes were also recoded such that homozygote genotypes matching the parental C57BL/6J genotype were designated AA, and heterozygote genotypes were designated AB. SNPs with unknown chromosomal coordinates in the original genotyping data (listed as chr 0) were identified in the genome. The chromosomal coordinates of each SNP were then converted to Sex-Averaged cM-G2F1 centimorgan positions using Mouse Map Converter (http://cgd.jax.org/mousemapconverter/). The final data matrix used for QTL analysis is included in Supplemental File 1 (sheet name-Converted Sample Genotypes).

Statistical analysis
A total of 79 mice were initially evaluated, and 3 removed for low-quality genotyping results. QTL linkage analysis of 76 mice (36 with the high and 40 with the low EPM phenotype, 37 males and 39 females) and 890 markers using EPM as a binary trait was performed using the R/qtl software (Rv.3.4.4,. One-dimensional scans were conducted without and with sex as an additive covariate using logistic regression with the EM algorithm (Dempster et al. 1976;Lander and Botstein 1989;Lander and Botstein 1989). Separate LOD significance thresholds were obtained from 1,000 permutations for autosomal SNPs and 18,850 permutations for X chromosome SNPs (see R/QTL documentation for explanation of the permutation counts applied within the scanone function). Two-dimensional scans were conducted with sex as an additive covariate using logistic regression with the EM algorithm. Separate LOD significance thresholds were obtained from 1,000 permutations for autosomal SNP pairs, 355,289 permutations for X chromosome SNP pairs, and 9,425 permutations for autosomal:X chromosome SNP pairs (see R/QTL documentation for explanation of the permutation counts applied within the scantwo function). For the QTL found significant at the 0.05 level in the two-QTL analyses, multiple-QTL analyses were performed with sex as an additive covariate using logistic regression with multiple imputation (Sen and Churchill 2001). The locations of the QTL were updated based on maximum likelihood (Zeng et al. 1999) using the refineqtl function within R/QTL.

Identification of candidate genes
Whole skin RNAseq data comparing C57BL/6J to DBA/1J was retrieved at NCBI GEO using the accession # GSE86315. This dataset included read counts previously generated by aligning RNAseq reads to the mouse genome (GRCm38/mm10) using TopHat2, assessment of mapping quality using RSeQC and RNA-SeQC, and read counting using HTSeq (Swindell et al. 2017). HTSeq reads from the C57BL/6J and DBA/1J control skin samples (2 males and 2 females per strain all treated with a non-toxic lanolin-derived occlusion cream) were used to generate normalized read counts (median ratio method) and compared to obtain differential expression values using DESeq2. The complete differential expression data including base mean (mean of the normalized counts), log2 fold change, and adjusted p values is provided in Supplemental File 2. Wildtype C57BL/6J McSC RNAseq data were retrieved at NCBI GEO using the accession # GSE102271. RNAseq reads were aligned to the Ensembl GRCm38.p5 primary DNA assembly using STAR (v2.5.2b) and normalized read counts (median ratio method) determined using DESeq2. Normalized read counts with a value of 0 were omitted from further analysis. The complete McSC expression data are provided in Supplemental File 3. Variants between C57BL/6J and DBA/ 1J were obtained from the Mouse Genome Project (REL-1505; ftp://ftpmouse.sanger.ac.uk/). The 'genome variants' function of PROVEAN (v1.1.3 and GRCm38 Ensembl 74; http://provean.jcvi.org) was used to score variants based on predicted protein function. PROVEAN summary and detailed results are provided in Supplemental File 4.

Reagent and Data Availability
All data associated with this manuscript is available within the manuscript or as supplemental files. Supplemental Files 1-4 are available via the GSA Figshare portal. Supplemental File 1 provides the SNP genotyping data of the N2 animals. Supplemental File 2 provides the differential mRNA expression data comparing DBA/1J and C57BL/6J skin from the reanalysis of GSE86315. Supplemental File 3 provides the mRNA expression data of wildtype C57BL/6J McSCs from the reanalysis of GSE102271. Supplemental File 4 provides the summary and detailed output from PROVEAN. Supplemental material available at Figshare: https://doi.org/10.25387/g3.7489625.

RESULTS
Distribution of EPM susceptibility in progeny derived from backcrossing D1B6F1-Tg(DctSox10)/ 0 to C57BL/6J To identify the genetic determinants from DBA/1J that promote EPM resistance we set out to map loci that modify the production of EPMs in progeny produced from backcrossing D1B6F1-Tg(DctSox10)/0 females to C57BL/6J males (these progeny are hereafter referred to as N2 mice). At approximately eight weeks of age, the hair was plucked along the lower back of N2 mice to induce and synchronize hair growth. One week later, skin from these mice was obtained from the plucked region and assessed for EPMs using histological methods ( Figure 1a). 122 N2 mice that carry the Tg(DctSox10) transgene were evaluated phenotypically, and exhibited a range of EPM measurements extending between the C57BL/6J-Tg(DctSox10)/0 and D1B6F1-Tg(DctSox10)/0 parental phenotypes (Figure 1b). A statistically significant gender effect (t-test, p-value = 0.009) is also observed in N2 mice with the phenotypic mean of the female N2 animals skewed toward more resistant to EPMs suggesting the need for including sex as covariate during QTL mapping.
QTL analysis provides evidence for three QTL loci associated EPM variability To map genetic modifiers that affect resistance to McSC differentiation, we used a selective genotyping approach where only the animals with the most extreme phenotypes are genotyped for linkage analysis (based on (Lander and Botstein 1989)). A total of 79 N2 mice were genotyped and represent those animals with the highest (. 50%, n = 39) and lowest (,20%, n = 40) percentage of hairs exhibiting EPMs (Figure 1b). These  Results from the single-QTL analysis. Cells shaded in gray highlight loci described in the Results. chr, chromosome; pos, centimorgan position; lod, LOD value; pval, p-value.
mice were genotyped using a panel of 1449, evenly-dispersed, mousespecific SNP loci assays (Illumina, Supplemental File 1). Among the 1449 SNPs evaluated, 890 were found to be reliable and informative between C57BL/6J and DBA/1J. 3 of the original 79 mice had lowquality genotyping scores and were removed prior to QTL analysis. Identification of individual QTL was first performed using a single-QTL genome scan approach (R/qtl; Broman et al. 2003). Only the genotyped animals mentioned above were included in these analyses and thus the high and low EPM percentage were treated as a binary trait. Results from the single-QTL analysis, without sex as a covariate and a 5% significance threshold, indicates the presence of one QTL on chr 13 (lod-3.56, p-value-0.02; Table 1). Using a 10% significance threshold, there is also weak support for an additional QTL on chr 7 when sex is included as an additive covariate (lod.add-2.91, p-value-0.09; Table 1). Seeing as there is a sex-dependent difference in the phenotypic mean of the N2 animals ( Figure 1b) we also evaluated for an interaction between QTL and sex and find that the LOD score for the chr 7 QTL increases if an interaction is allowed (lod.full-3.81, p-value 0.07; Table 1). Stratifying the data by analyzing males and females separately suggests that the chr 7 QTL is female-specific (lod.male-0.39, p-value-1.00; lod.female-3.54, p-value-0.03; Table 2), however, the interaction between sex and the chr 7 QTL is not significant (lod.int-0.90, p-value-0.40). The effects of each QTL were visualized by plotting the proportion affected as a function of genotype at the SNP markers nearest the chromosomal location with the highest LOD score (Figure 2). The effects at the chr 7 QTL are consistent with the DBA/1J allele conferring resistance to McSC differentiation with a low proportion of animals exhibiting the affected, high EPM phenotype in association with heterozygosity (AB) for the C57BL/6J and DBA/1J alleles (Figure 2a). In addition, females exhibit a noticeably larger effect than males, matching the above evidence suggesting that the chr 7 QTL is influenced by sex. The chr 13 QTL, on the other hand, exhibits an effect that suggests that this QTL may be involved in DBA/1J-mediated susceptibility to McSC differentiation. In this case, a high proportion of animals exhibit the affected, high EPM phenotype in association with heterozygosity for the C57BL/6J and DBA/1J alleles (Figure 2b). Single-QTL analysis for the chr X was also performed and no loci were identified that met the 10% significance threshold (results not shown).
In order to look for additional QTL that may contribute to DBA/1Jmediated EPM resistance, a two-QTL genome scan approach was performed. Results from the two-QTL analyses focusing on pairs of autosomal QTL (with or without sex as an additive covariate) does not identify any pairs of that significantly improve the two-QTL model (lod. full or lod.add) above that of a single QTL (lod.fv1 or lod.av1) when using a 5% significance threshold (Table 3). However, when assessing for pairs of QTL between the autosomes and the X chromosome (with sex as an additive covariate), there is significant evidence for a pair on chr 7 and the X chromosome (Table 3). The full model (lod.full) containing QTL at chr 7 and the X chromosome provides a better fit to the data than both the best single-QTL model (lod.fv1) and the additive model (lod. add). Interaction between these QTL is also significant (lod.int) suggesting that these QTL are epistasic. When the chr 7 and the X chromosome genotypes are considered together (Figure 2c), the effect of the chr 7 QTL in females is the same; a low proportion of individuals exhibit the affected, high EPM phenotype when heterozygous for the chr 7 QTL, independent of genotype at the X chromosome QTL. In males, however, the effect of the chr 7 QTL is opposite depending on whether the X chromosome is C57BL/6J-derived (A allele) or DBA/1Jderived (B allele). Specifically, the effect of the chr 7 QTL switches from EPM resistance in combination with the X chromosome B allele to EPM susceptibility in males with a the X chromosome A allele. No pairs of QTL on the X chromosome approach the criteria for significance (results not shown).
In summary, single-QTL and two-QTL analysis of the genotype and phenotype data of N2 backcross animals identifies three QTL loci that contribute to the EPM phenotype. These QTL reside on chr 7, 13 and X. The chr 7 QTL is associated with resistance to EPMs, but this effect can be modified in males by an additional epistatic QTL located on the X chromosome. The interaction of these two QTL help explain the sexspecific difference observed in the N2 distribution (Figure 1a). The chr n 13 QTL, on the other hand, has an effect opposite to that of the QTL at chr 7 and is associated with EPM susceptibility.
Identification of candidate genes that influence the EPM phenotype As the first step to prioritizing candidate genes for the QTL identified above, we used 1.5-LOD support intervals (the chromosomal region where the LOD score is within 1.5 of its maximum) to determine the bounds of the QTL linkage region (Figure 3, Table 4). In order to consider the chr 7 and chr X QTL simultaneously we performed multiple-QTL analysis to fit these QTL into one model, refine the QTL locations, and used these LOD estimates to define the 1.5-LOD support intervals (Figure 3a). The chr 13 intervals were derived from the single-QTL analysis (Figure 3b). For these QTL, the 1.5-LOD support intervals were relatively large and encompassed a significant number of genes; the chr 7 QTL interval covered 22 Mbp with 444 genes, the chr 13 QTL covered 58 Mbp with 603 genes, and the X chromosome QTL interval covered 44 Mbps with 377 genes.
As an approach to stratify the genes represented in these intervals we characterized them using publicly-available gene expression data, mouse resequencing data, an online program for predicting the effects of coding variants on protein function, and a curated list of genes known to be involved in pigmentation. First, any genes that are differentially expressed in skin from DBA/1J and C57BL/6J animals highlights cell-specific or systemic changes that may impact McSC function. For instance, genes that are upregulated in DBA/1J animals may be dominant drivers of EPM resistance that promote McSC stemness over differentiation, while those that are downregulated in DBA/1J may highlight proteins responsible for heightened melanogenesis in C57BL/6J animals. Using NCBI Gene Expression Omnibus (www.ncbi.nlm.nih.gov/geo/) we identified a study that compared the transcriptomes of skin derived from 11-week-old DBA/1J and C57BL/6J animals using RNAseq. This study evaluated strain-specific effects of Imiquimod on skin, and included expression data from control subjects (two males and two females for each strain) that were treated only with a non-toxic lanolin-derived occlusion cream (Swindell et al. 2017). Using DESeq2 to compare these DBA/1J and C57BL/6J control subjects we find a number of genes that exhibit a significant twofold difference in expression (p-adjusted , 0.05; Supplemental File 2).
Second, we were also interested to capture genes that are expressed at detectable levels in McSCs that are not up or downregulated in DBA/1J and C57BL/6J skins. We anticipate these genes represent candidates that may alter intrinsic McSC dynamics at the protein level, either through differential transcript expression or as the consequence of genetic variation in their coding sequence. Using Third, we evaluated the genes within the linkage intervals for those that exhibit genetic variation between C57BL/6J and DBA/1J as well as those that are known to participate in pigmentation. Mouse resequencing and associated variant call data were used to identify genes that contain coding variants between C57BL/6J and DBA/1J (Mouse Genomes Project, Sanger; GRCm38) (Yalcin et al. 2011;Keane et al. 2011). These include missense, insertion and deletion mutations that are predicted to have deleterious protein consequences by PROVEAN (Choi and Chan 2015) as well as nonsense and frameshift mutations (Supplemental File 4). Genes known to be involved in pigmentation were derived from a comprehensive and curated gene list (Baxter et al. 2018) Using the four criteria described above, we developed a short list of potential candidate genes. This short list includes any genes that exhibit a statistically-significant, twofold difference in gene expression between DBA/1J and C57BL/6J skins, any genes with a percent rank of expression greater than 50% in C57BL/6J McSCs that also contain a deleterious coding mutation, and any known pigmentation genes. While this approach may selectively filter out some candidate genes that are expressed nonautonomous to the McSC and those genes that are under the influence of differential transcript regulation or post-translational modifications, this abridged list provides a reasonable starting point for future follow-up (Table 5).

DISCUSSION
Modifier genes are recognized for their contribution to phenotypic variation observed in disease and it's been suggested that their identification may help predict "sensitive pathways and nodes for therapeutic intervention (Hamilton and Yu 2012)." With this in mind, we were interested in identifying naturally-occurring genetic variants in mouse that could contribute to variability in the phenotype of hair graying. Hair graying is the consequence of disrupting the regenerative activity of McSCs within the hair follicle or the function of their progeny. We developed a sensitized screen using mouse inbred lines as an unbiased approach to identify novel genes that participate in these mechanisms. Using the Tg(DctSox10) transgene as a condition to predispose mice to n lod.fv1-difference between the full LOD and the max single-QTL LOD for the chromosome pair; lod.add-max LOD score for the additive model for the chromosome pair; lod.av1-difference between the additive LOD and the max single-QTL LOD for the chromosome pair; lod.int-difference between the max full LOD and the max additive LOD for the chromosome pair (taken from r/qtl package documentation; see Broman et al., 2003 for additional details). Cells shaded in gray highlight models described in the Results. chr, chromosome; pos, centimorgan position.
hair graying via the mechanism of McSC differentiation, we found that the DBA/1J genetic background provides resistance to this cellular phenotype (Harris et al. 2015). In search of modifier genes that could mitigate McSC differentiation we performed QTL analysis on progeny derived from backcrossing (C57BL/6J x DBA/1J)F1-Tg(DctSox10)/ 0 animals to C57BL/6J. In brief, we identified three linkage regions across three chromosomes: chr 7, 13 and the X chromosome. While we were particularly interested in identifying QTL that could explain reduced McSC differentiation (the low EPM phenotype) associated with heterozygosity for DBA/1J, the chr 7 and 13 QTL have opposing effects consistent with promoting EPM resistance and EPM susceptibility, respectively (Figure 2a, b). Two-QTL analysis also revealed a novel, epistatic interaction between the QTL on chr 7 and the X chromosome. The effect of the X chromosome QTL is only observed in males and functions to toggle the effects of the chr 7 QTL from one of EPM resistance to EPM susceptibly when the allele at the X chromosome QTL is C57BL/6J-derived (Figure 2c). This unique interaction helps to explain the reduced effects of the DBA/1J allele at the chr 7 QTL in males (Figure 2a), as well as sex-specific skewing of the EPM phenotype in the N2 population as a whole (Figure 1).
Likely due to the relatively small number of N2 progeny evaluated for linkage, the 1.5-LOD intervals for these QTL were large and encompassed a number of genes. These were selectively filtered to include only those with differential gene expression in whole skin between C57BL/6J and DBA/1J, those with evidence to support expression within McSCs but with altered protein function, and genes associated with pigmentation phenotypes (Table 5). In considering the function of these candidate genes, we anticipated that all three of the QTL could influence McSC behavior at any number of timepoints, and include both autonomous and non-autonomous mechanisms. In general, McSCs reside in a specialized niche within the hair follicle and are activated in coordination with hair follicle stem cells. At each new hair cycle McSCs proliferate and produce progeny that will colonize the hair bulb. These progeny cells will differentiate into melanocytes and initiate the process of melanogenesis, which includes the synthesis of melanin within melanosome and the trafficking of these melanosomes for deposition into keratinocytes of the growing hair shaft (Osawa 2009). EPMs, like those in Tg(DctSox10) mice, are generated when McSCs within the stem cell niche do not self-renew properly and instead differentiate prematurely (Nishimura et al. 2005;Harris et al. 2013). Thus, resistance or susceptibility to EPMs could be the result of mitigating or exacerbating this process as early as when McSCs are deciding their fate or as late as the final steps of pigment production. In addition, because the initial EPM phenotype is driven by the Tg(DctSox10) transgene as a consequence of Sox10 overexpression, these QTL could also act as regulators of SOX10.
As one example, there are a number of genes within the linkage intervals for chr 7, 13 and the X chromosome that have known roles in cellular mechanisms associated with the function of the melanosome organelle. These include Oca2, Trpm1, Bloc1s5, Dtnbp1, Kif13a, and Atp7a. Since our phenotypic assessment of EPMs is dependent on the production of visual pigmentation, variability in the EPM phenotype may reflect changes in melanosome biogenesis or maturation. In particular, Trpm1 encodes a protein called transient receptor potential cation channel, subfamily M, member 1 (also known as melastatin) and is highly expressed in C57BL/6J McSCs (Table 5). TRPM1 localizes to non-melanosomal vesicles and its activity is associated with increased intracellular melanin content (Oancea et al. 2009). Analysis of mouse resequencing data demonstrates that the Trpm1 gene contains DBA/1J-related missense and frameshift coding mutations that n are predicted to alter protein function (Table 5). Thus it follows that the low EPM phenotype associated with the DBA/1J allele at the chr 7 QTL could be a consequence of malfunctional TRPM1 protein, and a subsequent reduction in pigment production in Tg(DctSox10) McSCs. Interestingly, Trpm1 knockout in mouse produces no apparent congenital pigmentation defect on its own (Koike et al. 2010) but is consistent with our hypothesis that, as a modifier gene, it may only act to modify previously existing coat color phenotypes in mouse. As a second example, quite a few of the candidates identified are involved in transcriptional regulation and could play a role in shaping McSC self-renewal or differentiation. These include Irx1, Irx2, Mef2c, Rsl1, Tfap2a, Zfp369, Cited1 and Med12. Mef2c and Tfap2a are contained within the chr 13 QTL interval, are expressed at relatively high levels in C57BL/6J McSCs (Table 5), and both have known function in promoting melanocyte differentiation (Agarwal et al. 2011;Seberg et al. 2017). The DBA/1J allele at the chr 13 QTL is associated with EPM susceptibility making Mef2c and Tfap2a plausible candidates for enhancing premature differentiation of Tg(DctSox10) McSCs.
As this study highlights, one benefit to using the natural variation of inbred mouse lines to evaluate QTL is that many of the lines have now been sequenced ). In addition, the structural variation between these lines and the reference genome is known and searchable (Yalcin et al. 2011). While additional fine mapping or the use of chromosome substitution mouse lines are helpful to shrink large linkage intervals, genomic tools, like computational prediction of protein function and whole genome gene expression data are now more readily available and can be applied to the identification of promising candidates. This is particularly helpful when limited resources prohibit the generation of additional crosses to refine the linkage region. With this approach we identified a list of gene candidates that may contribute to DBA/1J-specific changes in McSC differentiation and each will have to be tested individually for causality. DBA/1J and DBA/2J share sufficient homology that BAC transgenesis using available DBA/2J BAC clones (MM_DBa) may be a good way to begin confirming these candidates. For those genes that contain missense mutations of varying predicted consequence, assessing each for its effects on protein secondary structure in silico (e.g., using the RCSB Protein Data Bank), and its effects on stability and subcellular localization in vitro can help in the selection of those that are more likely to be detrimental to cell function in vivo. In addition, EPM phenotyping has been performed on other genetic backgrounds (Harris et al. 2015), and any variants that are common between strains that exhibit the low EPM phenotype and distinct between the high EPM strains could provide additional support for specific candidates. Finally, with the versatility of CRISPR/Cas9-mediated genome editing, creating mice with specific candidate gene defects will help to finalize the identification of genes involved in the variation of McSC differentiation, a cellular pathology associated with hair graying.

ACKNOWLEDGMENTS
We thank Arturo Incao for indispensable help with mouse colony maintenance and the remaining members of the Pavan Lab for their general support. In addition, we recognize Marypat Jones and Ursula Harper in the NHGRI Genomics Core for assistance with SNP genotyping. Candidate gene list for the QTL on chr 7, 13, and the X chromosome. The DBA/B6 columns refer to differential expression from RNAseq data comparing C57BL/6J and DBA/1J whole skin. Those values in bold represent genes that have both an adjusted p-value of ,0.05 and a twofold change in expression in either direction. The B6 McSC columns refer to the ranked RNAseq data from C57BL/6J McSCs. Deleterious mutations were determined using PROVEAN, and the type of mutation is provided. Known pigment genes are marked with 'yes'. ENS gene ID, Ensembl gene ID; TX START, bp position of the transcription start site (GRCm38); BASE MEAN, mean of the normalized RNAseq counts; LOG2FC, log2 fold change; PADJ, adjusted p-value; %RANK, rank of expression for the indicated gene as a percentage of the entire genes within the dataset.