Whole-genome association analyses of sleep-disordered breathing phenotypes in the NHLBI TOPMed program

Sleep-disordered breathing is a common disorder associated with significant morbidity. The genetic architecture of sleep-disordered breathing remains poorly understood. Through the NHLBI Trans-Omics for Precision Medicine (TOPMed) program, we performed the first whole-genome sequence analysis of sleep-disordered breathing. The study sample was comprised of 7988 individuals of diverse ancestry. Common-variant and pathway analyses included an additional 13,257 individuals. We examined five complementary traits describing different aspects of sleep-disordered breathing: the apnea-hypopnea index, average oxyhemoglobin desaturation per event, average and minimum oxyhemoglobin saturation across the sleep episode, and the percentage of sleep with oxyhemoglobin saturation < 90%. We adjusted for age, sex, BMI, study, and family structure using MMSKAT and EMMAX mixed linear model approaches. Additional bioinformatics analyses were performed with MetaXcan, GIGSEA, and ReMap. We identified a multi-ethnic set-based rare-variant association (p = 3.48 × 10−8) on chromosome X with ARMCX3. Additional rare-variant associations include ARMCX3-AS1, MRPS33, and C16orf90. Novel common-variant loci were identified in the NRG1 and SLC45A2 regions, and previously associated loci in the IL18RAP and ATP2B4 regions were associated with novel phenotypes. Transcription factor binding site enrichment identified associations with genes implicated with respiratory and craniofacial traits. Additional analyses identified significantly associated pathways. We have identified the first gene-based rare-variant associations with objectively measured sleep-disordered breathing traits. Our results increase the understanding of the genetic architecture of sleep-disordered breathing and highlight associations in genes that modulate lung development, inflammation, respiratory rhythmogenesis, and HIF1A-mediated hypoxic response.


Results:
We identified a multi-ethnic set-based rare-variant association (p = 3.48 × 10 −8 ) on chromosome X with ARMCX3. Additional rare-variant associations include ARMCX3-AS1, MRPS33, and C16orf90. Novel common-variant loci were identified in the NRG1 and SLC45A2 regions, and previously associated loci in the IL18RAP and ATP2B4 regions were associated with novel phenotypes. Transcription factor binding site enrichment identified associations with genes implicated with respiratory and craniofacial traits. Additional analyses identified significantly associated pathways. Conclusions: We have identified the first gene-based rare-variant associations with objectively measured sleepdisordered breathing traits. Our results increase the understanding of the genetic architecture of sleep-disordered breathing and highlight associations in genes that modulate lung development, inflammation, respiratory rhythmogenesis, and HIF1A-mediated hypoxic response.
Keywords: Sleep-disordered breathing, Sleep apnea, Whole-genome sequencing, WGS, Genome-wide association study, GWAS Background Sleep-disordered breathing (SDB) is a prevalent disorder associated with increased sleepiness, mortality, and morbidity from a wide range of cardiometabolic and other diseases [1,2]. The most common type of SDB is obstructive sleep apnea (OSA), characterized by repeated airway collapse leading to intermittent hypoxemia and sleep disruption, that is increased in prevalence with older age and male sex [2]. An estimated 936 million adults aged 30-69 have mild to severe OSA worldwide [3]. The disease is heritable and appears to be multifactorial, reflecting variable contributions of abnormalities in ventilatory control, craniofacial anatomy, and adiposity [2,[4][5][6][7]. Sleep-related hypoxemia can also be due to central sleep apnea, a less common disorder, due to a lack of respiratory drive [8]. OSA is typically measured clinically using the apnea-hypopnea index, which counts the number of total (apnea) and partial (hypopnea) breathing cessations per hour of sleep. Due to an incomplete understanding of its molecular basis, the standard OSA treatment of continuous positive airway pressure (CPAP) only addresses the downstream manifestations of airway collapse through nightly use of pressurized air to the nasopharynx, a therapy that often is poorly tolerated. Therefore, there is a critical need to identify molecular pathways that could provide specific therapeutic targets. The need for overnight studies to phenotype SDB traits has limited the available sample size for genetic analyses, and only several common-frequency genome-wide analysis studies have been reported [9][10][11]. Increased statistical power may increase the genetic resolution of regions that may not be adequately tagged by current genotyping arrays due to population differences and/or reduced linkage disequilibrium with biologically relevant regions.
The Trans-Omics for Precision Medicine (TOPMed) program is an NIH National Heart, Lung, and Blood Institute program designed to improve the understanding of the biological processes that contribute to heart, lung, blood, and sleep disorders [12]. TOPMed has generated whole-genome sequencing (WGS) data on over 100,000 individuals from multiple cohorts at > 30× depth, including seven studies with objective assessment of SDB. A variant imputation server using TOPMed data also allows for high-quality imputation of nonsequenced genotype chip data [13]. A complementary initiative sponsored by the Centers for Common Disease Genomics (CCDG) of the NIH National Human Genome Research Institute has generated sequencing data from additional individuals in two TOPMed cohorts. These initiatives provide the ability to examine the genetics of SDB at unprecedented detail in African-Americans (AA), Asian-Americans (AsA), European-Americans/Australians (EA), and Hispanic/Latino-Americans (HA).
In this first genome-wide sequencing analysis of SDB, we examine the apnea-hypopnea index (AHI), the standard clinic metric of SDB, and four complementary measurements of overnight hypoxemia: average and minimum oxyhemoglobin saturation (SpO 2 ) during sleep and the percent of the sleep recording with SpO 2 < 90% (Per90), and the average desaturation per hypopnea event. These indices were chosen because of clinical relevance, high heritability, or prior significant GWAS findings [9,11,14]. We examined 7988 individuals with objectively measured SDB and WGS data in conjunction with data from 13,257 individuals with imputed genotype data.

Methods
Each study had a protocol approved by its respective Institutional Review Board and participants provided informed consent. A study overview is provided in Additional file 2: Figure S1. There were two classes of data: "WGS studies" had WGS performed by the TOPMed program and, in some cases, in additional participants by the CCDG program (referred to as "WGS" studies); "Imputed studies" had array-based genotyping later imputed using the TOPMed imputation server (as described below). Some studies with WGS contributed imputed study data from additional array-based genotyped individuals. Ten studies were analyzed (Tables 1  and 2).

WGS studies
The Atherosclerosis Risk in Communities Study (ARIC), the Cardiovascular Health Study (CHS), and the Framingham Heart Study Offspring Cohort (FHS) included individuals who participated in the Sleep Heart Health Study (SHHS), who underwent polysomnography (PSG) between 1995 and 1998 using the Compumedics PS-2 system [15][16][17][18]. These samples included 1028 EAs from ARIC, 151 AAs and 557 EAs from CHS, and 478 EAs from FHS.
The Multi-Ethnic Study of Atherosclerosis (MESA) is investigating the risk factors for clinical cardiovascular disease [19]. PSG was obtained between 2010 and 2013 using the Compumedics Somte system [20]. This analysis includes data from 698 EAs, 486 AAs, 456 HAs, and 229 AsAs.
The Cleveland Family Study (CFS) was designed to investigate the familial basis of SDB, with four visits occurring from 1990 to 2006 [21]. Sleep was assessed either in a clinical research center using full PSG (Compumedics E series) (visit 4) or in the latest available prior examination using an in-home sleep apnea testing device (Edentrace). Data were analyzed from 505 AAs and 485 EAs (339 AAs and 234 EAs with full PSG data).
The Hispanic Community Health Study/Study of Latinos (HCHS/SOL) is studying multiple health conditions in HAs [22,23]. Home sleep apnea testing was performed during the baseline examination (2008-2011) using the ARES Unicorder 5.2, a validated device including a forehead-based reflectance oximeter, a nasal pressure cannula and pressure transducer, an accelerometer, and a microphone [24]. Two thousand three hundred thirty-nine individuals provided data.
The Jackson Heart Study (JHS) is investigating cardiovascular disease in AAs [25]. An in-home sleep study was performed from 2012 to 2016 using a validated type 3 sleep apnea testing device (Embla Embletta Gold) [26,27]. Five hundred seventy-five individuals contributed data.

Imputed genotype studies
The Osteoporotic Fractures in Men Study (MrOS) is a multi-center cohort study initially designed to examine the risk factors for osteoporosis, fractures, and prostate cancer in older males [28,29]. An ancillary study (MrOS Sleep;2003 focused on outcomes of sleep disturbances used PSG and nearly identical procedures as in MESA (Compumedics Safiro system) [30]. Two thousand one hundred eighty-one EA individuals were included, with genotyping performed using the Illumina Human Omni 1 Quad v1-0 H array.
The Starr County Health Studies (Starr) investigates the risk factors for diabetes in Mexican-Americans [31,32]. An in-home sleep apnea study occurred between 2010 and 2014 using a validated instrument that records finger pulse oximetry, actigraphy, body position, and peripheral arterial tonometry (Itamar-Medical WatchPAT-200) [33]. Seven hundred eighty-two HA individuals were studied, using Affymetrix 6.0 genotyping data. The Western Australian Sleep Health Study (WASH S) is a clinic-based study focused on the epidemiology and genetics of SDB [34]. PSG was obtained from 1508 European-ancestry patients (91% referred for SDB evaluation) from 2006 to 2010 (Compumedics Series E). Genotyping was performed using the Illumina Omni 2.5 array.

Phenotype and covariate definitions
We examined several SDB measures, including specific measures of OSA: AHI (number of apneas plus hypopneas per hour of sleep, with a minimum 3% desaturation per event) and average oxyhemoglobin desaturation per apnea or hypopnea, and measures of SDB severity [14]: average and minimum SpO 2 and the percentage of the night with SpO 2 < 90% (Per90). Apart from WASHS, all sleep data were scored by blinded scorers at one central Sleep Reading Center with high levels of scorer reliability using well-defined procedures [35]. The AHI reflected all events. We did not attempt to disentangle the apneahypopnea index from central versus obstructive sleep apnea events, due to the relatively low prevalence of central sleep apnea (< 2%) in these largely community-based studies [36,37] (some of which are enriched with snorers) and the complexities of classifying mixed events. We adjusted for age, age 2 , sex, age × sex, body mass index (BMI), and BMI 2 due to known age and sex effects, some of which are non-linearly associated with outcomes, and our goal of identifying obesityindependent loci. Age and BMI were obtained at the time of the sleep recording. We adjusted for BMI as over half of the AHI trait heritability is attributable to factors other than obesity as measured by the BMI and our goal was to identify associations with other mechanistic pathways (e.g., ventilatory control) that could indicate novel future targets. Phenotype analyses were pooled within populations to aggregate very rare variants for testing and therefore further adjusted for study. Population assignments were based on self-report, in accordance with other research from TOPMed and other consortia. AsA and EA-identifying individuals with population principal components > 5 standard deviations [38] from applicable 1000 Genomes and Human Genome Diversity Project super-populations were excluded. We used a two-stage procedure to rank-normalize the phenotypes adjusted for covariates [39]. Cryptic relatedness and population substructure were controlled for using linear mixed models. Genomic control was applied to populationspecific results (or cohort-specific imputed genotype results).

WGS and genotyping
Sequence data were derived from the TOPMed Freeze 6a release, jointly called by the TOPMed Informatics Research Center at the University of Michigan (http:// github.com/statgen/topmed_variant_calling). The methodology was described elsewhere [12]. In brief, WGS  [40]). WGS data were merged and normalized; inferred sequence contamination was identified; and SNPs and small indels were detected (structural variants are not currently available). Lower quality variants were excluded using Mendelian consistency checks. Variants were aligned to Build 38 and annotated using snpEff 4.3 t [41]. We excluded variants with < 10× depth or > 5% missingness, leaving 152.7 million polymorphic variants in 7988 individuals with SDB phenotypes. Up to 22,030,888 variants from individuals with sequencing were tested in the GWAS analyses, following filtering for quality control and minor allele frequencies.
Genotype data were imputed using the TOPMed Imputation Server [13] using a Freeze 5b (Build 38) template. Forward strand checks were performed using the Strand database and the Haplotype Reference Consortium imputation preparation script (https://www.well.ox. ac.uk/~wrayner/tools/) and confirmed using Ensembl variant allele checks and internal QC performed on the server. Study-level data were imputed separately. Analyses on variants with r 2 score > 0.5 were therefore performed separately for each study. Up to 22,105,437 variants from individuals with imputed data were tested in the GWAS analyses, following filtering for quality control, imputation r 2 , and minor allele frequencies.

Statistical analyses
Single and grouped variant analyses were performed using EMMAX and MMSKAT, both within the EPAC TS suite (v3.3) [42]. WGS genetic relatedness matrices (GRM) were constructed using autosomal variants (MAF > 0.1%) following a comparison of EPACTS point-wise heritability estimates of the AHI using different minimal MAFs. A grid search identified optimal GRM parameters with imputed data (MAF > 0.5%, r 2 > 0.90) using 929 ARIC individuals with imputation and WGS data. Log 10 P-values using identical association test parameters had a Spearman's ρ correlation of 0.951 between WGS and imputed data. Matrices were constructed separately for each study + population combination (due to potentially differential imputation coverage).
Gene-based group sets considered Ensembl-defined non-pseudogenes expressed in any GTEx v7 tissue. Variants needed to clear a series of frequency, regional, functional class, and presumed functionality score filters in order to test a gene using its most biologically plausible variants. Variants could have a maximum minor allele frequency of 5%. Regions were largely exon-based. We also included variants located within experimentally derived promoter regions and Ensembl-derived Tarbase miRNA binding sites; and regulatory variants located within 1000 bases of a particular gene, including ChIPseq determined transcription factor binding sites (TFBS), and Ensembl-derived CTCF, TFBS, and promoter sites [43][44][45]. Variants from a subset of 19 snpEff gene-based annotation functional classes (e.g., missense or nonsense, but not synonymous mutations) were considered. Finally, group set variants passing these prior filters were additionally filtered for the plausibility of biological function by requiring either a FATHMM-XF score > 0.5 or a CDTS < 1% constrained region score [46,47]. Exonic variants could alternatively have a PrimateAI score > 0.803 or a Havrilla et al. < 1% constrained coding region score [48,49].
Gene-based tests considered variants in WGS-only data. Pooled (across cohort) analyses were performed within each population in order to aggregate information on very rare variants across studies. Combined population results were obtained through meta-analysis of p-values weighted by sample size (due to potentially different MAF spectra driven by population demography). A significance level of p < 4.51 × 10 −8 was used, reflecting a Bonferroni adjustment for all genes tested across all phenotype and population configurations.
A second set-based analysis was designed to query for TFBS annotation enrichment [50]. We performed 250base pair sliding window analyses (to improve power by aggregating additional variants beyond an approximate ChIP-seq peak width of 100 base pairs). We filtered for variants with either a FATHMM-XF score > 0.5 or a CDTS 1% score with no MAF cut-offs and metaanalyzed MMSKAT results across the 4 populations, noting windows with p-values < 0.01. These intervals were tested for enrichment of ChIP-seq coordinates with at least 50% physical overlap for up to 437 transcription factors using ReMap 2018 v1.2 [51].
Single-variant EMMAX tests examined common variants (MAF > 0.5%). Meta-analysis across populations (and imputed genotype studies) used METAL with genomic control [52]. We performed bidirectional discovery and replication using the WGS and imputed samples (noting the high genomic resolution in the WGS samples and the higher sample size in the imputed data). We report results including at least 1000 individuals in discovery analyses, discovery association p-values < 1 × 10 −5 and replication association p-values < 0.05. Therefore, no population-specific discovery analyses of Asian-Americans were performed. Multi-ethnic analyses included a minimum of two populations where a variant cleared minimum MAF and imputation quality (for chip-based results) criteria. Significance was defined as p < 1 × 10 −8 in joint analyses, reflecting adjustment for five correlated phenotypes (Additional file 1: Table S3). We performed MetaXcan imputed GTEx gene expression analyses using joint EA results in selected tissues relevant to SDB and GIGSEA pathway analyses of MetaXcan output in whole blood (to maximize power), with empirical p-values incorporating 10,000 permutations [53,54]. Bioinformatics annotations of single-variant results (Additional file 1: Table S7) include significant eQTL associations from GTEx v7, and overlapping promoter and enhancer coordinates derived from Roadmap Epigenomics, BLUEPRINT, and Vermunt et al. brain tissues (enhancers only) [55][56][57][58]. Lookups of potentially druggable genes as defined within DGIdb, a database of 56,000 drug-gene interactions from over 30 literature sources, were performed using the GeneCards suite [59,60].

Study sample
A study overview is provided in Additional file 2: Figure  S1. Tables 1 and 2 provide a summary of the study samples and SDB traits analyzed using WGS and imputed genotypes, respectively. In total, there were 21,244 individuals (1942 AAs, 229 AsAs, 8341 EAs, and 10732 HAs). Median AHI levels ranged from mildly to moderately elevated, reflecting the age range and sex distribution of each cohort. Pairwise correlations of phenotypes and covariates are provided in Additional file 1: Table S3.

MetaXcan imputed gene expression and GIGSEA pathway analyses
We used joint WGS and imputed EA results to impute associations with gene expression levels using a MetaXcan framework for six tissues (subcutaneous and visceral omentum adipose, lung, monocytes, skeletal muscle, and whole blood). No individual tests reached Bonferroni significance (p < 2.60 × 10 −7 ; Additional file 1: Table S9). Genes that were observed in the top 10 results across the varied analyses (Additional file 1: Table S10) included ZNF83 (15 instances) and CHRNE (13 instances).
Whole blood MetaXcan results (with the largest sample size) were further evaluated in GIGSEA-based pathway analyses. KEGG pathway results are shown in Additional file 1: Table S11. The most significantly associated pathway was KEGG_STEROID_HORMONE_BIOSYNTHESIS (average SpO 2 empirical p-value = 7.00 × 10 −4 ). KEGG_ RIG_I_LIKE_RECEPTOR_SIGNALING_PATHWAY was observed in the top 10 results for four of the five phenotypes. Gene-centric transcription factor binding site (TFBS) enrichment analysis results are presented in Additional file 1: Table S12. V$PEA3_Q6 (ETV4) was the most significantly associated TFBS (average desaturation empirical p-value = 3.00 × 10 −4 ) and was the strongest association for AHI and minimum SpO 2 (empirical p-values 0.002 and 0.001, respectively). The most significant miRNA binding site enrichment analysis association was GCATTTG,MIR-105 (average SpO 2 p = 0.002; Additional file 1: Table S13). AGGCACT,MIR-515-3P (the strongest AHI association, p = 0.009) was observed in the top ten results for four phenotypes.

ChIP-seq transcription factor binding site interval enrichment
We performed a sliding window analysis to examine enriched intervals containing ChIP-seq derived coordinates for up to 437 transcription factors ( Table 6, Additional file 1: Table S14). FOXP2 TFBS were consistently the most enriched for all phenotypes. Other notable transcription factors in the top 5 included EGR1, KDM4B, KDM6B, and TP63. KDM4B and KDM6B are druggable [59,60]. Leading sliding window results are provided in Additional file 1: Table S15.

Discussion
Sleep-disordered breathing is associated with increased risk of a wide range of disorders, including cardiometabolic disease, cancer, cognitive impairment, and interstitial lung diseases, as well as premature mortality [2,61]. Treatment options, however, are limited by a lack of knowledge of molecular pathways, including those that may be "druggable." Recent analyses of SDB traits have focused on common variants and identified several preliminary genomelevel significant associations [9][10][11], but did not address gene-based or rare-variant effects. Ten studies and over 21, 000 individuals of multiple ancestries with WGS data at unprecedented resolution from the NHLBI TOPMed program combined with densely imputed data from other sources contributed to these results. We identified several variant, gene-based, and pathway-level associations. Analyses adjusted for obesity, a major SDB risk factor, identified loci and genes implicated in pulmonary, inflammatory, and craniofacial pathways. Some associations were populationspecific, while others were sex-specific, consistent with population differences and strong sex differences for SDB [20,62]. Notably, across multiple ancestral groups, we identified a set-based rare-variant association (p = 3.48 × 10 −8 ) on chromosome X with ARMCX3. Lead MMSKAT gene-based results meta-analyzed across populations within one order of magnitude of significance (p < 4.51 × 10 −8 ) are shown. Populationspecific information for each gene is displayed in the latter columns for AA, AsA, EA, and HA, respectively. Individual populations varied in the number of polymorphic variants available for testing (e.g., due to singletons or excessively common variants). ARMCX3-AS1 is a RNA gene that is anti-sense to the proteincoding ARMCX3 gene. Full results for genes with p < 0.01, including Ensembl-derived gene biotypes and descriptions, are provided in Additional file 1: Table S4. A list of individual variants comprising each gene is provided in Additional file 1:  Lead MMSKAT gene-based population-specific associations within one order of magnitude of significance (p < 4.51 × 10 −8 ) are shown. The Variants column indicates the number of filtered polymorphic variants with minor allele frequency < 5% available for testing, a portion of which were singletons. *Druggable gene [59,60]. Full results for genes with p < 0.01, including descriptions, are provided in Additional file 1: Table S5. A list of individual variants comprising each gene is provided in Additional file 1:

Gene-based results
Across multiple populations, ARMCX3 (ALEX3) and the RNA anti-sense gene ARMCX3-AS1 were associated with apnea-hypopnea triggered intermittent hypoxia. ARMCX3 regulates mitochondrial aggregation and trafficking in multiple tissues and facilitates neuronal survival and axon regeneration [63][64][65]. Wnt signaling regulates reactive oxygen species (ROS) generation and ARMCX3-associated mitochondrial aggregation [64,66]. Potential mechanisms for further study include sensitized carotid body chemoreflexes, interaction with inflammatory mechanisms, and neuronal dysfunction within respiratory centers. Sleep apnea and reduced ventilatory drive are enriched in individuals with a primary mitochondrial disorder [67]. Mitochondria are an important source of ROS, which modulate the acute hypoxic ventilatory response. Mitochondria impact HIF1A signaling and may contribute to oxygen sensing [68,69].
ROS are required for intermittent hypoxia-induced respiratory long-term facilitation [70]. These effects may mitigate the level of hypoxia resulting from recurrent apneas, or conversely, lead to ventilatory instability, promoting apnea occurrence. Mitochondrial ROS also activate the NLRP3 inflammasome in multiple pulmonary diseases, consistent with an inflammation model that includes our IL18-pathway and HK1 results, ROSrelated proinflammatory responses to lung capillary pressure, and evidence of alveolar epithelial injury/SDB interactions [10,69,[71][72][73]. Our findings suggest value in investigating the mechanisms by which ARMCX3 predisposes to SDB, and whether these associations are mediated by neuronal dysfunction and/or ROS and carotid body sensitization, and interact with the inflammasome. Additional genes were significantly associated in population-specific analyses, including the mitochondrial ribosomal gene MRPS33. Mitoribosomes are responsible Two-hundred-fifty-base pair sliding window coordinates with association p < 0.01 were queried for interval enrichment of ChIP-seq-derived transcription factor binding sites using the ReMap annotation tool. ChIP-seq coordinates were required to have >50% overlap with a sliding window interval. ReMap-derived expected overlaps are obtained from the equivalent number of similarly sized random regions. E-value indicates the expected value, with a higher logtransformed value indicating greater enrichment. Full results are provided in Additional file 1: Table S14 for the expression of the 13 essential components of the oxidative phosphorylation system, and a majority of the small subunit proteins have been implicated in disease [74]. The expression of several small and large subunit proteins are altered in a hypoxic environment [75]. MRPS33 expression varies with oxygen treatment in COPD [76].

Single-variant results
We identified four common frequency associated loci, including multiple-population associations with the IL18RAP region. The IL18RAP region has been associated with minimum SpO 2 [10], and here we further identify an association with average event desaturation, highlighting a role in an OSA-specific trait. Multiple variants in this region are also GTEx eQTL variants for both interleukin-18 receptor subunits IL18RAP and IL18R1 (Additional file 1: Table S7) and experimental studies support a role for IL18 signaling in mediating this association, possibly through effects of pulmonary inflammation on gas exchange (reviewed in [10]). We identified three population-specific loci, including two novel associations in individuals of European ancestry ( Figs. 1 and 2). Sixty-five variants in the NRG1 region were associated with the AHI (p < 1.0 × 10 −8 , Additional file 1: Table S7). This region was suggestively associated with sleep apnea in a Korean population [77]; however, the lead signals appear to be independent (rs10097555 Korean p = 2.6 × 10 −6 , EA p = 0.91). NRG1 is associated with lung development and acute lung injury and mediates inflammasome-induced alveolar cell permeability [78][79][80]. NRG1 promotes accumulation of HIF1A and has increased expression in vascular smooth muscle cells following exposure to intermittent hypoxia [81,82]. The lead SLC45A2 region variant rs28777 (average SpO 2 p = 8.08 × 10 −10 ) has been associated with multiple traits and is in a splicing regulatory element with extreme population differentiation [83]. An association in the ATP2B4 region with average SpO 2 in HAs [9] has been extended to a second hypoxemia trait at the same variant (Per90 p = 3.31 × 10 −10 ). This gene is the main cellular membrane calcium pump in erythrocytes and also regulates vascular tone [84,85].

Pathway analyses
Several gene pathways were identified in EA individuals using imputed gene expression in whole blood (Additional file 1: Table S11). KEGG_RIG_I_LIKE_RECEPTOR_SIG-NALING_PATHWAY (retinoic acid-inducible gene Ilike) was the most commonly observed, occurring in the top 10 results for 4 of the 5 phenotypes. This pathway initiates the immune response to RNA virus infection [86], consistent with a role for inflammation at the NRG1 and IL18RAP loci. Steroid hormone biosynthesis (the most significantly associated pathway), PPAR signaling, and metabolism (via "starch and sucrose metabolism") suggest the importance of biological pathways modulating energy homeostasis and balance and metabolic function [87]. In the gene-centric GIGSEA TFBS analysis, V$PEA3_Q6 (ETV4) was the lead association for three phenotypes. ETV4 influences branching in the developing lung and regulates hypoxia-inducible factor signaling [88,89], a major mechanism influencing ventilatory control.

Transcription factor binding site enrichment
Several transcription factors were identified through interval enrichment of observed TFBS across the genome ( Table  6). FOXP2 was consistently the most enriched transcription factor and is known to regulate gene expression in epithelial lung tissue and response to lung injury through an inflammatory mechanism [90,91]. FOXP2 is also expressed in brainstem respiratory areas including the pre-Bötzinger complex (which is essential for respiratory rhythmogenesis) and impacts airway morphology [92,93]. Two lysine demethylases (KDM4B and KDM6B) were also identified. KDM6B (JMJD3) is required for a functional pre-Bötzinger complex [94,95] and reduced KDM6B protein expression was reported in hypoxic OSA patients [96]. Kdm6b also plays roles in immune function and lung development [97][98][99]. Drosophila Kdm4b knock-outs have increased sleep [100]. KDM4B (JMJD2B) and KDM6B are both members of the JmjC protein domain family and are regulated by HIF1A, require oxygen as a cofactor, and act as oxygen sensors for chromatin in hypoxia [101,102]. EGR1 mediates hypoxia-induced pulmonary fibrosis [103]. TP63 is associated with cleft palate in Tp63 deficient mice, which is associated with an increased prevalence of OSA [104,105], suggesting that its relationship to OSA may be through pathways influencing craniofacial development. Among the leading 250-base pair sliding window results (Additional file 1: Table S15), 4:105708751-105709001 (Per90 HA p = 2.72 × 10 −9 ) is of note due to regional associations with lung function and expression in the human lung [106].

Strengths and weaknesses
This study is the first genome-wide analysis of objectively measured SDB traits using deep sequencing. Together with improved imputation quality, the TOPMed resource has enabled unprecedented genetic resolution. We examined clinically relevant phenotypes measured using rigorous methodology [2,14]. We analyzed data from 10 studies of individuals from four population groups that used different ascertainment strategies, which may potentially improve the generalization of our results. While this analysis is among the largest performed for SDB traits to date, our moderate sample size has lower power to detect weaker associations, and data were not available to replicate these first rare-variant associations. We did not specifically study the central apnea-hypopnea index due to the relatively low prevalence of central sleep apnea (< 2%) in these largely community-based studies [36,37]. While there are multiple lines of evidence in the literature to support our findings, additional experimental follow-up analyses are required.

Conclusions
We have identified the first rare-variant and additional common-variant associations at genome-level significance with objectively measured SDB traits in humans. The results point to biologically relevant pathways for further study, including a novel X-linked association (ARMCX3), and a number of associations in genes that modulate lung development, inflammation, respiratory rhythmogenesis, and HIF1A-mediated hypoxic-response pathways. These associations will motivate future sample collection and follow-up in cell-line and animal validation studies, with potential therapeutic benefit for sleep-disordered breathing and related comorbidities.

Acknowledgements
We acknowledge our TOPMed Consortium and TOPMed Sleep Traits Working Group collaborators, who are listed in Additional file 1: Tables S1 and S2. The authors wish to thank the participants and study staff of all of our cohorts for their important contributions. We gratefully acknowledge the studies and participants who provided biological samples and data for TOPMed. The authors thank the staff and participants of the ARIC study for their important contributions. A full list of principal CHS investigators and institutions can be found at CHS-NHLBI.org. The Framingham Heart Study thanks the study participants and the multitude of investigators who over its 70-year history continue to contribute so much to further our knowledge of heart, lung, blood, and sleep disorders and associated traits. The authors thank the staff and participants of HCHS/SOL for their important contributions. The views expressed in this manuscript are those of the authors and do not necessarily represent the views of the National Heart, Lung, and Blood Institute; the National Institutes of Health; or the US Department of Health and Human Services. This manuscript was not approved by the HCHS/SOL publications committee. Investigator's website-http://www.cscc. unc.edu/hchs/. The authors also wish to thank the staff and participants of the JHS. We thank the field staff in Starr County for their careful collection of these data and are especially grateful to the participants who so graciously cooperated and gave of their time. had full access to the study data and take responsibility for the integrity of the data and accuracy of analyses.

Funding
Whole-genome sequencing (WGS) for the Trans-Omics in Precision Medicine (TOPMed) program was supported by the National Heart, Lung, and Blood Institute (NHLBI). WGS for "NHLBI TOPMed: Trans-Omics for Precision Medicine Whole Genome Sequencing Project: ARIC" (phs001211.v1.p1) was performed at Baylor College of Medicine Human Genome Sequencing Center (HHSN268201500015C and 3U54HG003273-12S2) and the Broad Institute of MIT and Harvard (3R01HL092577-06S1). WGS for "NHLBI TOPMed: The Cleveland Family Study (WGS)" (phs000954.v2.p1) was performed at the University of Washington Northwest Genomics Center (3R01HL098433-05S1). WGS for "NHLBI TOPMed: Cardiovascular Health Study" (phs001368.v1.p1) was performed at Baylor College of Medicine Human Genome Sequencing Center (HHSN268201500015C). WGS for "NHLBI TOPMed: Whole Genome Sequencing and Related Phenotypes in the Framingham Heart Study" (phs000974.v3.p2) was performed at the Broad Institute of MIT and Harvard (3R01HL092577-06S1). WGS for "NHLBI TOPMed: Hispanic Community Health Study/Study of Latinos (HCHS/SOL)" (phs001395) was performed at the Baylor College of Medicine Human Genome Sequencing Center (HHSN268201500015C and 3U54HG003273-12S2). WGS for "NHLBI TOPMed: The Jackson Heart Study" (phs000964.v3.p1) was performed at the University of Washington Northwest Genomics Center (HHSN268201100037C). WGS for "NHLBI TOPMed: NHLBI TOPMed: MESA" (phs001416.v1.p1) was performed at the Broad Institute of MIT and Harvard (3U54HG003067-13S1). Centralized read mapping and genotype calling, along with variant quality metrics and filtering were provided by the TOPMed Informatics Research Center (3R01HL-117626-02S1). Phenotype harmonization, data management, sample-identity QC, and general study coordination were provided by the TOPMed Data Coordinating Center (3R01HL-120393-02S1). The Genome Sequencing Program (GSP) was funded by the National Human Genome Research Institute (NHGRI); the National Heart, Lung, and Blood Institute (NHLBI); and the National Eye Institute (NEI). The GSP Coordinating Center (U24 HG008956) contributed to cross-program scientific initiatives and provided logistical and general study coordination. The Centers for Common Disease Genomics (CCDG) program was supported by NHGRI and NHLBI, and CCDG-funded whole-genome sequencing of the ARIC and HCHS/SOL studies was performed at the Baylor College of Medicine Human Genome Sequencing Center (UM1 HG008898 and R01HL059367