A gene-based recessive diplotype exome scan discovers FGF6, a novel hepcidin-regulating iron-metabolism gene

Genome-wide association studies (GWASs) are well designed to detect additive effects of modest effect sizes. We hypothesized that gene-based tests sensitive to recessive diplotypes, including recessive single-site effects and compound heterozygosity, may reveal additional genes underlying complex diseases. Carrying variants conferring a compromised function on both homologous chromosomes is likely to impact molecular physiological states. Deep-sequencing studies have conclusively shown a vast reservoir of rare variants segregating in human populations.¹  Rare variants in functional categories (eg, missense, regulatory motifs) may generate pathogenic effects through recessively acting diplotypes, and such effects are apt to remain concealed from standard GWAS analyses. Simple power calculations show that recessive diplotype inheritance produces signals that are difficult for standard GWAS methods to discover (supplemental Figure 1, available on the Blood Web site). Compound heterozygosity disease models also enjoy a high degree of biological plausibility, particularly if the alleles confer compromised protein function.^2-5  Recessive diplotype modes of inheritance are well established in Mendelian diseases, such as cystic fibrosis,⁶  mevalonic aciduria,⁷  β-thalassemia,⁸  and Niemann-Pick disease.⁹  Although not systematically examined in population-based studies, there is a sizable repository of genes underlying complex diseases with recessive, loss-of-function effects.^10-14  Hence, we posited that an exome-wide, gene-based screen of recessive diplotypes using putative functional variants in both oligogenic and complex diseases may expand our knowledge of disease genes.

Iron-metabolism disorders, including adult hereditary hemochromatosis, collectively are common conditions with considerable public health implications.^15,16  Importantly, the hepatic hormone hepcidin is a key regulator of iron homeostasis by controlling iron flux from enterocytes and macrophages to plasma through degradation of the cellular iron exporter ferroportin (SLC40A1). Within cells, ferritin is the iron-storage protein that can be used for indirect iron quantification. To investigate the inheritance of hemochromatosis, several segregation analyses were initially conducted, concluding that a recessive mode of inheritance is highly plausible.^17,18  Several human studies have investigated the genetics of iron overload, revealing several critically important genes. Notably, HFE, encoding the membrane-bound hereditary hemochromatosis protein, was mapped 2 decades ago through family-based linkage^19-22  and association approaches.^23,24  Additional studies have definitively placed the missense polymorphism C282Y (rs1800562) in HFE as the major susceptibility factor in adult-onset, type 1 hereditary hemochromatosis.^25,26  Additional genes have been identified through pathway-based genetic association studies and GWASs, including BMP2, BMP4, HJV, TF, TMPRSS6, NAT2, FADS2, and TFR2.^27-29 

The homogenous population in rural central Wisconsin is the source population for the Personalized Medicine Research Project (PMRP), a biobank linked to electronic health records (EHRs) housed by the Marshfield Clinic Research Institute.³⁰  Samples from over 20 000 individuals comprise the PMRP. The study was conducted in accordance with the Declaration of Helsinki. All samples were collected following written informed consent. All investigators using the PMRP samples had obtained Research Ethics and Compliance Training certification through the Collaborative Institutional Training Initiative (CITI) program. The study protocol was reviewed and approved by the Marshfield Clinic Institutional Review Board (details in “Acknowledgments”). The Central Wisconsin population is largely stationary and primarily derived from Bavarian migrants in the late 1800s. The population carries high utility for disease gene mapping through reduction in confounding by population stratification and lower expected levels of allelic and locus heterogeneity. Additionally, environmental exposures are thought to be relatively uniform across this population. For these reasons, the PMRP has been effectively used in numerous human genetics studies.^31-34  PMRP DNA samples were collected and stored ∼14 years ago and all individuals have longitudinal EHR information housed at the Marshfield Clinic, averaging >30 years. The EHR is composed of International Classification of Diseases, Ninth Revision (ICD-9) diagnostic codes, laboratory test results, clinical procedure data, prescription information, and physician notes. Hemochromatosis cases and controls were selected from the PMRP population. Hemochromatosis cases were selected on the basis of the percent transferrin saturation laboratory values (the ratio of serum iron to transferrin iron-binding capacity) exceeding 48% and having 2 or more instances of ICD-9 codes indicating the diagnosis of hemochromatosis: 275.0 (iron metabolism disorder, excluding anemia), 275.01 (hereditary hemochromatosis), 275.03 (unspecified hemochromatosis), and/or 275.09 (other iron-metabolism disorders). To reduce confounding by population stratification, a principal components analysis on the exome-genotyping data was implemented using all samples, blinded to disease status. Individuals considered genetic background outliers (>3 standard deviations from the centroid of the first 2 principal components) were excluded from the study. Following the removal of outliers, the resulting set of individuals was highly homogeneous based on the first 3 principal components. Exhaustive pairwise kinship coefficients were calculated and 1 individual from pairs of individuals exhibiting third-degree or closer relatedness were removed. Of the ∼10 000 individuals previously subjected to the exome-genotyping array and quality control procedures, the phenotype algorithm identified 18 individuals who were selected as hemochromatosis cases. Controls (n = 6896) were individuals without abnormal saturation values and without any instances of hemochromatosis ICD-9 codes.

Of the full PMRP cohort, ∼10 000 DNAs were interrogated by high-density genotyping on the Illumina HumanCoreExome beadchip. This exome array of >500 000 markers has ancestry informative markers, a panel of identity-by-descent single nucleotide polymorphisms (SNPs), coverage of markers found to be genome-wide significant in GWASs, and excellent coverage of exonic variants. The version of the beadchip was designed and used in the AMD Consortium.³⁴  Rare variants (<1% frequency) represented 47.8% of the markers, moderately common variants (1% to 10% frequency) were 8.1% of the variants, and 44.1% of the variants interrogated were common alleles (>10% frequency). The genotyping quality control measures were previously described (call rates for each variant or individual >0.985).³⁴  Variants exhibiting departure from Hardy-Weinberg equilibrium (P < 1 × 10⁻⁶) were excluded from subsequent analyses. Additional recent studies have used data generated from this genotyping to discover susceptibility genes for common diseases.³⁵  Following quality control procedures, 413 701 variants remained for analysis. The site frequency spectrum of the resulting variants is displayed in supplemental Figure 2 and supplemental Table 1.

In general, gametic phasing is necessary to directly determine compound heterozygous individuals at a particular gene. Using all 10 000 exome-genotyped samples from the PMRP, the software package Beagle was applied to infer phased haplotypes from unphased genotype data using a localized haplotype-cluster model algorithm.³⁶  The calculations were performed on a high-performance computing cluster housed at the Marshfield Clinic. As the subsequent analyses were gene-based and the genotyping data were concentrated on exonic variants, each gene in the exome was phased separately using this approach. Although rare variants can present difficulties in phasing, the use of a large sample size from a highly homogeneous population aids in mitigating the error rate. Notably, Beagle has been shown to have error rates in phasing between 0.77% to 0.94% for medium (n = 1000) to large (n = 5000) sample sizes using a 500K GWAS array.³⁶  Recently, the switch error rate was calculated for the Beagle applied to 2 sequencing data sets. Beagle attained a switch error rate of 1.525% and 0.488% for the 1000 Genomes Project and Haplotype Reference Consortium, respectively.³⁷ 

Following the phasing of the genotype data, putative functional variants were identified. The putative functional variants included in the analyses satisfied the quality control criteria as described previously.³⁴  Markers used in the analyses were either GWAS-significant as of June 2015 and/or annotated as missense, nonsense, 3′ untranslated region, 5′ untranslated region, or occurring within a splice site region. Additionally, on the resulting set of variants, annotation was performed using the ANNOVAR software.³⁸  Only variants that were also annotated as pathogenic by at least 2 ANNOVAR predictions were included in the scan. Following filtering for putative functional variants, 129 556 SNPs remained for use in our gene-based recessive diplotypes scan. The putative functional variants of FGF6 are shown in supplemental Table 2.

At each gene, individuals were classified as having a recessive diplotype configuration if they carried at least 1 putative functional allele on each homolog (PF). Individuals carrying at least 1 homolog free from putative functional alleles were deemed as having a wild type (WT) diplotype (W). The total number of case/control individuals carrying a recessive diplotype was denoted by PF_cs and PF_ct, respectively. Similarly, the total number of case/control individuals carrying a WT diplotype was denoted by W_cs and W_ct, respectively. Following the determination of these counts, a Fisher exact test was applied. As the hypergeometric null density holds for all sample sizes, the Fisher exact test is robust to imbalance between case and control sample sizes. Simulations have recapitulated this finding showing that the Fisher exact test does not inflate type I error rates under unbalanced designs.³⁹  Individuals carrying 1 or more homozygous genotype(s) at a single site for a putative functional allele were included in the PF category. Genes without any high-quality, putative functional alleles across all samples were removed from the analyses. Across all genes with analyzable data, a conservative experiment-wise multiplicity correction was calculated using 15 900 gene-based tests. To compare the recessive diplotype analysis procedure to a standard rare variant gene-based test, the RVTESTS software,⁴⁰  which implements the sequence kernel association test, was also applied to the genotype data.⁴¹  Additionally, to investigate the sex-specific effects, the Haldane odds ratio (OR) was calculated separately for female and male strata. Lastly, the Mantel-Haenszel test of homogeneity was calculated to determine the level of statistical evidence for sex-specific differences in effects.⁴² 

To explore the efficacy of the proposed approach, we performed analytic power calculations under the alternative model of compound heterozygosity/recessive inheritance of disease at 2 sites, each segregating 2 alleles. By doing so, we sought to compare the power of a standard GWAS analysis (Armitage trend test) to a log-likelihood ratio G test.⁴²  Supplemental Figure 1 shows the power of each of these tests across different sets of penetrances and haplotype frequencies. To consolidate the different sets of haplotype frequencies, the results are plotted as a function of linkage disequilibrium between the 2 sites. The power of the G test for recessive diplotypes exceeded the power for the Armitage trend test across virtually all of the parameter space. Additional work in this area was recently performed by Sanjak et al showing similar results.⁴³ 

Amino acid sequencing of the core iron-metabolism genes were collected, including transferrin receptor 1 (TFRC), FTH1, IREB1, SKP1, SKP1, ACO1, TFR2, TF, HMOX1, ACO2, HAMP, FGF6, and FGFR1. The alignments were derived from National Center for Biotechnology Information (NCBI) BLASTn database. Phylogeny for different genes were compared, showing the earliest evolutionary time point; then, occurrence for each gene was mapped to the phylogenetic tree.⁴⁴  Protein-protein interaction network inference was conducted to fibroblast growth factor 6 protein (FGF-6) and main iron-metabolism proteins. The finial network was tuned after removing nonnecessary nodes between FGF-6 and key iron molecules including hereditary hemochromatosis protein (HFE) and SLC40A1.

Colon cancer cell lines (HCT8 and HCT116), a kidney cancer cell line (786-O), a liver cancer cell line (HepG2), and a fibroblast cell line (HFF-1) were cultured in Dulbecco’s modified Eagle medium (DMEM) supplemented with 10% fetal bovine serum (FBS) at 37°C in a 5% CO₂ humidified incubator. To investigate the change in iron uptake under different protein treatments or plasmid transfection, the ferric ammonium citrate (FAC) cell culture method was used. Cells were cultured in normal DMEM and FBS medium with the presence of 10 μM FAC and 500 μM ascorbate for 48 hours during detection of cellular iron concentration. Cells cultured in normal medium exhibited very low iron concentrations. Total intracellular iron concentration in cells cultured with FAC for 48 hours was dramatically increased over cells cultured in normal medium. Recombinant human FGF-6 protein (active) and anti-ferritin were purchased from Abcam, Flag tag antibody was purchased from Proteintech Group, and anti–glyceraldehyde-3-phosphate dehydrogenase (GAPDH) antibody was purchased from Shanghai Yeasen Biotechnology. FAC (1 mM) and ascorbate (50 mM) were dissolved in distilled water. NaOH, HCl, KMnO₄, ferrozine, neocuproine, ammonium acetate, ascorbic acid, and FeCl₃ were purchased from Beijing Oka Biological Technology. Plasmid with raw FGF6 sequence was purchased from PPL-Shanghai Co, Ltd, which was constructed in an N-terminal p3XFLAG-CMV vector, whereas 3 different FGF6 mutations (E127X, D174V, R188Q) were synthesized with overlapping polymerase chain reaction (PCR).

Total intracellular iron content was measured by the ferrozine assay.⁴⁵  Cells were cultured in 12- or 24-well plate for 48 hours and washed 3 times with cold phosphate-buffered saline (PBS). After being lysed for 2 hours with 50 mM NaOH, 100 μL of cell lysates was mixed with 10 mM HCl, and 100 μL of the iron-releasing reagent (a freshly mixed solution of equal volumes of 1.4 M HCl and 4.5% [wt/vol] KMnO₄ in H₂O). The mixtures were incubated for 2 hours and 30 μL of iron-detection reagent (6.5 mM ferrozine, 6.5 mM neocuproine, 2.5 M ammonium acetate, and 1 M ascorbic acid) was added; after 30-minute incubation, 280 μL of solution was added to a 96-well plate and read 550 nm on a microplate reader. In addition, FeCl₃ (0-100 μM) was used as iron standards and protein quantification was determined by a Lowry protein assay.

Cell lysates were harvested when incubated with iron for 48 hours, then equal amounts of protein from every sample were subjected to 12% sodium dodecyl sulfate–polyacrylamide gel electrophoresis and then transferred to polyvinylidene difluoride membranes. After blocking with 5% bovine serum albumin (BSA), the membranes were incubated with GAPDH (1:10000), ferritin (1:1000), and Flag (1:2000) at 4°C overnight. Then, membranes were washed 3 times with Tris-buffered saline plus polysorbate 20, incubated with anti-rabbit or anti-mouse secondary antibody. The bands were visualized using Image QuantTL software.

Cells were washed with PBS 3 times, fixed with 4% glutaraldehyde for 10 minutes, and incubated at 37°C for 60 minutes with 2 mL of Prussian blue solution comprising equal volumes of 2% hydrochloric acid aqueous solution and 2% potassium ferrocyanide (II) trihydrate. After the cells were stained with 0.5% neutral red for 3 minutes, iron staining was visualized by Nikon microscope. Iron-positive high positive staining cells divided by total cell number was used to evaluate the iron-deposition levels.

Total RNA was extracted from the cells using TRIzol (Invitrogen). One microgram of total RNA was subjected to complementary DNA synthesis using the High Capacity cDNA Reverse Transcription kit (Applied Biosystems) according to the manufacturer’s instructions. The specific primers for each gene were designed using Primer 5 and synthesized by Generay Biotech Co, Ltd. The reverse transcription polymerase chain reaction (RT-PCR) amplification was conducted using a SYBR Green I PCR kit (TaKaRa) according to the manufacturer^’s instructions. The reaction was carried out on an ABI Prism 7900 Detector System (Applied Biosystems). RT-PCR conditions were 95°C for 3 minutes, followed by 40 cycles of 95°C for 15 seconds, 60°C for 40 seconds; the conditions for obtaining the dissociation curve were 95°C for 15 seconds, 60°C for 15 seconds, 95°C for 15 seconds. The data obtained from the assays were analyzed with SDS 2.3 software (Applied Biosystems). For each sample, the relative gene expression was calculated using a relative ratio to GAPDH. RT-PCR primers can be found in the supplemental Table 3.

The primary antibody used was anti-FGF-6 (1:200; BBI) and anti-ferritin (1:100; abcam). Liver and skin tissues from 4 liver cancer patients and 6 systemic sclerosis (SSc) patients, respectively, and normal controls were formalin-fixed and paraffin-embedded. Sections were deparaffinized and incubated with 5% BSA for 60 minutes. Cells positive for FGF-6 were detected by incubation with the primary antibody for 2 hours at room temperature followed by incubation with 3% hydrogen peroxide for 10 minutes. Rabbit anti-rabbit immunoglobulin G labeled with horseradish peroxidase was used as secondary antibody. The expression of FGF-6 or ferritin was visualized with 3,3-diaminobenzidinetetrahydrochloride (DAB-4HCl). The expression of FGF-6 in SSc and tumor tissues was quantitated by the average optical density (AOD) of the positive signal in each sample using NIH ImageJ software (Windows and Java-1.8.0).

To discover novel iron-overload–predisposing genes, we conducted a gene-based scan for recessive diplotypes composed of putative functional alleles across the exome using biobanked samples linked to electronic medical records obtained from a rural, genetically homogeneous population in central Wisconsin. Of the 10 000 samples evaluated, our transferrin saturation and diagnostic code-based phenotype algorithm identified 18 case individuals and 6896 controls. We estimated gametic phase on all individuals and restricted our analyses of diplotypes to putative functional variants. Our recessive diplotype scan identified 2 exome-wide significant genes (Figure 1; Table 1; supplemental Figure 3), HFE (P = 1.29 × 10⁻⁸; OR = 28.7) and FGF6 (P = 1.99 × 10⁻⁶; OR = 22.8). For comparison, the SKAT/rvtest procedure on the FGF6 genotype data yielded an asymptotic P = 3.86 × 10⁻⁵ and permuted P = 1.0 × 10⁻⁴. Notably, the recessive diplotype scan result exceeded exome-wide significance, whereas the rare variant test did not. Additionally, there was no statistical evidence of effect differences between females and males for the FGF6 data (Mantel-Haenszel test of homogeneity P = .728). These results motivated our investigation of FGF-6 function and the impact of specific FGF6 variants on iron metabolism.

To explore the involvement of FGF6 in iron metabolism, we found evidence for FGF-6 interactions with FGFR1, MAPK1/3, INS, FN1, and involvement in the iron-metabolism subnetwork involving transferrin (TF), HFE, hepcidin antimicrobial peptide (HAMP), and SLC40A1 (supplemental Figure 4) by investigating FGF-6 protein-protein interactions. FGF-6, also known as heparin secretory-transforming protein 2 or heparin-binding growth factor 6, has multiple heparin-binding sites (HBSs). Three known nonsynonymous variants located in the HBSs (R188Q) or flanking sites (D174V and E172X) were speculated to be important for FGF-6 function. Furthermore, D174V and E172X are located in the regions between FGFR-binding region (FGFR-BR-3) and HBS-1 (Figure 2). Hence, we studied these 3 variants in functional studies to further investigate the involvement of FGF-6 in iron metabolism.

To investigate the potential mechanism linking FGF-6 to iron metabolism, the effects of FGF6 on iron uptake and the expression of iron-metabolism genes in HepG2, HCT8, HCT116, 786-O, and HFF1 cells were evaluated. Using cultured cells and a ferrozine assay to detect iron, total intracellular iron concentration was significantly decreased in HepG2, 786-O, HCT8, HCT116, and HFF-1 cells when treating with active FGF-6 protein in a dose-dependent manner (Figure 3; supplemental Figure 5). Testing the effect of FGF6 on the expression of a set of genes involved in iron metabolism (HAMP, HDAC2, HMOX1, TFRC, and HEPH), HepG2 cells were subjected to treatment from control or FGF-6 protein and FGF6 messenger RNA (mRNA) or control and mRNA expression relative to GAPDH was measured in the 5 iron-metabolism genes. RT-PCR analysis revealed that HAMP and HDAC2 mRNA levels were significantly increased after the FGF-6 active protein introduction in HepG2 cells compared with treatment with PBS as control (Figure 4A). FGF6 plasmid transfection significantly increased HAMP, HDAC2, and HMOX1 levels, whereas TFRC levels significantly decreased in HepG2 compared with a vector without FGF6 (Figure 4B). HEPH expression did not change with either FGF6 plasmid or FGF-6 protein, suggesting that the effect of FGF-6 may be independent of HEPH (Figure 4A-B).

To investigate the effects of the FGF6 alleles on FGF-6 function, we transfected plasmids carrying either the WT FGF6 or variant FGF6 with each of the 3 point mutations described in Figure 2 (supplemental Figure 6). The M1 (E172X) and M3 (R188Q) variants exhibited a significant downregulation of HAMP mRNA compared with WT in HepG2 cells (Figure 4C), HCT-116 cells (Figure 4D), 786-O (supplemental Figure 7), A498 (supplemental Figure 7), but not HCT-8 cells (supplemental Figure 7). Evaluating the effect of M2 (D174V) on HAMP expression compared with WT only yielded a significant reduction in HepG2 (Figure 4C), but not in any of the other cell lines (Figure 4D; supplemental Figure 7). Furthermore, we noted HAMP mRNA levels in M1 and M3 transfections were comparable to control levels, which illustrated a strong attenuation of FGF-6 function for M1 and M3 variants (Figure 4C-D; supplemental Figure 7). Examining the impact of specific variants on intracellular iron concentration in HepG2 and HCT-116 cells, M1 and M3 produced significantly elevated iron deposition (Figure 4E-F; supplemental Figure 8) and ferritin expression (Figure 4G-H; supplemental Figure 8), indicating a deficiency in M1/M3 FGF6-mediated iron homeostasis. In addition, the intracellular iron-accumulation pattern was confirmed by immunohistochemistry (IHC) using Perls stain (supplemental Figure 9). In contrast, M2 did not produce a significant departure from WT in iron concentration and ferritin expression in HepG2 and HCT-116 (Figure 4). Furthermore, TFRC expression was significantly upregulated in the presence of M3 compared with WT (Figure 4C). The functions mentioned were also validated in HFF-1 (supplemental Figure 8).

We hypothesized that FGF-6 might be involved in human autoimmune diseases and cancers because abnormal iron metabolism has been reported in numerous studies.^46-50  More specifically, decreased hepcidin has been implicated in the anemia of chronic disease, which frequently accompanies these systemic inflammatory states. To explore the relationship between FGF6 expression and iron deposition in autoimmune tissues, FGF6 expression and iron deposition in the skin lesions from SSc patients and healthy controls were examined. We found significantly decreased FGF-6 protein by immunohistochemistry assay (Figure 5A) and elevated iron deposition in SSc skin tissue by ferrozine assay (Figure 5B), especially in the epidermis. Increased iron deposition was confirmed by Perls stain in SSc skin tissues (supplemental Figure 10A). We also investigated the relationship between FGF6 protein expression with iron deposition in liver cancer tissues. We found that FGF-6 was significantly decreased in nonmetastatic cancer lesion tissues (Figure 5C) and the increased iron deposition (Figure 5D; supplemental Figure 10B). However, increased FGF-6 expression was observed in metastatic liver carcinoma tissue (supplemental Figure 11), suggesting that FGF-6 plays different roles in oncogenesis and metastasis, analogous to transforming growth factor β.^51,52 

Iron homeostasis results from a combination of pathways and 4 main cell types: enterocyte, hepatocyte, macrophage, and erythroblast. The epidermal growth factor/epidermal growth factor receptor signaling pathway, heme production, STAT signaling, cyclic adenosine monophosphate signaling, ferritin storage, and bone morphogenetic protein–SMAD signaling are all involved in iron regulation. We conducted an exome-wide, gene-based recessive diplotype scan using putative functional variants to reveal additional genes underlying hemochromatosis susceptibility, an approach that can be widely applied to investigate complex disease susceptibility generated by compound heterozygosity and recessive single-site effects using existing exome-wide association genotype and sequencing data. Although the case sample size was very small, this novel scan identified FGF6 as being significantly associated with hemochromatosis following correction for multiple testing. FGF6 belongs to the paracrine FGF gene family and is largely expressed in skeletal muscle, which plays an important role in iron metabolism as it contains 10% to 15% of iron stores. We conducted the evolutionary analysis of FGF6 and known iron-metabolism genes including FGFR1, TFRC, FTH1, IREB1, TF, HMOX1, ACO2, and HAMP (encoding hepcidin). The appearance of iron-metabolism genes can be separated into 2 stages. TF and HMOX1, which are found in animals from Caenorhabditis elegans to Homo sapiens, indicate an origin in early Bilateria evolution (∼635 million years ago [Mya]). FGF6, FGFR1, ACO2, and HAMP can be found from Danio rerio to H sapiens, but are not present in C elegans and Drosophila, indicating emergence in early Vertebrata (∼485 Mya). The coappearance of these genes suggests possible coregulatory functions (supplemental Figure 4A). Functional experiments demonstrated that FGF-6 strongly impacted hepcidin expression, thereby playing a role in regulation of iron homeostasis. These results suggest FGF-6 mediates its effect on iron metabolism via hepcidin. The induction of hepcidin expression by FGF-6 leads to degradation of ferroportin through binding and internalization of ferroportin by hepcidin as shown in Figure 6. We additionally found that 3 FGF6 nonsynonymous variants increased intracellular iron concentrations and reduced hepcidin levels compared with WT FGF6, indicating loss of function. Interestingly, a genome-wide RNA interference–profiling study reported that knockdown of FGF6 increased transferrin-mediated endocytosis.⁵³  Rs12368351, ∼8 kb downstream of FGF6, has been associated with phosphorus levels⁵⁴ ; and 2 SNPs, rs140668749 and rs10849061, within 20 kb downstream of FGF6, are associated with migraine.^55,56  Previous studies have indicated that iron plays a role in autoimmunity and a study examining pulmonary arterial hypertension in SSc noted iron deposition in lung elastin fibers and giant cells,⁵⁷  however, epidermal iron deposition in SSc has not been previously investigated. We observed that FGF-6 is involved with iron deposition in SSc and liver cancer. Together, these results demonstrate that FGF receptor (FGFR) signaling through FGF-6 is a critically important mechanism in iron metabolism.

On 15 April 2014, the study was reviewed and approved by the Marshfield Clinic Research Institute Institutional Review Board (FWA00000873, IRB00000673; title: Two Allele Loss of Function Genotype Array Study). Marshfield Clinic received a Certificate of Confidentiality from the National Institutes of Health.

This work was supported by the National Institutes of Health Clinical and Translational Science Award program, previously through National Center for Research Resources grant 1UL1RR025011 and National Center for Advancing Translational Sciences grant 9U54TR000021, and now by National Center for Advancing Translational Sciences grant UL1TR000427. This work was also funded by the National Natural Science Foundation of China (31521003, 81770066) and the 111 Project (B13016). Additional funding was provided by Marshfield Clinic Research Institute grant SCH10218 and generous donors to the Marshfield Clinic.

The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Contribution: S.G. performed analyses, interpreted results, designed the functional experiments, and aided in drafting the manuscript; S.J. conducted molecular and cell biology experiments; N.E. interpreted results, provided hematology and pathway expertise, and reviewed and edited the manuscript; M.M. aided in the analyses, and with the experimental design, and reviewed the manuscript; M.W., Y.M., and W.W. provided clinical and biochemistry advice and aided in drafting the manuscript; Z.Y. performed initial genetic and statistical analyses, managed data, and reviewed the manuscript; B.O. implemented and refined the phenotyping algorithms; T.K. and J.J. aided in the regulatory paperwork and reviewed the manuscript; P.A. and F.W. generated data and aided with experiments; R.S. performed data management tasks; J.J.M. provided clinical advice and reviewed the manuscript; J.K.M. supervised the management of biological samples for genotyping and reviewed the manuscript; L.J. reviewed the manuscript and provided general scientific advice; J.A.S. provided molecular and cellular biology advice as well as clinical advice and reviewed and edited the manuscript; J.W. supervised the functional experiments, reviewed the manuscript, and provided biological advice; and S.J.S. designed the study, supervised the genetic analyses, developed phenotyping algorithms, developed analysis methods and power calculations, interpreted results, and aided in drafting and editing the manuscript.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: Jiucun Wang, School of Life Sciences, Fudan University, Shanghai, China; e-mail: jcwang@fudan.edu.cn; or Steven J. Schrodi, Center for Human Genetics, Marshfield Clinic Research Institute, 1000 N Oak Ave–MLR Marshfield, WI 54449; e-mail: schrodi.steven@mcrf.mfldclin.edu or schrodi@wisc.edu.