Identification of novel rare sequence variation underlying heritable pulmonary arterial hypertension

Identification of novel rare sequence variation underlying heritable pulmonary arterial hypertension Stefan Gräf* [1,2,3], Matthias Haimel [1,2,3], Marta Bleda [1], Charaka Hadinnapola [1], Laura Southgate [4,5], Wei Li [1], Joshua Hodgson [1], Bin Liu [1], Richard M. Salmon [1], Mark Southwood [6], Rajiv D. Machado [7], UK PAH Cohort Study Consortium, NIHR BioResource – Rare Disease Consortium, Jennifer M. Martin [1,2,3], Carmen M. Treacy [1,6], Katherine Yates [1,2,3], Louise C. Daugherty [2,3], Olga Shamardina [2,3], Deborah Whitehorn [2,3], Simon Holden [8], Micheala Aldred [9], Harm J. Bogaard [10], Colin Church [11], Gerry Coghlan [12], Robin Condliffe [13], Paul A. Corris [14], Cesare Danesino [15,16], Mélanie Eyries [17], Henning Gall [18], Stefano Ghio [16], Hossein-Ardeschir Ghofrani [18,19], J. Simon R. Gibbs [20], Barbara Girerd [21], Arjan C. Houweling [10], Luke Howard [19], Marc Humbert [21], David G. Kiely [13], Gabor Kovacs [22,23], Robert V. MacKenzie Ross [24], Shahin Moledina [25], David Montani [21], Michael Newnham [1], Andrea Olschewski [22], Horst Olschewski [22,23], Andrew J. Peacock [11], Joanna Pepke-Zaba [6], Inga Prokopenko [19], Christopher J. Rhodes [19], Laura Scelsi [16], Werner Seeger [18], Florent Soubrier [17], Dan F. Stein [1], Jay Suntharalingam [24], Emilia Swietlik [1], Mark R. Toshner [1], Anton Vonk Noordegraaf [10], David A. van Heel [26], Quinten Waisfisz [10], John Wharton [19], Stephen J. Wort [27,19], Willem H. Ouwehand [2,3], Nicole Soranzo [2,28], Allan Lawrie [29], Paul D. Upton [1], Martin R. Wilkins [19], Richard C. Trembath [5], Nicholas W. Morrell* [1,3]


Introduction
Idiopathic and heritable pulmonary arterial hypertension (PAH) are rare disorders characterised by occlusion of arterioles in the lung 1 , leading to marked increases in pulmonary vascular resistance 2 . Life expectancy from diagnosis averages 3-5 years 3 , with death ensuing from failure of the right ventricle.
Mutations in the gene encoding the bone morphogenetic protein type 2 receptor (BMPR2), a receptor for the transforming growth factor-beta (TGF-β) superfamily 4, 5 account for over 80% of families with PAH, and approximately 20% of sporadic cases 6 . Mutations have been identified in genes encoding other components of the TGF-β/bone morphogenetic protein (BMP) signalling pathways, including ACVRL1 7 and ENG 8 . On endothelial cells, BMPR2 and ACVRL1 form a signaling complex, utilizing ENG as a co-receptor. Case reports of rare sequence variation in the BMP signalling intermediaries, SMAD1, SMAD4 and SMAD9 9, 10 , provide compelling evidence for a central role of dysregulated BMP signalling in PAH pathogenesis.
Analysis of coding variation in BMPR2-negative kindreds revealed heterozygous mutations in genes not directly impacting on the TGF-β/BMP pathway, including CAV1 11 , and the potassium channel, KCNK3 12 . Deletions and loss of function mutations in TBX4, an essential regulator of embryonic development, were identified in childhood onset PAH 13 . A clinically and pathologically distinct form of PAH, known as pulmonary veno-occlusive disease or pulmonary capillary haemangiomatosis (PVOD/PCH), was shown recently to be caused by biallelic recessive mutations in EIF2AK4 14, 15 , a kinase in the integrated stress response.
The purpose of the present study was to identify additional rare sequence variation contributing to the genetic architecture of PAH, and to assess the relative contribution of rare variants in genes implicated in prior studies. A major finding is that rare likely causal heterozygous variants in several previously unidentified genes (ATP13A3, AQP1 and SOX17) were significantly overrepresented in the PAH cohort, and we provide independent validation for GDF2 as a causal gene.

Description of the PAH cohort
In total, 1048 PAH cases (1038 index cases and 10 related affected individuals) were recruited for WGS. Of these, 908 (86.7%) were diagnosed with idiopathic PAH, 58 (5.5%) gave a family history of PAH and 60 (5.7%) gave a history of drug exposure associated with PAH 16 . Twenty two cases (2.1%) held a clinical diagnosis of PVOD/PCH (Figure 1a). Demographic and clinical characteristics of the PAH cohort are provided in Supplementary  Table 1. An additional UK family was recruited separately for novel gene identification studies. Briefly, the proband was diagnosed at 12 years with a persistent ductus arteriosus and elevated pulmonary arterial pressure. Explant lung histology following heart-lung transplantation revealed the presence of plexiform lesions. Two of the proband's offspring were also diagnosed with childhood-onset PAH, one of which had an atrial septal defect.
The proband's parents, siblings and a third child showed no evidence of cardiovascular disease.

Pathogenic variants in previously reported PAH disease genes
Our filtering strategy detected rare deleterious variation in previously reported PAH genes in 19.9% of the PAH cohort. For BMPR2, rare heterozygous mutations were identified in 160 of 1048 cases (15.3%). The frequency of BMPR2 mutations in familial PAH was 75.9%, in sporadic cases 12.2%, and 8.3% in anorexigen-exposed PAH cases. Forty-eight percent of BMPR2 mutations were reported previously 17 , and the remainder were newly identified in this study. Fourteen percent of BMPR2 mutations resulted in the deletion of larger proteincoding regions ranging from 5 kb to 3.8 Mb in size. Supplementary Table 2 provides the breakdown of BMPR2 SNVs and indels, and the larger deletions are shown in Figure 2a-c with a detailed summary in Supplementary Table 3.
In a case-control comparison of the frequencies of deleterious variants confined to the previously reported PAH genes, we observed significant overrepresentation of rare variants in BMPR2, TBX4, ACVRL1 and biallelic variants in EIF2AK4 only (P < 0.05) (Supplementary Table 6).

Identification of novel PAH disease genes
The strategy to identify novel causative genes in PAH employed a series of case-control analyses (Figure 1b). To identify signals that might be masked by variants in previously reported PAH genes, we excluded subjects with rare variants and deletions in BMPR2, EIF2AK4, ENG, ACVRL1, TBX4, SMAD9 and KCNK3. A genome-wide comparison of protein truncating variants (PTVs), representative of high impact variants, identified a higher frequency of PTVs in ATP13A3 (6 cases) (P adj = 0.0346). Moreover, we identified additional PTVs in several putative PAH genes, including EVI5 (5 cases, 1 control) and KDR (4 cases, 0 controls) (Figure 3a), that require further validation to evaluate their contribution to PAH pathogenesis (Supplementary Table 7).
We next analysed rare missense variants overrepresented in the PAH cohort, again excluding subjects with variants in the previously reported PAH genes. This revealed significant overrepresentation of rare variants in GDF2 after correction for multiple testing (P adj = 0.0023), followed by AQP1 (Figure 3b and Supplementary Table 8). Next, in a combined analysis of rare missense variants and PTV, only GDF2 remained significant (P = 0.001). Rare variants in additional putative genes occurred at higher frequency in cases compared to controls, including AQP1, ALPPL2, ATP13A3, OR8U1, IFT74, FLNA, SOX17, ATP13A5, C3orf20 and PIWIL1 (uncorrected P value < 0.0005), but were not significant after correction for multiple testing (Figure 3c and Supplementary Table 9).
In order to increase power to detect rare associations, we deployed SKAT-O on filtered rare PTVs and missense variants. Excluding previously reported genes, this analysis revealed an association with rare variants in AQP1 (P adj = 4.28x10 -6 ) and SOX17 (P adj = 6.7x10 -5 ) ( Figure  3d). AQP1 and SOX17 were both also nominally significant in the combined burden tests, described above. Association was also found with rare variants in MFRP (P adj = 1.3x10 -5 ). However, we consider MFRP a false-positive finding for reasons given in the Discussion. Supplementary Table 10 shows the top 50 most significant genes identified by SKAT-O, providing further candidates to be evaluated in future studies. Details of rare variants in novel PAH genes (GDF2, ATP13A3, AQP1, SOX17) identified in cases are provided in Supplementary Table 11.
Notably, a genome-wide assessment of larger structural variation did not identify any additional large deletions after exclusion of subjects harbouring deletions in BMPR2 ( Figure  2d-e).
The proportion of PAH cases with mutations in the new genes was 3.5%. The clinical characteristics of PAH cases with mutations in these genes are provided in Supplementary  Table 5b. Of note, cases with mutations in SOX17 and AQP1 were significantly younger at diagnosis (32.8 ± 16.2 years [P = 0.002] and 36.9 ± 14.3 years [P = 0.013], respectively) compared to cases with no mutations in the previously established genes (51.7 ± 16.6 years).

Independent validation and familial segregation analysis
To provide further validation of the potentially causal role of mutations in the new genes identified, we examined whole-exome data from an independent UK family with three affected individuals across two generations. Microsatellite genotyping across chromosome 2q33 had previously demonstrated non-sharing of haplotypes in affected individuals, consistent with exclusion of linkage to the BMPR2 locus. No pathogenic variants were identified in the protein-coding regions of the BMPR2 gene or other TGF-β pathway genes. Analysis of exome sequence data from individual II-1 identified a novel heterozygous c.411C>G (p.Y137*) PTV in the SOX17 gene. Segregation analysis in the extended family demonstrated that the mutation had arisen de novo in the affected father (II-1) and was transmitted to the affected offspring (III-1). All unaffected family members were confirmed as wild-type ( Figure 4a).
Three HPAH subjects harbouring rare variants in AQP1, identified in the NIHR BR-RD WGS study, were also selected for familial co-segregation analysis . No pathogenic variants in any of the previously reported genes were identified in these families. The first pedigree comprised three affected individuals across two generations. Sanger sequencing confirmed the presence of the heterozygous AQP1 c.583C>T (p.R195W) missense variant in the proband (E011942), the affected father (E011942.f) and the healthy younger paternal uncle (E011942.u1). An additional unaffected uncle did not carry the AQP1 variant. These results indicate likely incomplete penetrance in the unaffected carrier, as observed in BMPR2 families 20 . No additional clinical information was available for the deceased grandparents ( Figure 4b). The remaining two families comprised affected parent-offspring individuals. By Sanger sequencing we independently confirmed a heterozygous AQP1 c.527T>A (p.Val176Glu) missense variant in proband (E012415) and his affected father (Figure 4c), as well as a heterozygous AQP1 c.583C>T (p.R195W) missense variant in proband (E010634) and her affected father (Figure 4d). These results highlight recurrent AQP1 variation across unrelated families and demonstrate co-segregation with the phenotype.

Predicted functional impact of variants in novel PAH genes
To evaluate the potential functional impact of rare variants identified in the likely causative new genes we performed structural analysis of GDF2, ATP13A3, AQP1, and SOX17. In addition we undertook a functional analysis of the GDF2 variants identified.
Heterozygous mutations in GDF2 exclusive to PAH cases comprised 1 frameshift variant and 7 missense variants. GDF2 encodes growth and differentiation factor 2, also known as bone morphogenetic protein 9 (BMP9), the major circulating ligand for the endothelial BMPR2/ACVRL1 receptor complex 21 . Amino acid substitutions were assessed against the published crystal structure 22 of the prodomain bound form of GDF2 ( Figure 5). Variants clustered at the interface between the prodomain and growth factor domain. Since the prodomain is important for the processing of GDF2, it is likely that amino acid substitutions reduce the stability of the prodomain-growth factor interface. In keeping with these predictions, HEK293T cells transfected with GDF2 variants exclusive to PAH cases, demonstrated reduced secretion of mature GDF2 into the cell supernatants (Figure 5d), compared with wild type GDF2.
We identified 3 heterozygous frameshift variants, 2 stop gained, 2 splice region variants in ATP13A3, which are predicted to lead to loss of ATPase catalytic activity ( Figure 6a). In addition, we identified 4 heterozygous likely pathogenic missense variants in PAH cases, two near the conserved ATPase catalytic site and predicted to destabilise the conformation of the catalytic domain (Figure 6b-d). The distribution of variants ( Figure 6a) suggests that these mutations impact critically on the function of the protein.
The majority of rare variants identified in AQP1, which encodes aquaporin 1, are situated within the critical water channel (Figure 7). In particular the p.Arg195Trp variant, identified in 5 PAH cases, locates at the hydrophilic face of the pore. This arginine at position 195 helps define the constriction region of the AQP1 pore structure and is conserved across the water specific aquaporins 26 . Rare variants in SOX17, included 4 nonsense variants (including the PTV identified in the additional UK family) predicted to lead to loss of the beta-catenin binding region, and 6 missense variants predicted to disrupt interactions with Oct4 and betacatenin 27,28 (Figure 8).
GDF2 is known to be secreted from the liver, but the cellular localization of proteins encoded by the other novel genes is less well characterised. Thus we employed immunohistochemistry to examine localization in the normal and hypertensive human pulmonary vasculature. Figure 9 shows that AQP1, ATP13A3 and SOX17 are predominantly localised to the pulmonary endothelium in normal human lung and to endothelial cells within plexiform lesions of patients with idiopathic PAH. In addition, we determined the relative mRNA expression levels of AQP1, ATP13A3 and SOX17 in primary cultures of pulmonary artery smooth muscle cells (PASMCs), pulmonary artery endothelial cells (PAECs) and blood outgrowth endothelial cells (BOECs) 31 . AQP1 was expressed in PASMCs and endothelial cells, with a trend towards higher levels in PASMCs (Figure 10a). ATP13A3 was highly expressed in both cell types (Figure 10b), whereas SOX17 was almost exclusively expressed in endothelial cells (Figure 10c). Although AQP1 and SOX17 are known to play roles in endothelial function, the function of ATP13A3 in vascular cells is entirely unknown. Thus, we determined the impact of ATP13A3 knockdown on proliferation and apoptosis of BOECs. Loss of ATP13A3 led to marked inhibition of serum-stimulated proliferation of BOECs, and increased apoptosis in serum-deprived conditions (Figure 10d-f).

Discussion
We report the most comprehensive analysis to date of rare genetic variation in a large cohort of index cases with idiopathic and heritable forms of PAH. Whilst we utilised WGS, the main goal was the identification of novel rare causal variation underlying PAH in the protein coding sequence. The approach involved a rigorous case-control comparison using a tiered search for variants. First, we searched for high impact PTVs overrepresented in cases, having excluded previously established PAH genes. This revealed PTVs in ATP13A3, a poorly characterised P-type ATPase of the P5 subfamily 32 . There is little information regarding the function of the ATPase, ATP13A3, which appears widely expressed in mouse tissues 32 .
Although, the precise substrate specificity is unknown, ATP13A3 plays a role in polyamine transport 33 . Based on available RNA sequencing data, ATP13A3 is highly expressed in human pulmonary vascular cells and cardiac tissue (https://www.encodeproject.org). We confirmed that ATP13A3 mRNA is expressed in primary cultured pulmonary artery smooth muscle cells and endothelial cells, and provide preliminary data that loss of ATP13A3 inhibits proliferation and increases apoptosis of endothelial cells. These findings are consistent with the widely accepted paradigm that endothelial apoptosis is a major trigger for the initiation of PAH 34, 35 . It will be of considerable interest to determine the role of ATP13A3 in vascular cells and whether it is functionally associated with BMP signalling, or represents a distinct therapeutic target in PAH.
Analysis of missense variation, and a combined analysis of all predicted deleterious variation, revealed that mutation at the GDF2 gene is also significant determinant of predisposition to PAH. Of the new genes identified, GDF2 provides further evidence for the central role of the BMP signalling pathway in PAH. GDF2 encodes the major circulating ligand for the endothelial BMPR2/ACVRL1 receptor complex 21 . Taken together, the genetic findings suggest that a deficiency in GDF2/BMPR2/ACVRL1 signalling in pulmonary artery endothelial cells is critical in PAH pathobiology. The majority of GDF2 variants detected in our adult-onset PAH cohort were heterozygous missense variants, in contrast to a previous case report of childhood onset PAH due to a homozygous nonsense mutation 36 . The finding of causal GDF2 variants in PAH cases, associated with reduced production of GDF2 from cells, provides further support for investigating replacement of this factor as a therapeutic strategy in PAH 37 .
To maximise the assessment of rare variation in a case-control study design, we deployed the SKAT-O test. This approach revealed a significant association of rare variation in the aquaporin gene, AQP1, and the transcription factor encoded by SOX17. Of note, both AQP1 and SOX17 were within the top 8 ranked genes in our combined PTV and missense burden test analysis (Supplementary Table 9), providing further confidence in their causative contribution to PAH.
Aquaporin 1 belongs to a family of membrane proteins that function as channels facilitating water transport in response to osmotic gradients 26 , and AQP1 is known to promote endothelial cell migration and angiogenesis 38 . Thus, approaches that maintain or restore pulmonary endothelial function could offer new therapeutic directions in PAH. Conversely, AQP1 inhibition in pulmonary artery smooth muscle cells ameliorated hypoxia-induced pulmonary hypertension in mice 39 , suggesting that further studies are required to determine the key cell type impacted by AQP1 mutations in human PAH, and the functional impact of these AQP1 variants on water transport. The demonstration of familial segregation of AQP1 variants with PAH provides further support for the potentially causal role of these mutations in disease. However, we also identified an unaffected AQP1 variant carrier consistent with reduced penetrance, which is well described for other PAH genes, including BMPR2.
Although functional studies are required to confirm the mechanisms by which mutations in SOX17 cause PAH, this finding provides additional support for the vascular endothelium as the major initiating cell type in this disorder. SOX17 encodes the SRY-box containing transcription factor 17, which plays a fundamental role in angiogenesis 40 and arteriovenous differentiation 41 . Moreover, conditional deletion of SOX17 in mesenchymal progenitors leads to impaired formation of lung microvessels 42 . The demonstration of familial segregation of the SOX17 p.Y137* PTV with early onset PAH provides additional evidence for a causal role for these variants in PAH. The co-existence of a patent ductus arteriosus in the index case and an atrial septal defect (ASD) in one of the affected offspring is of interest and suggests an association with congenital heart disease. Small ASDs are not uncommon in idiopathic PAH, and a more detailed clinical phenotyping of SOX17 mutation carriers will be required to determine whether the presence of ASDs and other congenital heart abnormalities are more common in carriers of these mutations.
Whilst the SKAT-O analysis also provided support for the MFRP gene, recessive bi-allelic mutations in MFRP cause retinal degeneration and posterior microphthalmos 43 . The expression of MFRP transcripts is largely confined to the central nervous system 44 and the majority of variants were present in the Genome Aggregation Database (GnomAD, http://gnomad.broadinstitute.org). On the basis of these considerations, variants in MFRP are unlikely to contribute to PAH aetiology.
This analysis provides new insights on the frequency and validity of previously reported genes in PAH. We confirmed that mutations in BMPR2 are the most common genetic cause and validated rare causal variants in ACVRL1, ENG, SMAD9, TBX4, KCNK3 and EIF2AK4. Although our findings question the validity of CAV1, SMAD1 and SMAD4 as causal genes, previous reports might represent private mutations occurring in very rare families. The use of WGS in this study allowed closer interrogation of larger deletions around the BMPR2 locus than has been possible previously. Nevertheless, additional analyses are required to determine the full impact of structural variation (inversions, duplications, smaller deletions) at this and other loci.
The non-PAH cohort used in the case-control comparisons for this study comprised individuals, or relatives of individuals, with other rare diseases recruited to the NIHR Bioresource for Rare Diseases (NIHR BR-RD) in the UK (see Methods). In general, for very rare causal variants, the comparison between PAH cases and non-PAH rare disease controls should not reduce our ability to detect overrepresentation of rare variants in a particular gene in the PAH cohort, if mutations in that gene are specific to PAH. However, if rare variants in a gene were responsible for more than one phenotype, it is possible that this would reduce the power to detect overrepresentation in the PAH cohort. For example, if mutations occurred in different functional domains of the expressed protein, this might lead to PAH if mutations affected one domain, but other phenotypes if they affected another domain. Overcoming this potential limitation will require additional analysis of the functional impact of variants and their distribution within a gene, and more detailed information on the phenotypes of subjects in the non-PAH group.
Taken together, this study identifies rare sequence variation in new genes underlying heritable forms of PAH, and provides a unique resource for future large-scale discovery efforts in this disorder. Mutations in previously established genes accounted for 19.9% of PAH cases. Including new genes identified in this study (GDF2, ATP13A3, AQP1, SOX17), the total proportion of cases explained by mutations increased to 23.5%. It is likely that independent confirmation of the expanded list of putative genes identified in this study will increase further the proportion of cases explained by mutations, but this will require larger international collaborations. The results suggest that the genetic architecture of PAH, beyond mutations in BMPR2, is characterised by substantial genetic heterogeneity and consists of rare heterozygous coding region mutations shared by small numbers of cases. The contribution of rare variation within non-coding regulatory regions to PAH aetiology remains to be determined. This will require functional annotation of regulatory and other noncoding regions specific for relevant cell types, further case-control analyses of these regions and ultimately functional studies of gene regulation to assess the pathogenicity of noncoding variants. Our findings to date provide support for a central role of the pulmonary vascular endothelium in disease pathogenesis, and suggest new mechanisms that could be exploited therapeutically in this life-limiting disease.
centre. An additional UK family diagnosed with HPAH was ascertained as described previously 45 . Blood and saliva samples were collected under written informed consent of the participants or their parents for use in gene identification studies. The non-PAH cohort for the case-control comparison comprised 6385 unrelated subjects recruited to the NIHR BR-RD study.
High-throughput sequencing DNA extracted from venous blood underwent whole genome sequencing using the Illumina TruSeq DNA PCR-Free Sample Preparation kit (Illumina Inc., San Diego, CA, USA) and Illumina HiSeq 2000 or HiSeq X sequencer, generating 100 -150 bp reads with a minimum coverage of 15X for ~95% of the genome (mean coverage of 35X). Whole-exome sequencing was conducted for individual II-1 ( Figure 4a) using genomic DNA extracted from peripheral blood. Paired-end sequence reads were generated on an Illumina HiSeq 2000.
(HumanCoreExome-12v1.1, HumanCoreExome-24v1.0, HumanOmni2.5-8v1.1), do not overlap quality control excluded regions or multiallelic sites in the 1000 Genomes (1000G) Phase 3 dataset 49 , do not have any missing genotypes in NIHR BR-RD, had a MAF of 0.3 or above and LD pruning was performed using PLINK 57 with a window size of 50 bp, window shift of 5 bp and a variance inflation factor threshold of 2. The 2,110 samples from the 1000G Project including the European (EUR), African (AFR), South Asian (SAS) and East Asian (EAS) populations (excluding the admixed American population) were filtered for the selected SNPs and the filtered data were used to perform a principal component analysis (PCA) using PC-Air. We modelled the scores of the leading five principal components as data generated by a population specific multivariate Gaussian distribution and estimated the corresponding mean and covariance parameters. Genotypes from the NIHR BR-RD samples were projected onto the loadings for the leading five principal components from the 1000G PCA and we computed the likelihood that each sample belonged to each subpopulation under a mixture of multivariate Gaussians models. Each sample was allocated to the population with the highest likelihood, unless the highest likelihood was similar to likelihood values for other populations, as might be expected for example under admixed ancestry or if the sample came from a population not included in 1000G. Such ambiguous samples were labeled as "other". PC-Relate was used to to identify related individuals in NIHR BR-RD. We used the first 20 PCs from PC-Air to adjust for relatedness and extracted the pairwise Identity-By-State distances and kinship values. The pairwise information was used by Primus to infer family networks and calculate the maximum set of unrelated samples.

Cohort definition and allele frequency calculation
Based on the relatedness analysis, we defined the following sample subsets: (a) the maximum number of unrelated non-PAH controls (UPAHC, n=6385), (b) all affected PAH cases (PAHAFF, n=1048), and (c) all unrelated PAH index cases (PAHIDX, n=1038). These subsets were used to annotate the variants in the multi-sample VCF file with calculated minor allele frequencies using the fill-tags extension of BCFtools 58 .

Rare variant filtering
Filtering of rare variants was performed as follows: 1) variants with a MAF less than 1 in 10,000 in UPAHC subjects, UK10K and ExAC were retained (adjusted for X chromosome variants to 1 in 8,000); 2) variants with a combined annotation dependent depletion deleteriousness (CADD) score of less than 15 were excluded. CADD scores were calculated using the CADD web service (http://cadd.gs.washington.edu) for variants lacking a score; 3) premature truncating variants (PTVs) or missense variants of the canonical transcript were retained; 4) missense variants predicted to be both tolerated and benign by SIFT and PolyPhen-2, respectively, were removed.
To identify likely causative mutations (as reported in Supplementary Table 5), variants in previously reported and putative genes, identified in this study, were examined in more detail to exclude variants that did not segregate in families (where data available). Furthermore, variants shared between cases and non-PAH controls, as well as variants of uncertain significance that co-occurred with previously reported causative mutations or high impact PTVs were also excluded.
Burden analysis of protein-truncating and missense variants Filtered variants were grouped per gene and consequence type (predicted PTV / missense) and subjects with at least one variant were counted (no double counting) per group and tested for association with disease. We applied a one-tailed Fisher's exact test with post hoc Bonferroni correction to calculate the P value for genome-wide significance.

Rare variant analysis using SKAT-O
To further investigate the aggregated effect that rare variants contribute to PAH aetiology, we applied a Sequence Kernel Association test (SKAT-O). SKAT-O increases the power of discovery under different inheritance models by combining variance-component and burden tests. Variants were filtered based on MAF as specified above, and only PTV and missense variants were included. For the analysis we implemented SKAT-O in RvTests v1.9.9 59 with default parameters and weights being Beta (1,25), and applying a correction for read length, gender and the first five principal components of the ethnicity PCA. Variants were collapsed considering only the protein-coding region in the canonical transcript of the protein-coding genes in the genome assembly GRCh37.

Analysis of large deletions
Copy number variation was identified using Canvas 60 and Manta 61 . Deletions called by both Manta and Canvas with a reciprocal overlap of ≥ 20% were retained. Of these, deletions were excluded if both failed standard Illumina quality metrics or overlapped with known benign deletions in healthy cohorts 62 . Deletions with a reciprocal overlap of ≥ 50% between samples were merged and filtered for a frequency of less than 1 in 1,000 in WGS10K and overlapping exonic regions of protein coding genes (GRCh37 genome assembly). The number of subjects with deletions were added up by gene (no double counting of subjects) and tested for association with the disease. We applied a one-tailed (greater) Fisher's exact test with Bonferroni post hoc correction for multiple testing to determine the P values for genome-wide significance.

Confirmation of variants
Variant sequencing reads for SNVs, indels and deletions were visualised for validation on Integrative Genomes Viewer (IGV) 18 , and were confirmed by diagnostic capture-based highthroughput sequencing, if the IGV inspection was not satisfactory. For the familial segregation analysis, linkage to the BMPR2 locus was first examined by microsatellite genotyping analysis. Mutation screening of the BMPR2, ACVRL1, ENG, AQP1 and SOX17 genes was conducted by capillary sequencing using BigDye Terminator v3.1 chemistry. All DNA fragments were resolved on an ABI Fragment Analyzer (Applied Biosystems). The family trees were drawn using the R package FamAgg 63 . 1:100 dilution). A section of the pulmonary arteriole was collected, fixed in formalin and embedded in paraffin, and sections analysed to ensure that the vessel was of pulmonary origin. Papworth Hospital ethical review committee approved the use of the human tissues (Ethics Ref 08/H0304/56+5) and informed consent was obtained from all subjects.
Human blood outgrowth endothelial cells (BOECs) were generated from 40-80 ml peripheral blood isolated with informed consent from healthy controls. The study was approved by the Cambridgeshire 3 Research Ethics Committee, reference number 11/EE/0297. BOECs were cultured in EGM-2MV with 10% standard grade FBS (Life Technologies, Carlsbad, CA) and were used for experiments between passages 4 and 8 65 . The endothelial nature of BOECs was assessed by flow cytometry for endothelial surface marker expression as described previously 31 . All cell lines were routinely tested for mycoplasma contamination. For experiments involving BOEC generation, Human pulmonary artery endothelial cells (PAECs) were purchased from Lonza (Cat. No. CC-2530; Basel, Switzerland) and maintained in EGM-2 with 2% fetal bovine serum (FBS; Lonza). PAECs were used between passages 4 and 8. For experimental studies, cells were treated with EBM-2 containing Antibiotic-Antimycotic (A/A; 100U ml -1 penicillin, 100 mg ml -1 streptomycin and 0.25 mg ml -1 amphotericin B, Invitrogen, Renfrewshire, UK) and FBS concentrations as stated. Cell lines were routinely tested for mycoplasma contamination and only used if negative.

RNA preparation and quantitative reverse transcription-PCR
Total RNA was extracted using the RNeasy Mini Kit with DNAse digestion (Qiagen, West Sussex, UK). cDNA was prepared from 1 µg of RNA using the High Capacity Reverse Transcriptase kit (Applied Biosystems, Foster City, CA), according to the manufacturer's instructions. All quantitative PCR reactions were prepared in MicroAmp optical 96-well reaction plates (Applied Biosystems) using 50 ng µl -1 cDNA with SYBR Green Jumpstart Taq Readymix (Sigma-Aldrich), ROX reference dye (Invitrogen) and custom sense and antisense primers (all 200 nmol l -1 ). Primers for human ACTB (encoding β-actin), AQP1, ATP13A3, B2M, HPRT and SOX17 were designed using PrimerBLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) (Supplementary Table 12). Reactions were amplified on a Quantstudio 6 Real-Time PCR system (Applied Biosystems). The relative abundance of each target gene in different cell lines was compared using the equation 2 -(CtGOI-Ct3HK) , where Ct3HK corresponded to the arithmetic mean of the Cts for ACTB, B2M and HPRT for each sample. For expression analysis of siRNA knockdown, the 2 -(ΔΔCt) method was used and fold expression determined relative to the DH1 control.

siRNA transfection
Prior to transfection, cells were preincubated in Opti-MEM-I reduced serum media (Invitrogen) for 2h before transfection with 10nM siRNA that had been lipoplexed for 20 min at RT with DharmaFECT1 (GE Dharmacon, Lafayette, CO). Cells were then incubated with the siRNA/DharmaFECT1 complexes for 4h at 37 o C before replaced by full growth media. Cells were kept in growth media for 24h before further treatment. Knockdown efficiency was confirmed by mRNA expression or immunoblotting. For proliferation assays, parallel RNA samples were collected both on day0 and day6, confirming that ATP13A3 expression was reduced by >90% on Day 0 and still reduced by >70% at Day 6. For all other assays, parallel RNA samples were collected on the day of the experiment to confirm knockdown, which was >90%. The siRNAs used were oligos targeting ATP13A3 (SASI_Hs02_00356805) from Sigma-Aldrich and ON-TARGETplus non-targeting Pool (siCP; GE Dharmacon).

Flow cytometric apoptosis assay
BOECs were plated 150,000/well into 6-well plates and transfected with siATP13A3 or siCP lipoplexed with DharmaFECT1. Cells were then serum-starved in EBM-2 (Lonza) containing 0.1% FBS and A/A for 8 hours before treating with EBM-2 and A/A containing either 0.1%FBS or 5%FBS for another 24 hours. Cells were then trypsinized and after washing with PBS, stained using the FITC Annexin V Apoptosis Dectection Kit I (BD Biosciences). For each condition, dual-staining of 5µl FITC conjugated Annexin V and 5µl propidium iodide (PI) were added and incubated at room temperature for 15 minutes. For the single staining controls for compensation, either 5µl FITC Annexin V or 5µl PI was added into nontransfected cells. All samples were analysed on BD Accuri™ C6 Plus platform (BD Biosciences). Data were collected and analysed using FlowJo software, with AnnexinV + /PIcells defined as early apoptotic (Treestar).

Caspase-Glo 3/7 assay
BOECs were seeded at a density of 150,000/well into 6-well plates and transfected with siATP13A3 or siCP lipoplexed with DharmaFECT1. For each condition, cells were trypsinized from 6-well plates and reseeded in triplicates into a 96-well plate at a density of 15,000-20,000/well and left to adhere overnight. Cells were quiesced in EBM-2 containing 0.1%FBS for 24h before treating with or without EBM-2 and A/A containing either 0.1%FBS or 5%FBS for 16 hours. For measuring caspase activities, 100ul Caspase-Glo® 3/7 Reagent (G8091 Promega) was added into each well, incubated and mixed on a plate shaker in the dark for 30 minutes at room temperature. The lysates were transferred to a white-walled 96well plate and luminescence was read in a GloMax® luminometer (Promega).

Accession codes
Data of PAH cases included in this manuscript and eligible for public release according to the UK Research Ethics rules have been deposited in the European Genome-phenome Archive (EGA) at the EMBL -European Bioinformatics Institute under accession number EGAD00001003423.            processing. The pre-pro-protein is processed into the mature growth factor domain (GFD) bound to the prodomain upon secretion 23 . (b) Plot of GDF2 mutations found only in PAH cases superimposed on the structure of prodomain bound GDF2. (pdb code 4YCG) 22 . The GDF2 growth factor domain is shown in green and the prodomain in cyan. (c) Magnified view of the Arg110 and Glu143 mutations. The wild type amino acids make double salt bridges to stabilise the prodomain conformation at the interface between the growth factor domain and prodomain. The E143K and R110W mutations both disrupt these interactions, destabilising the interaction between the growth factor domain and prodomain. (d) GDF2 levels secreted into supernatants of HEK293T cells transfected with likely pathogenic variants found in PAH cases, compared with wild type GDF2 and cells transfected with an empty vector. *** P < 0.001 by ANOVA.
. CC-BY-NC-ND 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted January 26, 2018. ; https://doi.org/10.1101/185272 doi: bioRxiv preprint . CC-BY-NC-ND 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted January 26, 2018. ; https://doi.org/10.1101/185272 doi: bioRxiv preprint Figure 6 Structural analysis of ATP13A3 mutations. (a) Topology of ATP13A3, plotted according to UniProtKB Q9H7F0. Frameshift and stop-gained mutations identified in PAH cases are shown as khaki circles, and missense mutations as red circles. Frameshift/stopgained mutations are predicted to truncate the protein prior to the catalytic domain and essential Mg binding sites, leading to loss of ATPase activity. (b) Sequence alignment of ATP13A3 with ATP1A1 (P05024), of which the high resolution structure was used for the structural analysis in (c). The conserved regions of ATP13A3 and ATP1A1, essential for ATPase activity 24 , show good alignment (data not shown). Only regions containing the missense PAH mutations are shown, with positions of the four missense mutations highlighted in yellow above the sequences. (c) Structural analysis of the 4 PAH missense mutations plotted on the ATP1A1 crystal structure based on the sequence alignment in (b) (pdb code 3wgu) 25 . Green: α subunit (P05024), cyan: β subunit (P05027), grey: γ-subunit transcript variant a (Q58k79). Y535, Y677, R685 and I787 are the numbering in ATP1A1. Positions of the four missense mutations found in PAH are labeled and highlighted by red circles. (d) Magnified view of the cytoplasmic region of the ATPase, showing the presence of ADP at the active site. The conserved regions essential for ATPase activity are shown in light pink. The L675V and R858H mutations are located close to the ATP catalytic region.
. CC-BY-NC-ND 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted January 26, 2018.  with seven other mammals. The bovine AQP1 has the high resolution (2.2Å) published structure. Mutations identified in PAH cases are highly conserved and highlighted in yellow. (b) Crystal structure of bovine AQP1 (pdb code 1j4n) 26 . Left: side view; right: top view from the extracellular direction. AQP1 is shown as a semi-transparent cartoon and five water molecules in the water channel are shown as red spheres. Key residues lining the water channels are represented with stick structures. (c) Magnified view of the water channel, with H-bonds connected to water molecules in the channel highlighted. Two asparagine-proline-alanine (NPA) motifs, essential for the water transporting function of AQP1, are shown in magenta. Conserved His180 that constricts the water channel is shown in yellow. Mutations found in PAH cases, Arg195Trp and Val176Glu, are labelled and shown as orange stick structures. Arg195 and His180 are highly conserved in the known water channels and are strong indicators of water channel specificity. Arg195Trp and Val176Glu mutations are predicted to disrupt the conformation of this conserved water channel.