The study of molecular genetics has identified more than 6000 different monogenetic diseases, mostly through familial linkage analysis. Autosomal dominant, autosomal recessive, and X-linked genetic diseases with complete penetrance were the first and easiest to identify by linkage analysis. Later, genome-wide association studies (GWAS) were developed to identify genes contributing to incompletely penetrant and poly-genetic diseases. Both methods have successfully identified key genetic contributors to various diseases. However, these genetic analyses have also had many failures. Studies of thousands of individuals have failed to uncover significant genetic heritability for numerous diseases (the “missing heritability” problem), leading to a re-evaluation of the basic tenants of GWAS (reviewed in Marjoram et al., 2014).

Family linkage analysis is very good at identifying rare mutations with large effects that confer susceptibility to common diseases. Examples include the BRCA1 and BRCA2 genes in breast cancer (Easton et al., 1995; Ford et al., 1998). By comparison, GWAS studies excel at analyzing complex diseases and disease genes with a weak effect (Risch and Merikangas, 1996). The most common target in GWAS studies has been the single nucleotide polymorphism (SNP) (Palmer and Cardon, 2005). SNPs identified in GWAS studies have largely been shown to be representative of loci of interest and thus have an indirect association with the actual causal variant(s). For example, the FTO gene has been shown through GWAS studies to associate with obesity. However, its effect on obesity was eventually shown to be through an enhancer mechanism on downstream genes, including IRX3 and IRX5 (Smemo et al., 2014; Claussnitzer et al., 2015). A major advantage of GWAS studies is the unbiased inclusion of the entire genome, thus allowing for the identification of novel loci of interest and the development of original hypotheses to be tested. Once statistically significant associations have been established with one or more loci, they can be confirmed with data from other GWAS populations or through determination of the causal relationship by experimental means (Sekar et al., 2016).

While GWAS studies have produced important discoveries, their experimental designs impose hurdles that can present significant challenges. The first hurdle is population size. GWAS studies are, by definition, very large, requiring scans of large numbers of SNPs. In order to reach statistical significance, sample sizes in the thousands are required. For common diseases, very large sample sizes are attainable; however, for rarer diseases, sample size becomes problematic and is prohibitive in some cases. The second hurdle is that GWAS studies focus on identifying one or sometimes a small group of mutations in loci that are associated with a disease phenotype. It is then left to interpretation what, if any, relationship there is between these loci and the condition being studied. The third hurdle is that GWAS studies largely operate on the common disease–common variant hypothesis. This hypothesis states that the genetic contributions to the susceptibility of common diseases are attributable to a small number of variants present in more than 1–5% of the population (Manolio, 2010). While this hypothesis has proven useful for common diseases, it is not as applicable to rare diseases and more complex diseases whose mutations may not reach the 1% threshold. The fourth and final hurdle is that GWAS is not well suited to diseases in which the genetic impact on a disease state is largely dependent upon epistatic interactions. Epistatic interactions are described as the interaction between two loci whereby the phenotype of one locus depends on the genotype of the second locus (Carlborg and Haley, 2004). For quantitative traits, epistasis describes a situation where the sum of the single locus effects cannot be used to predict the phenotype of a particular genotype (Carlborg and Haley, 2004). When considering the effect of genetic interaction at the pathway level, one or two small effect mutations may not cause a noticeable phenotype, but 10 or 20 mutations could combine to create a significant “load” on the pathway and cause deterioration in functionality, leading to a diseased state. Epistatic interactions could also occur within a single gene where one SNP functionally affects another. In this model, neither mutation alone causes an effect, but together they combine to act synergistically to cause an effect greater than would be predicted by an additive model.

The GWAS method has identified a number of genetic associations that are important to numerous complex diseases. Within these predictions, often a large portion of the heritability remains to be uncovered. Understanding this “missing heritability” has become one of the major goals of clinical genetic research. Some of the tactics that have been proposed include combining linkage and association analysis (Ott et al., 2011), correcting an overestimation in heritability (Zuk et al., 2012), and considering gene–gene interactions (Liu et al., 2012).

A particularly intriguing approach to resolving the missing heritability problem was recently presented by Huebinger et al., who used a PGL risk score (Huebinger et al., 2010) to predict the outcome of patients exposed to severe burn trauma. The PGL risk score was defined as the “sum total of mutant alleles present at risk loci within a common biological pathway.” The PGL risk score is calculated using a number of assumptions that require certain concessions to be made in order to be both useful and accurate. Others have used similar pathway analysis methods to study genetic risks for coronary infarction (Yiannakouris et al., 2006; Trichopoulou et al., 2008) and Parkinson's disease (Lesnick et al., 2007). The PGL method has the potential to offer an alternative method to GWAS for identifying relevant genetic contributors, especially for lower incidence diseases.

Kernicterus or CBE is a rare disease in affluent countries because of preventive measures. However, neonatal jaundice is very common, occurring in 60–80% of full term (≥37 weeks' gestational age) otherwise healthy newborns. About 7% of jaundiced infants receive phototherapy (Bhutani et al., 2004a,b, 2013; Keren et al., 2009). Classical kernicterus is a well-described clinical tetrad of (i) abnormal motor control, movements, and muscle tone, (ii) an auditory processing disturbance with or without hearing loss, (iii) oculomotor impairments, especially impairment of upward vertical gaze, and (iv) dysplasia of the enamel of deciduous (baby) teeth, an abnormality leading to pitting, flaking, and chipping of the enamel (Shapiro et al., 2006). In addition to classical kernicterus, subclasses of kernicterus are characterized as auditory predominant, motor predominant, and subtle kernicterus (BIND) (Shapiro, 2010).

While efforts to prevent dangerous hyperbilirubinemia have been successful in the developed world, a high incidence of glucose-6-phosphate dehydrogenase deficiency, delays in seeking care, lack of TSB screening resources, and ineffective treatments in low and middle income countries ensure that bilirubin encephalopathy remains a serious world health problem (Olusanya et al., 2014; Slusher et al., 2015). Ogunlesi et al. (2007) reported a 2.7% incidence of bilirubin encephalopathy among neonates admitted to two hospitals in Nigeria between 2002 and 2005. This high prevalence of kernicterus, in conjunction with new reports of the more subtle predominantly auditory and predominantly motor forms of kernicterus (Shapiro and Popelka, 2011), underscores the need to better understand why hyperbilirubinemia causes permanent brain damage in some infants and not in others. The relative rarity of kernicterus in the West and, to a lesser extent, milder forms of bilirubin encephalopathy, makes it exceedingly difficult to perform traditional genetic analysis through large population GWAS studies. In addition, the apparent wide range of susceptibility to a common bilirubin exposure points toward polygenetic determination of risk for neurotoxicity. The current understanding of the genetics of bilirubin encephalopathy (described below) also points toward a polygenetic cause for increased risks to bilirubin neurotoxicity.

Together, these factors indicate a need for an alternative method, such as the PGL risk score, to understand genetic factors contributing to risk for CBE. The PGL method is a good fit for bilirubin encephalopathy because bilirubin encephalopathy results from a one-time insult, similar to the severe burns examined in Huebinger et al.'s study. For bilirubin encephalopathy that insult is in the form of moderate to extreme hyperbilirubinemia at birth. Also, similar to the total body surface area measure for a serious burn, the body's response to bilirubin toxicity is variable and not well explained by its key diagnostic test, peak TSB levels.

While significant progress has been made in understanding the various causes of neonatal hyperbilirubinemia, few studies have been undertaken to understand the mechanism of bilirubin toxicity. Multiple cellular and molecular cascades likely underlie bilirubin-induced neuronal injury, including plasma membrane perturbations, excitotoxicity, neuroinflammation, oxidative stress, and cell cycle arrest. It is likely that variable effectiveness of the cellular response to bilirubin exposure contributes to the range of observed clinical outcomes. To explain part of this variability, we hypothesize that a combination of SNPs in key pathways are responsible for creating a vulnerable state to bilirubin neurotoxicity. The development of this vulnerable state likely occurs through a combination of increased CNS bilirubin levels and impaired cellular response to bilirubin.

On the basis of the current understanding of bilirubin encephalopathy, we suspect that the genetic impact is the result of the combined effect of mutations with relatively small effects that, if present alone, would otherwise be clinically benign. Alternative hypotheses that could also explain the TSB's poor predictive value for kernicterus include disrupted albumin/bilirubin binding capacity, reductions in available albumin, and use of TSB levels as the sole diagnostic test for bilirubin exposure. A number of studies have been published that examine the issue of TSB's reliability in both predicting and improving outcomes. These studies have compared the effectiveness of alternatives to TSB such as albumin/bilirubin binding capacity (Kapitulnik et al., 1974), unbound (free) bilirubin levels (Cashore and Oh, 1982; Ahlfors et al., 2009; Morioka et al., 2015), and bilirubin-albumin ratios either in place of or in addition to TSB testing (Ahlfors, 1994; Amin et al., 2001, 2016).

While Nakamura et al. showed over 30 years ago that free bilirubin levels predict neurotoxicity better than TSB levels do (Nakamura et al., 1985), technical issues, reliability, and lack of availability of these tests have called their usefulness into question (McDonagh and Maisels, 2006; Hulzebos and Dijk, 2014). The use of the bilirubin/albumin ratio has also been shown to be as strong a predictor of neurotoxicity as TSB (Iskander et al., 2014) but has not been shown to offer an advantage in predicting or improving outcomes over TSB alone (Iskander et al., 2014) or when used in addition to TSB (Hulzebos et al., 2014).

Bilirubin is largely a breakdown product of hemoglobin. Bilirubin rises after birth in all mammals, and in humans usually peaks at 3–5 days after birth and then falls to very low values through the rest of life (Kivlahan and James, 1984; Maisels et al., 2014). It is bound in the blood to proteins, especially albumin. Very little bilirubin normally reaches brain tissue. There is evidence that at low levels bilirubin is beneficial, acting as a natural antioxidant (Adin et al., 2005; Bakrania et al., 2016), but very high levels can exceed the buffering capacity of the blood and allow bilirubin to move into the brain, causing a unique pattern of brain injury (Shapiro et al., 2006; Shapiro, 2010).

Hyperbilirubinemia results from the excessive breakdown of red blood cells (RBC) (hemolysis) or impaired elimination of bilirubin. Most genetic abnormalities linked to the presentation of bilirubin encephalopathy are related to production and elimination. Hemolysis may occur as a result of immune mechanisms (e.g., Rh disease, ABO blood group incompatibilities), increased fragility of RBC [Glucose-6-phosphate dehydrogenase (G6PD) deficiency], or from bleeding (e.g., cephalohematoma). Decreased elimination is enhanced in infants with genetic variations in the enzyme or the promotor regions of UDP-glucuronosyltransferase 1A1 (UGT1A1), (Watchko and Tiribelli, 2013). Other factors that contribute to hyperbilirubinemia include decreased oral intake and stooling. It is believed that hyperbilirubinemia due to decreased oral intake may account for a substantial number of readmissions of newborns to hospitals after discharge (Seidman et al., 1995; Farhat and Rajab, 2011; Maisels et al., 2014).

Possible sources of outcome variability include factors such as the source of hyperbilirubinemia, the gestational age of the infant, the duration of exposure, the kinetics of bilirubin's movement between miscible compartments of the body (e.g., blood, tissue, CNS, extracellular, and intracellular), the type of treatment administered and other known risk factors such as acidosis or sepsis (Shapiro, 2010; Hulzebos and Dijk, 2014). Furthermore, a significant genetic contribution to the susceptibility of an individual to neurological damage due to hyperbilirubinemia may be related to the disruption of mitochondria, which may in turn affect cellular and molecular cascades including plasma membrane perturbations, excitotoxicity, neuroinflammation, oxidative stress, and cell cycle arrest referred to previously (Rodrigues et al., 2002a,b; Vaz et al., 2010; Brites, 2011; Barateiro et al., 2012).

Clinical observations have shown that in cases where TSB reaches extremely high levels (>40 mg/dL or >680 μM¹), the likelihood of developing at least some form of kernicterus approaches 100%, although exceptions to this generalization have been published (Newman et al., 2003, 2006). These clinical observations likely indicate that for most neonates, there is a maximum level of bilirubin or “tipping point” where the body's natural defense mechanisms for clearing unbound bilirubin from susceptible neurons are overwhelmed and deposition occurs. However, because kernicterus has been known to develop at lower TSB levels, it is clear that the maximum level of bilirubin exposure varies from person to person (Ritter et al., 1982; Odutolu and Emmerson, 2013). Even allowing for external factors and the inexact nature of TSB as a measure of exposure to bilirubin in brain tissue, there is still a strong likelihood that genetic variability contributes to the overall susceptibility of a child to bilirubin neurotoxicity and subsequent chronic adverse neurodevelopmental sequelae.

Here we propose using the PGL method in order to determine the extent to which a multifactorial genetic response contributes to the variable susceptibility to bilirubin toxicity and the resultant variability in disease outcome. Due to the lack of a CBE GWAS to guide our pathway selection, we propose creating three custom “pathways” of genes that are organized into “tiers” based on the mechanism of physiological contributions and the relative strength of evidence that these genes contribute to bilirubin susceptibility and sensitivity. Each of these three gene tiers will be analyzed individually and independently using the PGL method.

The PGL method differs from traditional association analysis in a number of important ways. First, the PGL score uses a biased approach by leveraging what is already known about a disease to focus on a specific set of genes (pathway) and SNPs. Second, unlike GWAS and linkage analysis, the PGL score method involves the calculation of a risk score. Third, the method can be performed using a smaller sample size than in GWAS studies because the number of variables is dependent on candidate pathway size, thus reducing the required sample size to achieve adequate statistical power. Finally, the PGL method allows detection of combinatorial epistatic interactions that would otherwise have been missed by both GWAS and linkage analysis.

The decision to use a biased approach comes with the understanding that the available knowledge will very likely present only a partial picture of the entire genetic heritability. However, the hope is that if enough of the biological background is known, then the analysis will prove to be useful and actionable. By accepting a higher rate of false negatives, we will find it is then possible to drastically reduce the number of genes and thus variants being surveyed. Reducing the number of SNPs being investigated allows for a corresponding reduction in sample size, opening up analysis to rarer diseases. Huebinger et al. used their knowledge of injury and infection to select the Toll-like receptor signaling pathway to assess whether an increased PGL risk score within this pathway could predict the outcome following burn trauma (Huebinger et al., 2010). Similarly, using the PGL method, we will perform a targeted genetic analysis on a small number of samples and will determine whether a common mutational “load” exists in our affected individuals that points toward the development of a state of bilirubin hypersensitivity.

After the pathway of interest is selected, the next step is to select the candidate SNPs from within the pathway. For this crucial step, Huebinger et al. chose to create their candidate SNP list from variants that had previously published associations with burn trauma. Once these SNPs are determined, the PGL risk score can be calculated by summing the numerical scores for genotypes across all SNPs that are homozygous protective (assigned a score of “0”), heterozygous (score of “1”) and homozygous detrimental (score of “2”). The initial calculation results in the raw risk score where each SNP is treated equally. The risk score is then weighted to better represent the impact of each of the SNPs. For this step Huebinger et al. chose to use the adjusted odds ratio from previously reported studies for each SNP between each individual locus and either death or sepsis. For loci without previously published data, the value was left unweighted. The data could then be analyzed further by Mann-Whitney U Test to investigate whether the PGL risk score alone can predict developing sepsis. In addition, the multivariate logistic regression test could be used to investigate the simultaneous effect of multiple variables as risk factors for various outcomes. A synopsis of the PGL workflow can be seen in Figure 1.

For their analysis Huebinger et al. chose to adjust for age, gender, race, total body surface area burned, and inhalation injury. They were able to show that each unit increase in weighted PGL risk scores was associated with increased probability of both complicated sepsis and death. By analyzing data from only 155 burn trauma patients and using only 6 SNPs at 6 different loci, they showed that a larger PGL risk score indicated increased odds of both sepsis and mortality in response to burn trauma. This work, along with other similar pathway load work (Lesnick et al., 2007), demonstrates that by using a targeted approach it is possible to show that a genetic load placed on a crucial pathway can have predictive capabilities even with a relatively small sample size.

To develop a SNP set for each tier, we will collect whole exome sequence data from patients with moderate to high TSB levels (10–25 mg/dL and >25 mg/dL respectively) both with and without evidence of classic auditory or motor dysfunction. We will identify positive “hits” through a pipeline utilizing the Partek Genomics Suite and Flow (Partek Inc., St. Louis, MO, USA). Causal variants in this study will then be identified through the use of QIAGEN's Ingenuity® Variant Analysis™ software (www.qiagen.com/ingenuity) from QIAGEN Redwood City. Ingenuity® Variant Analysis™ allows for multiple types of SNP filtering to create a list of screened causal variants that are present in a significant number of the case group vs. controls. The causal variant filter is based on the American College of Medical Genetics and Genomics (ACMG) categorization guidelines (Richards et al., 2015). Alternatively, or in combination with ACMG categorization, it is possible to select variants that cause a gain of function or a loss of function based on information in the Ingenuity® Knowledge Base. The Ingenuity® Knowledge Base is a “repository of expertly curated biological interactions and functional annotations created from millions of individually modeled relationships between proteins, genes, complexes, cells, tissues, drugs, and diseases” manually created from the literature (http://www.ingenuity.com/science/knowledge-base). The causal variant list can also be filtered to include only those genes present in each of the gene set tiers.

Once the key SNPs predicted to be deleterious are identified, the PGL risk score can be calculated. The raw score will be calculated from the individual genotypes as described above as homozygous protective = 0, heterozygous = 1 and compound heterozygote/homozygous deleterious = 2. The raw score will then be summed and the adjusted odds ratio, confidence interval, and p-value will be calculated. Multivariate logistic regression will then be used to evaluate the simultaneous effects of multiple variables as risk factors for the development of one of three outcomes. Possible outcomes include predominantly auditory dysfunction, predominantly motor dysfunction, and auditory plus motor dysfunction (classic kernicterus) (Shapiro, 2010). Outcomes will be further analyzed with adjustments for TSB, albumin, gender, race, sepsis, and gestational age at birth.

For weighted analysis, the SNPs will be weighted based on the Sorting Intolerant from Tolerant (SIFT) score (Ng, 2003) or the Combined Annotation-Dependent Deletion (CADD) score (Kircher et al., 2014). By utilizing the ACMG score for initial selection and the SIFT or CADD scores for weighting, we are able to evaluate the direct effect of a SNP, a crucial difference from the GWAS based method. A synopsis of the proposed PGL workflow to be used for bilirubin encephalopathy can be seen in Figure 2.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.