Multiplex alcohol dependence families were identified based on the presence of a same-sex sibling pair, each with a diagnosis of alcohol dependence by DSM-III criteria, the prevailing DSM at study initiation. One member of the proband pair was recruited from a treatment facility in the Pittsburgh area. Those providing permission to contact relatives for preliminary screening were provisionally included. Following screening, only families who were free of other psychiatric comorbidity (other than alcohol dependence) in the first and second-degree relatives of the proband pair were included. Families were excluded if recurrent major depression, bipolar disorder, schizophrenia, or a primary substance use disorder (SUD) other than alcohol dependence was present in either the proband pair of adult alcohol dependent individuals or their first degree relatives. Alcohol dependence was required to have occurred at least 1 year before the onset of the SUD in all individuals studied. Accordingly, those cases with a comorbid SUD present were considered to be secondary disorders.

All participants in the family study (the same sex proband pair of siblings, other siblings, and parents of the siblings) including those with sequencing were administered a structured, in-person, psychiatric interview [Diagnostic Interview Schedule (DIS)] (2) by a trained Masters-level clinician, allowing for a diagnosis of alcohol dependence by DSM-III (3) (the diagnostic system in place when the study was initiated) to be determined. Each participant was administered questions necessary to determine whether the participant met Feighner criteria (4) for alcohol dependence. The Feighner criteria diagnosis requires at least one symptom in three out of four symptom categories. The categories were: Category 1 – medical consequences (e.g., withdrawal symptoms or health problems); Category 2 – attempt to control drinking (e.g., limit use to certain times of the day, drinking before breakfast, and use of non-beverage alcohol such as mouth wash); Category 3 – legal/social problems (e.g., arrests, fights, and DUI); and Category 4 – the participant’s report of excessive drinking (e.g., family/friends object to participant’s drinking and the participant felt guilty about their own drinking). The earliest age the participant had at least one symptom in each category was determined. The age at which three problem categories were positive was considered to be the age of alcohol dependent use.

The interview data was supplemented with an open-ended clinical interview by a second clinician in order to provide reliability of the DIS and Feighner criteria interviews. A best-estimate diagnosis of alcohol dependence was made using the information from the DIS, the Feighner criteria, the second open-ended interview, and any family history information provided by the participating relatives. For the present analysis, DIS items were rescored to be consistent with DSM-V and the average number of symptoms accumulated to provide current validity of the diagnoses obtained using the earlier DSM-III nomenclature.

Although each family included a proband pair with an AUD, for the purpose of the present analysis, only one member of the pair was included to insure that unrelated cases could be compared with unrelated controls. The demographic characteristics for these individuals are presented here. The probands included 26 women and 27 men. Self-reported ancestry was 94% Caucasian and 6% African-American. They averaged 34.3 ± 6.3 years at interview.

Based on self-reported ancestry, the full control group of 8,893 individuals consisted of 58.1% Caucasian, 14.3% African-American, 12.1% Hispanic, with East Asian, Middle Eastern, and South Asian representing 0.7% of the total. Another 14.7% were of unknown ancestry (Table 1).

Whole exome sequencing of blood-extracted DNA was performed at the Institute for Genomic Medicine (IGM) at Columbia University for 282 individuals from multiplex AUD families along with IGM control subjects using the same bioinformatic pipeline. Our goal was to compare variants observed in one member of 53 proband pairs from separate families to a large control sample maintained at IGM. Sequencing was performed using the SureSelect Human All Exon (65 MB; Agilent Technologies, Santa Clara, CA, USA) or the NimbleGen SeqCap EZ version 2.0 or 3.0 exome enrichment kit (Roche NimbleGen, Madison, WI, USA) on HiSeq 2000 or 2500 sequencers (Illumina, San Diego, CA, USA) according to standard protocols. The IGM pipeline maintained a goal of achieving 10-fold sequencing coverage. The quality of sequencing was monitored using software designed to filter the raw sequencing data. Data were excluded based on several considerations including the presence of duplicate reads, and being among a list of known sequencing artifacts. Data were included only if they were among the consensus coding sequence public transcripts (CDDS Release 14). Raw sequence data was taken in as an Illumina lane level fastq files along with the reference genome (Human Reference Genome-NCBI Build 37) to which it is aligned using the Burrows-Wheeler Alignment (BWA) tool. The new data that was generated contained the reads and the genetic location resulting in Sequence Alignment and Map (SAM) files that were then processed as a Binary Alignment and Map (BAM) files for variant calling. Identification of the variants was performed with GATK software with analyses performed using the Analysis Tool for Annotated Variants (ATAV).

We determined both the percentage of case subjects and the percentage of control subjects who had at least 10-fold coverage at each site. At least 10-fold sequencing read coverage was achieved for ≥95% of the megabase pairs (Mbp) of the consensus coding sequence (CCDS; Release 14) for the case and control subjects. The percentage of cases and controls meeting this criterion did not differ significantly insuring that results were not influenced by differential coverage in the CCDS sequence.

Individual CCDS sites were excluded from analysis if the absolute difference in percentages of case subjects compared with control subjects who achieved at least 10-fold coverage differed exceeded a predetermined threshold (<10%). Harmonization of coverage between cases and controls is an important step. Where coverage is highly imbalanced, the less well represented group can show an inability to call a particular variant. This can lead to an enrichment bias toward the better represented group (5). All collapsing tests were then performed on the pruned 30.67 Mbp of CCDS sites. This insured that both cases and controls had similar opportunity to call variants. Absence of preferential inflation of background variation was further confirmed by conducting an exome-wide tally of rare autosomal synonymous (i.e., presumed neutral) variants per individual. This tally did not find a significant difference between the case and control groups. Autosomal read depth and the sequencing coverage were consistent between the cases and control samples with no significant difference being seen.

The search for genes potentially influencing risk for AUD was implemented using genetic collapsing tests (5). We used the conventional gene based collapsing method in which the protein coding boundaries of the gene is the unit of analysis and the criteria for comparing cases and controls is based on whether there is at least one qualifying variant in the gene. This method differs from rare variant burden tests which define a genetic region and then aggregate the information within the defined region to describe the summary dose or burden. Because genetic relatedness can distort the contribution of any particular variant to the test statistic, we first pruned our sample of 282 members of densely affected multiplex families to include only one case per family resulting in 53 cases. We chose one of the members of the proband pair that had been identified through participation in a substance use treatment program and who had led us to the remaining family members. We focused our analyses on CCDS protein-coding sites with minimal variability in coverage between the case and control populations. We analyzed qualifying variants, a term originally introduced to refer to the subset of genetic variations within the sequence data (6) that meets specific population allele frequency and predicted variant effect criteria. We defined nine different qualifying variant models in all (Table 2) for additional analysis to determine if deleterious variants might be found. The primary model of interest was one in which a search for “ultra-rare” non-synonymous variants was conducted in order to capture a category of genetic variation expected to be most enriched for variants of high effect. Ultra-rare variants were identified using an external sources (Exome Variant Server and Exome Aggregation Consortium release 0.3) (7), and sequence data from our combined case and control test populations to find variants with a minor allele frequency (MAF) of less than 0.05%. Qualifying variants were restricted to indels and single nucleotide variants annotated as having either a loss-of-function (LoF) effect, an inframe indel, or a “probably damaging” missense prediction by Polymorphism Phenotyping version 2 (PolyPhen, HumDiv¹) (8). These analyses relied on the predicted effects of the LoF and missense annotated variants whose functions have not been individually confirmed in the laboratory. We subsequently performed analyses of CCDS genes using the nine alternative qualifying variant models as defined in Table 2, including an autosomal recessive model and a synonymous variant negative control model.

For each of the nine models, we tested the list of 18,666 CCDS genes. For each gene, an indicator variable (1/0 states) was assigned to each individual on the basis of presence of at least one qualifying variant in the gene (state 1) or no qualifying variants in that gene (state 0). A two-tailed Fisher’s exact test (FET) was then performed for each gene to compare the rate of case subjects carrying a qualifying variant compared with the rate of control subjects. For our study-wide multiplicity adjusted significance threshold, a Bonferroni correction was made using the number of genes tested across the non-synonymous models. The Bonferroni threshold for a single model is 2.67 × 10^–6. Because our primary interest was in ultra-rare variants this threshold is appropriate for the primary analysis. For the complete analysis of the nine models, the appropriate threshold is 2.97 × 10^–7. We did not correct for the synonymous (negative control) model due to differing sampling rates between cases and controls for genes on the X chromosome. Collapsing analyses were performed using an in-house package, ATAV.² Additional binomial analyses, logistic regression analyses, and FETs were completed using the “stats” package in R version 3.2.2 (R Foundation for Statistical Computing, Vienna, Austria).

All probands met DSM-III criteria for alcohol dependence, the diagnostic scheme in place at the time of the participant ascertainment. Additionally, they met Feighner criteria for alcohol dependence, and DSM-V criteria for AUD (Table 3). Using individual responses to the DIS interview retained in file, all probands were diagnosed with DSM-V criteria. All probands met DSM-V criteria (Table 3). Eight of the 11 conceptual areas included in DSM-V were available for scoring using responses to questions that mapped to DSM-V concepts. The mean and standard deviation for the sample was 5.88 ± 1.01 symptoms, indicating a severe form of AUD overall.

We searched 18,666 protein coding genes for a statistically significant excess of rare deleterious variation in 282 members of multiplex AUD families. The present report is based on a comparison of 53 unrelated probands (one per family) from these families in comparison to a data set of 8,893 unrelated controls similarly sequenced at IGM. Our results are presented in three sections: (1) primary analyses using all available controls (8,893); (2) secondary analyses using ancestry filtered controls (N = 2,814) providing the closest match to the cases with AUD; and (3) analyses directed by a previous search of genes identified in a GWAS analysis conducted by Peng et al. (9) that included analysis for rare variants.

Although no gene achieved genome wide significance, case enrichment for 16 genes was significant at 10^–4 and one gene at 10^–5 indicating they are plausible candidates for follow-up studies. A total of six gene variants qualify as ultra rare with MAF of 0.005 with two of these also considered potentially damaging: SELENBP1 and GMPPB. Four gene variants are ultra rare but not previously identified as damaging: CDSN, CHRNA9, IFT43, and TLR6. Two of the genes, CDSN and TLR6 are involved in the immune response. This finding is of interest because toll-like receptors are integral components of neuroimmune adaptation with immune signaling in the brain reported to influence development of AUD (12). Extensive pre-clinical work has suggested the importance immune signaling especially for TLR4, a member of the toll-like receptor family, and results of clinical trials using agents that alter immune functioning to modify alcohol use (46). IFT43 plays a role in cilia assembly and maintenance (13) but has not been reported to be related to alcohol use. CHRNA9 is of interest because of previous reports of its association with AUD (15), nicotine dependence (14), and drug dependence (16). In a study conducted by Zuo et al. data for over 26,000 participants were analyzed that included 9 different psychiatric disorders and 15 independent cohorts. All CHRN genes, with exception of those that are muscle-type genes, were associated with nicotine dependence and alcohol dependence independently and in participants with both disorders (15). Due to the large sample size and the inclusion of multiple rare variants from the CHRN constellation p-values were exceptionally strong with p-values as low as 10^–39.

The present study also found variation between probands and controls in eight genes with MAF of 0.001 identified with a Flexible PolyPhen or Flexible REVEL algorithm: ZNF514, OXGR1, DIEXF, TMX4, MTBP, PON2, ANKRD46, and CRHBP. The possible relevance to AUD is most evident in the case of CRHBP, a gene involved in the stress response and a component of the HPA axis. Both animal studies (29) and clinical investigations (34, 35) support the role of HPA axis genes in alcohol consumption. In addition to these eight genes with MAF of 0.001, differences between cases and controls were seen in an additional three genes with MAF of 0.001 with variants that also qualify as protein truncating variants: AGTRAP, ANKRD46, and PPA1. However, no direct association between these genes and alcohol consumption have appeared previously in the clinical or animal literature to our knowledge.

Using the controls filtered to provide a better match to the AUD cases, we found nine genes with statistical significance of 10^–4. Two of these nine genes were also significant at 10^–4 in the comparison with the larger unfiltered control set. One gene of special interest, CHRNA9, which was significant at 0.0018 in the ancestry controlled comparison was significant at 10^–4 in the comparison with the larger control sample of 8,983 individuals.

Gene set analysis was run on 48 genes previously reported by Peng et al. (9) with MAF between 0.03 and 0.005. Four of these genes showed significance at 0.05 in our study. These genes showing marginal significance in our sample are in agreement with the Peng et al. report and appear to hold promise for follow-up.

The major strength of this study is the sample that was utilized for assessing rare variants. The sample was ascertained using a two proband ascertainment method that required the presence of two same sex siblings with a confirmed diagnosis of AUD to include them and their family members in our family study. All probands and their family members were administered a structured psychiatric interview and diagnoses made using a multimodal approach that included multiple sources of information including results of the structured interview and a follow-up unstructured interview allowing for a best estimate diagnosis to be made. This two proband requirement for inclusion of a family resulted in selection of a set of densely affected families. Such families are ideal for uncovering rare variants. However, this requirement was labor intensive, requiring interviews with approximately 100 potential probands to be able to select one family that met all of the requirements for inclusion. Because the goal was to search for genes conferring greater susceptibility for AUD and provide a sample with minimal comorbidity, families were excluded where a first-degree relative of the proband pair had bipolar disorder, schizophrenia or major depression by DSM criteria.

Inclusion of women with AUD along with the men with AUD from high density families is a strength. An additional strength of the study is the control sample of 8,983 individuals that were sequenced at IGM using the same methodology as the case sample drawn from families with a high density of AUD. Although the control sample contains a wide variety of individuals with varied ethnicity, a secondary analysis was performed utilizing ancestry matched cases and controls that offered additional information. This analysis determined if variants identified in the cases differed from control subjects controlling for ethnic background.

There are limitations in the conclusions that can be drawn based on the small number of cases that could be analyzed though a large set of similarly sequenced controls could be included. Among these is the fact that confirmation in the ancestry filtered analysis could only be performed in a Caucasian (one of the European clusters) set of cases and controls. Analysis of within-family variation using data from all 282 cases is planned and will provide an important confirmation of the strength of conclusions about specific genes. Nevertheless, the genes identified appear to be plausible candidates for future follow-up.

A comment regarding where this report and other rare variant analyses using NGS fits into the search for genetic variation associated with AUD is needed, especially in view of the ever larger GWAS studies currently being published. It has been noted by Povysil et al. (5) that because GWAS searches among common variants, that as sample sizes become ever larger, the newly identified variants have had smaller effects on risk. Additionally, it has been argued that if a sufficient number of individuals are genotyped that the identified variants reported to be associated with complex traits or diseases would be spread broadly and densely across the genome. This omnigenic model though not the intent of individual GWAS studies may have implications for finding mechanistic explanations for disease etiology and medication development. In a review of results for GWAS studies designed to uncover genetic variation for AUD, Hart and Kranzler (47) concluded that findings for GWAS studies of AUD have been largely inconsistent with the exception of variants encoding the alcohol metabolizing enzymes. A large scale meta-analysis of 274,424 individuals (48) combined results of GWAS studies across populations finding 10 variants within genes that appear to confer increased risk for AUD. Among these 10, 4 were variants that encode alcohol metabolizing enzymes. More recently a large scale study of SUD involving over a million participants has been reported (49). This report found some variants common to many specific disorders but others that were unique. Because of the large amount of comorbidity across various SUD including AUD, it may be necessary to focus more attention on acquisition of phenotypically well characterized samples and in the collection of family data where familial comorbidity can be assessed. These steps may reduce the phenotypic heterogeneity making it possible to identify rare variants that provide a mechanistic explanation for the etiology of AUD.