Rare variant analysis of 4241 pulmonary arterial hypertension cases from an international consortium implicates FBLN2, PDGFD, and rare de novo variants in PAH

Pulmonary arterial hypertension (PAH) is a lethal vasculopathy characterized by pathogenic remodeling of pulmonary arterioles leading to increased pulmonary pressures, right ventricular hypertrophy, and heart failure. PAH can be associated with other diseases (APAH: connective tissue diseases, congenital heart disease, and others) but often the etiology is idiopathic (IPAH). Mutations in bone morphogenetic protein receptor 2 (BMPR2) are the cause of most heritable cases but the vast majority of other cases are genetically undefined. To identify new risk genes, we utilized an international consortium of 4241 PAH cases with exome or genome sequencing data from the National Biological Sample and Data Repository for PAH, Columbia University Irving Medical Center, and the UK NIHR BioResource – Rare Diseases Study. The strength of this combined cohort is a doubling of the number of IPAH cases compared to either national cohort alone. We identified protein-coding variants and performed rare variant association analyses in unrelated participants of European ancestry, including 1647 IPAH cases and 18,819 controls. We also analyzed de novo variants in 124 pediatric trios enriched for IPAH and APAH-CHD. Seven genes with rare deleterious variants were associated with IPAH with false discovery rate smaller than 0.1: three known genes (BMPR2, GDF2, and TBX4), two recently identified candidate genes (SOX17, KDR), and two new candidate genes (fibulin 2, FBLN2; platelet-derived growth factor D, PDGFD). The new genes were identified based solely on rare deleterious missense variants, a variant type that could not be adequately assessed in either cohort alone. The candidate genes exhibit expression patterns in lung and heart similar to that of known PAH risk genes, and most variants occur in conserved protein domains. For pediatric PAH, predicted deleterious de novo variants exhibited a significant burden compared to the background mutation rate (2.45×, p = 2.5e−5). At least eight novel pediatric candidate genes carrying de novo variants have plausible roles in lung/heart development. Rare variant analysis of a large international consortium identified two new candidate genes—FBLN2 and PDGFD. The new genes have known functions in vasculogenesis and remodeling. Trio analysis predicted that ~ 15% of pediatric IPAH may be explained by de novo variants.


Background
Pulmonary arterial hypertension (PAH) remains a progressive, lethal vasculopathy despite recent therapeutic advances. The disease is characterized by pulmonary vascular endothelial dysfunction and proliferative remodeling giving rise to increased pulmonary artery pressures and pulmonary vascular resistance. These pathological changes of the lung vasculature strain the right ventricle of the heart, leading to right ventricular hypertrophy, right heart failure, and high mortality [1][2][3]. Dysregulated vascular, inflammatory, and immune cells contribute to these pathological processes [3]. PAH can present at any age, but the~3:1 female to male ratio in adult-onset disease is not observed in pediatric-onset disease, in which the disease incidence is similar for males and females. The estimated prevalence of PAH is 4.8-8.1 cases/million for pediatric-onset [4] and 5.6-25 cases/million for adult-onset disease [5]. Early genetic linkage and candidate gene studies indicated an autosomal dominant mode of inheritance for PAH risk. However, the known susceptibility variants are incompletely penetrant, many individuals who carry monogenic risk variants never develop PAH, and a subset of patients have deleterious variants in more than one risk gene. For example, bone morphogenetic protein receptor type 2 (BMPR2) mutations are observed in 60-80% of familial (FPAH) cases, but data from population registries indicate that penetrance of the disease phenotype ranges from 14 to 42% [6]. These data suggest that additional genetic, epigenetic, environmental factors, and gene × environment interactions contribute to disease.
Together, these data indicate that rare genetic variants underlie~75-80% of FPAH [6], at least 10% of adultonset idiopathic PAH (IPAH) [7,8], and up to~36% of pediatric-onset IPAH [31]. A substantial fraction of nonfamilial PAH cases remains genetically undefined. The low frequency of risk variants for each gene, except BMPR2, indicates that large numbers of individuals are required for further validation of rare risk genes and pathways, and to understand the natural history of each genetic subtype of PAH. Towards this end, we analyzed 4175 PAH cases from an international consortium with ES or GS. The National Biological Sample and Data Repository for PAH (aka PAH Biobank) was comprised of 2570 PAH cases (1110 IPAH and 1239 APAH) and the UK NIHR BioResource -Rare Diseases Study was comprised of 1144 cases, almost entirely IPAH. Thus, the increased power of the combined cohort was a 2-fold increase in the number of IPAH cases, and we focused our association analyses on this PAH subclass. The cohort size precluded testing of the oligogenicity hypothesis suggested by the incomplete penetrance of known PAH risk genes. Non-inherited de novo mutations could also contribute to genetically unexplained non-familial cases but require access to parental sequencing data. We previously showed that pediatric-onset PAH cases were enriched with damaging de novo variants. Here, we expand the analysis to a cohort of 124 pediatric childparent trios.

Patient cohorts and control datasets
A total of 4175 PAH cases from the National Biological Sample and Data Repository for PAH (PAH Biobank, n = 2570 exomes) [7], UK NIHR BioResource -Rare Diseases Study (UK NIHR BioResource, n = 1144 genomes) [8], and the Columbia University Irving Medical Center (CUIMC, n = 461 exomes) [17,18,27] were included in a combined analysis of rare inherited variants. The subset of 124 affected child-unaffected parents trios (n = 111 CUIMC, n = 8 UK NIHR BioResource, n = 5 PAH Biobank) were included in an analysis of de novo variants. An additional 65 BMPR2 mutation-positive cases from CUIMC without exome sequencing data were previously reported [17,18] and included in the overall cohort counts (total of 4241 cases). As previously described, cases were diagnosed by medical record review including right heart catheterization and all were classified as World Symposium on Pulmonary Hypertension (WSPH) Group I [32]. Written informed consent for publication was obtained at enrollment. The studies were approved by the institutional review boards at CCHMC, individual PAH Biobank Centers, the East of England Cambridge South national research ethics committee (REC, ref. 13/EE0325) or CUIMC.

ES/GS data analysis
PAH Biobank, CUIMC, and SPARK cohort samples were all sequenced in collaboration with the Regeneron Genetics Center as previously described [7,8,17,18,27]; the UK NIHR BioResource sequence data were also previously described [8]. For case and SPARK control data, we used a previously established bioinformatics procedure [34] to process and analyze exome and genome sequence data. For the UK NIHR BioResource data, we extracted reads from GS data by the following procedure: (1) obtained all reads that were mapped to the human genome regions that overlapped with the target regions of xGEN exome capture intervals (Exome Research panel 1.0); (2) the mate pairs of these reads. We then processed the extracted GS data using the same pipeline as the ES data. Specifically, we used BWA-MEM [35] to map and align paired-end reads to the human reference genome (version GRCh38/hg38, accession GCA 000001405.15), Picard v1. 93 MarkDuplicates to identify and flag PCR duplicates, and GATK v4.1 [36,37] HaplotypeCaller in Reference Confidence Model mode to generate individual-level gVCF files from the aligned sequence data. We then performed joint calling of variants from all three datasets using GLnexus [38]. We used the following inclusion rules to select variants for downstream analysis: AF < 0.05% in the cohort, < 0.01% in gnomAD exome_ALL (all ancestries); > 90% target region with dp ≥ 10; mappability = 1; and allele balance ≥ 0.25. We also ran DeepVariants [39,40], a new tool based on machine learning, for all cases and SPARK controls. We used the ES mode for ES data and GS mode for GS data, and then filtered by "PASS" DeepVariants. Inclusion criteria for variants observed in multiple carriers was ≥ 50% of all calls PASS DeepVariants. For gnomAD data, only variants located in xGen-captured proteincoding regions were used; filtering was based on GATK metrics obtained from gnomAD and only "PASS" variants were included. SNVs with VQSR <− 20 and indels with VQSR <− 5 were excluded. Variants used for downstream analyses were restricted to the subset called by both GLnexus and DeepVariants.
De novo variants were defined as a variant present in the offspring with homozygous reference genotypes in both parents. We used a series of filters to identify de novo variants: VQSR tranche ≤ 99.7 for SNVs and ≤ 99.0 for indels; GATK Fisher Strand ≤ 25; quality by depth ≥ 2. We required the candidate de novo variants in probands to have ≥ 5 reads supporting the alternative allele, ≥ 20% alternative allele fraction, Phred-scaled genotype likelihood ≥ 60 (GQ), and population AF ≤ 0.01% in ExAC and required both parents to have ≥ 10 reference reads, < 5% alternative allele fraction, and GQ ≥ 30.
We used Ensembl Variant Effect Predictor (VEP; Ensemble 93) [41] to annotate variant function and ANNOVAR [42] to aggregate variant population frequencies and in silico predictions of deleteriousness. Rare synonymous variants were further evaluated with SpliceAI [43] to identify cryptic splice site variants (score ≥ 0.5). Rare variants were defined as AF ≤ 0.01% in gnomAD exome_ALL (all ancestries). A total of 18,939 protein-coding genes were identified containing ≥ 1 rare variant, excluding mucin and major histocompatibility complex genes due to low sequence complexity. Deleterious variants were defined as likely gene-disrupting (LGD, including premature stop-gain, frameshift indels, canonical splicing variants, cryptic splice site variants, and exon deletions) or predicted damaging missense (D-Mis) based on gene-specific REVEL score thresholds [18,44] (see below). All rare inherited and de novo variants in candidate genes were manually inspected using Integrative Genome Viewer (IGV) [45]. Indels were confirmed independently by Sanger sequencing.

Statistical analysis
To identify novel risk genes for IPAH, we performed a rare variant association test in unrelated participants of European ancestry. Genetic ancestry and relatedness of cases and SPARK controls were checked using Peddy [46], and only unrelated cases (n = 2789) and controls (18,819: 11,101 SPARK parents and 7718 gnomAD individuals) were included in the association test. The gno-mAD controls were confined to non-Finnish Europeans (NFE). We performed a gene-based case-control test comparing the frequency of rare deleterious variants in PAH cases with unaffected controls. To reduce batch effects in combined datasets from different sources [47], we limited the analysis to regions targeted by xGen and with at least 10× coverage in 90% of samples. We then tested for similarity of the rare synonymous variant rate among cases and controls, assuming that most rare synonymous variants do not have discernible effects on disease risk.
To identify PAH risk genes, we tested the burden of rare deleterious variants (AF ≤ 0.01%, LGD or D-Mis) in each protein-coding gene in cases compared to controls using a variable threshold test [48]. Specifically, we used REVEL [44] scores to predict the deleteriousness of missense variants, searched for a gene-specific optimal REVEL score threshold that maximized the burden of rare deleterious variants in cases compared to controls, and then used permutations to calculate statistical significance as described previously [7] to control the type I error rate. We checked for inflation using a quantilequantile (Q-Q) plot and calculated the genomic control factor, lambda, using QQperm (https://cran.r-project. org/web/packages/QQperm/QQperm.pdf).
Lambda equal to 1 indicates no deviation from the expected distribution. We performed two association tests, one with LGD and D-Mis variants combined and the other with D-Mis variants alone. We defined the threshold for genome-wide significance by Bonferroni correction for multiple testing (n = 40,000, 18,939 protein-coding genes containing rare variants times two tests for each gene, yielding a threshold p value = 1.25e−6). We used the Benjamini-Hochberg procedure to estimate false discovery rate (FDR) by p.adjust in R.
To test whether recurrent variants in individual genes represented independent mutational events or were due to founder events, we first tested for relatedness among samples using KING [49], in addition to Peddy [46]. None of the cases with recurrent variants had any evidence of relatedness. Second, we assessed shared haplotypes of recurrent variant carriers using SHAPEIT2 [50] and the HapMap genetic map [51]. Since all of the recurrent variant carriers were of European ancestry, we restricted the HapMap data to the European population.
To estimate the burden of de novo variants in cases, we calculated the background mutation rate using a previously published tri-nucleotide change table [52,53] and calculated the rate in protein-coding regions that are uniquely mappable. We assumed that the number of de novo variants of various types (e.g., synonymous, missense, LGD) expected by chance in gene sets or all genes followed a Poisson distribution [52]. For a given type of de novo variant in a gene set, we set the observed number of cases to m1, the expected number to m0, estimated the enrichment rate by (m1/m0), and tested for significance using an exact Poisson test (poisson.test in R) with m0 as the expectation.

Protein modeling
Homology structures of conserved protein domains in FBLN2 and PDGFD were built using EasyModeller 4.0 [54]. Template structures were downloaded from the protein database (PDB) for endothelial growth factor (EGF, PDB ID 5UK5) and CUB (PDB ID 3KQ4) domains. The template structure for platelet-derived growth factor (PDGF)/vascular EGF (VEGF) was downloaded directly from PrePPI [55,56].

Gene expression
Single-cell RNA-seq data of aorta, lung, and heart tissues were obtained from Tabula Muris, a transcriptome compendium containing RNA-seq data from~100,000 single cells from 20 adult-staged mouse organs [57]. We chose 14 tissue/cell types including endothelial, cardiac muscle, and stromal cells from the three tissues, restricting the analysis to tissues/cell types for which there was RNA-seq data from at least 70 individual cells (Additional file 1, Supplementary Figure 1). Relative gene expression was based on the fraction of cells with > 0 reads in each cell type. PCA of cell type-specific gene expression profiles was performed using a script available through GitHub [58].
A comparison of the clinical characteristics and hemodynamic data for pediatric-versus adult-onset PAH cases is shown in Additional file 2 (Supplementary Table 1). Notably, the female:male ratio among pediatric-onset cases was significantly lower (1.65:1) compared to adult-onset cases (4:1, p < 0.0001 by Fisher's exact test), and children had higher mPAP and mPCWP, decreased cardiac output and increased pulmonary vascular resistance compared to adults at diagnosis (all differences p < 0.0001 by Student's t test).

Identification of novel risk genes: FBLN2 and PDGFD
To perform a combined analysis of US and UK sequencing data, we reprocessed the UK data using our inhouse pipeline, including predictions of missense variant deleteriousness [7]. Quality control procedures included detection of cryptic relatedness among all PAH participants. We performed a gene-based case-control association analysis to identify novel PAH risk genes using only unrelated cases. To control for population stratification, we confined the association analysis to individuals of European ancestry (2789 cases, 18,819 controls) and then screened the whole cohort, including nonEuropeans, for rare deleterious variants in associated genes. As a quality control check for the filtering parameters employed, we compared the frequencies of rare synonymous variants, a variant class that is mostly neutral with respect to disease status, in European cases vs controls. We observed similar frequencies of synonymous variants in cases vs controls (enrichment rate = 1.0, p value = 0.28) (Additional file 2, Supplementary Table  2). Furthermore, a gene-level burden test revealed no enrichment of rare synonymous variants in cases (Additional file 1, Supplementary Figure 2). We then proceeded to test for gene-specific enrichment of rare deleterious variants (AF < 0.01%, LGD and D-Mis, or D-Mis only) in cases compared to controls. We note that to improve power, we empirically determined the optimal REVEL score threshold to define deleterious missense variants in a gene-specific manner using a variable threshold test [7]. To account for potential different modes of action for different risk genes, we tested the association twice for each gene: one with LGD and D-Mis variants and the other with D-Mis variants alone. In this approach, LGD and D-Mis together is optimized for complete or partial loss of function; D-Mis alone is optimized for gain of function or dominant negative variants. We set the total number of tests at twice the number of protein-coding genes for multiple test adjustment, a conservative approach considering that the data used in these two tests per gene are not independent. The Q-Q plot of p values from tests in all genes shows negligible genomic inflation (Additional file 1, Supplementary Figure 3). Rare deleterious variants in eleven genes were significantly associated (false discovery rate, FDR < 0.1) with PAH. Among these, seven are known or previously reported candidate PAH risk genes: BMPR2, TBX4, GDF2, ACVRL1, SOX17, AQP1, ATP13A3, and KDR. Three are new candidate genes: COL6A5 (collagen type VI alpha 5 chain), JPT2 (Jupiter microtubuleassociated homolog 2), and FBLN2 (fibulin 2).
The increased power inherent to the combined cohort over the PAH Biobank or UK NIHR BioResource alone is due to a twofold increase in the number of IPAH cases, including the number of European cases used for association analysis. Power analyses indicated that the study had ample power to detect risk genes with large effect size and modest variant allele frequency, or large variant allele frequency and modest effect size, relative to IPAH risk genes identified in smaller cohorts (Additional file 1, Supplementary Figure 4). To take advantage of the increased number of European IPAH cases in the combined cohort, we then restricted the analysis to IPAH. Again, testing for association across all proteincoding genes for 1647 IPAH cases compared to 18,819 controls was generally consistent with expectation under the null model (Fig. 1). Rare predicted deleterious variants in seven genes were significantly associated (FDR < 0.1) with IPAH, including three known genes (BMPR2, GDF2, and TBX4), two recently identified candidate genes (SOX17 and KDR), and two new candidate genes (FBLN2, and PDGFD, platelet-derived growth factor D). More than 95% of samples for both cases and controls had at least 10× depth of sequence coverage across the target regions for FBLN2 and PDGFD (Additional file 1, Supplementary Figure 5), excluding the possibility that the associations were driven by coverage differences between cases and controls. We also tested for gene-level associations restricting the analysis to European APAH cases (n = 998). The Q-Q plot of p values from all gene tests is shown in Additional file 1, Supplementary Figure  6. Known PAH gene ACVRL1 showed association with APAH, consistent with its role in APAH-HHT, but no genes were significantly associated at FDR < 0.1.
KDR has recently been implicated as a causal gene for PAH based on a small familial study [60] and our population-based phenotype-driven (SKAT-O) analysis of the UK NIHR BioResource cohort with replication in the PAH Biobank [61]. Both of those analyses were based on protein-truncating variants. Herein, we provide additional statistical evidence based on a burden test including both LGD and D-Mis variants using our variable threshold method. Six cases (5 IPAH, 1APAH-CHD) carry D-Mis variants with empirically determined REVE L > 0.86; details of the variants are provided in Supplementary Table 3. All of the variants are located in the conserved tyrosine kinase domain of the encoded protein (www.uniprot.org). One of the variants, c.3439C>T is recurrent in three cases. There was no evidence of relatedness for these cases, and the relatively short shared haplotype length and common population frequency (Additional file 2, Supplementary Table 4) indicate that the variant occurrences represent independent mutational events rather than being derived from a founder event. None of these cases have variants in other known PAH risk genes. The age-of-onset for the six cases is 57 ± 20 years (mean ± SD, range 25-75 years) and all are of European ancestry. Statistically significant association following Bonferroni correction for multiple testing provides confirmation of the association of KDR with PAH using an alternative burden-based statistical method.
The associations of FBLN2 and PGDFD were both driven by D-Mis variants. We next screened the entire combined cohort, including participants of non-European ancestry, for rare deleterious missense variants in FBLN2 and PDGFD. In total, seven cases carry FBLN2 variants (6 IPAH, 1 APAH) and ten cases carry PDGFD variants (9 IPAH, 1 PAH associated with diet and toxins) ( Table 2). Most of the carriers are of European ancestry; one FBLN2 carrier is of East Asian ancestry and one PDGFD carrier is of African ancestry. One    residues in an EGF domain which may change the conformation of the domain and impact protein function (Fig. 2b). Recurrent FBLN2 p.(Asp982Tyr) disrupts a Ca ++ binding site [62] in another EGF domain (Fig. 2b), which may reduce the affinity and frequency of Ca ++ binding. PDGFD p.(Asp148Asn) disrupts a Ca ++ binding site within the CUB domain [63] (Fig. 2c) and recurrent PDFGD p. (Tyr321Asn) is predicted to disrupt a hydrogen bond within the PDGF/VEGF domain (Fig. 2c). In addition, PDGFD p.(Arg295Cys) is located in close proximity to Cys356 and Cys358, potentially introducing new disulfide bonds within the PDGF/VEGF domain.

Clinical phenotypes of FBLN2 and PDGFD variant carriers
The clinical phenotypes of all FBLN2 and PDGFD variant carriers are provided in

Gene expression patterns of PAH candidate risk genes
We hypothesized that PAH risk genes are highly expressed in certain cell types relevant to the disease etiology and that joint analysis of cell type-specific expression data with genetic data could inform cell types associated with disease risk [66]. We obtained single-cell RNA-seq data of aorta, lung, and heart tissues available through the Tabula Muris project, a transcriptome compendium containing RNA-seq data from adultstaged mouse organs [57]. We chose 14 tissue/cell types including endothelial, cardiac muscle, and stromal cells as a proxy for the cell types of the pulmonary artery (unavailable). A list of the tissues, cell types, and the number of cells sequenced per tissue/cell type is provided in Additional file 1 (Supplementary Figure 1a). We queried gene expression for twelve known PAH risk genes (ACVRL1, BMPR2, CAV1, EIF2AK4, ENG, KCNK3, KDR, NOTCH1, SMAD4, SMAD9, SOX17, TBX4) and the two new candidate risk genes (FBLN2, PDGFD). A heat map with hierarchical clustering of relative gene expression is shown in Fig. 3a. The majority of known risk genes (7/12) have relatively high expression in endothelial cells from the three tissues; most others have high expression in tissue-specific cardiac muscle, stromal cells, or fibroblasts. PDGFD is located in the same cluster as BMPR2, SOX17, and KDR; these genes are specifically and highly expressed in endothelial cell types.
FBLN2 is highly expressed in both endothelial and fibroblast cell types. We then randomly selected a set of 100 genes without reported associations with PAH and performed PCA of cell type-specific expression profiles of known risk genes and random genes. The second component (PC2) largely separates known risk genes and random genes (Fig. 3b, c). Consistent with hierarchical clustering, endothelial expression in all three tissues was positively correlated with PC2 (Additional file 1, Supplementary Figure 1b). Projecting all protein-coding genes onto PC2, seven of twelve known risk genes are ranked in top 5% among all genes (Fig. 3d) (binomial test: enrichment =20, p = 1.6E−05). Two new candidate genes, FBLN2 and PDGFD, are ranked in the top 1.8% of PC2.

Identification of novel candidate pediatric PAH risk genes by de novo variant analysis
We next focused on pediatric-onset disease, a subpopulation in which genetic factors likely play a larger causal role compared to adults. The study was underpowered to carry out a gene-based case-control association analysis due to the relatively small number of pediatric patients (n = 442); however, 124 pediatriconset PAH probands with child-parent trio data were available for de novo variant analysis. The trio cohort consisted mostly of IPAH (55.6%, n = 66) and APAH-CHD (37.9%, n = 45) cases. We performed a burden test for enrichment of exonic de novo variants among all trio probands by comparing the number of variants observed vs expected based on the background mutation rate. Similar rates of de novo mutations were observed for synonymous, LGD alone, and total missense variants (Table 4). However, there was a significant burden of D-Mis and LGD + D-Mis variants among cases over that expected by chance ( LGD + D-Mis variants in genes highly expressed in developing heart and lung, still significantly more than expected (enrichment rate = 1.86, p = 0.008, 10 expected risk variants). We tested the burden of de novo variants among IPAH cases and observed enrichment of D-Mis and LGD + D-Mis variants similar to that of the overall trio cohort (Additional file 2,  Table 6). Similar to other early-onset severe diseases, including CHD and bronchopulmonary dysplasia, the genes identified fit a general pattern for developmental disorders-genes intolerant to loss of function variants (pLI > 0.5 for 40% of the genes) and with known functions as transcription factors, RNAbinding proteins, protein kinases, and chromatin modification. Three of the genes are known CHD risk genes (NOTCH1, PTPN11, and RAF1), and 37% of the genes are known causal genes for a variety of developmental syndromes. Case variant PTPN11 p.(Asp61Gly) is a known causal variant for Noonan syndrome [68], and RAF1 p.Pro261 is a hotspot for multiple gain-of-function mutations, including p.(Pro261Thr), causing Noonan syndrome [69].

Clinical phenotypes of pediatric de novo variant carriers
Among the 36 patients who carry LGD or D-Mis de novo variants (Additional file 2, Supplementary Table 7), there is a 1.8:1 ratio of females to males, a mean age-ofonset of 5.4 ± 4.6 years, 50% of the cases (n = 18) have a diagnosis of IPAH, 33.3% (n = 12) APAH-CHD and an overlapping but distinct 36.1% of cases have other congenital or growth and development anomalies. NOTCH1 variant carrier, JM1357, has a diagnosis of APAH-CHD with tetralogy of Fallot, and a recent exome sequencing study of~800 tetralogy of Fallot cases identified NOTCH1 as the top association signal [70]. PTPN11 variant carrier, JM155, has a diagnosis of APAH-CHD associated with Noonan syndrome and the c.182A>G variant is known to be pathogenic in Noonan syndrome. Variants in PSMD12 cause Stankiewicz-Isidor syndrome, sometimes associated with congenital heart defects, and variant carrier 06-095 has a diagnosis of APAH-CHD. Hemodynamic data for the de novo variant carriers (Additional file 2, Supplementary Table 7) was similar to that of all pediatric cases in the cohort (Additional file 2, Supplementary Table 1).

Discussion
Combined analysis of a large US/UK cohort enriched in adult-onset IPAH cases enabled identification of five known and two new IPAH candidate risk genes with FDR < 0.1: FBLN2 and PDGFD are the new genes. The association was based on a gene-level case-control analysis of 1647 unrelated European IPAH cases. The variants contributing to FBLN2 and PDGFD associations are D-Mis variants predicted to alter highly conserved protein conformation, Ca ++ binding sites, or intramolecular binding sites within conserved protein domains, likely leading to important structural changes in critical domains. The non-founder nature of recurrent FBLN2 p.(Asp982Tyr) (n = 4 cases), and two PDGFD variants recurrent in two unrelated cases each, adds further support for pathogenicity of these alleles. In addition, we confirmed the recent association of KDR with PAH [60,61] based on an alternative statistical approach. We further show that all three of these candidate genes have high expression in lung and heart endothelial cell types, similar to other well established risk genes (BMPR2 and SOX17), further supporting the plausibility of these genes contributing to PAH risk. De novo variant analysis of pediatric-onset PAH (124 trios) showed a 2.45× enrichment of rare deleterious exonic variants, indicating that de novo variants contribute to~15% of pediatric cases across PAH subtypes. The de novo variants implicate new candidate risk genes likely unique to pediatric PAH, but some of the molecular pathways may inform both pediatric-and adult-onset PAH. FBLN2 encodes an extracellular matrix protein important for elastic fiber formation and regulation of cell motility, proliferation, and angiogenesis. FBLN2 is expressed in the lung vasculature but most studies have focused on gene expression in the heart vasculature. In mice, Fbln2 is expressed in epithelial-mesenchymal transformation during embryonic heart development and is upregulated postnatally throughout coronary vasculogenesis and angiogenesis when transformed mesenchymal cells migrate to the extracellular matrix [71,72]. Fbln2 -/mice are viable, fertile, and have intact elastic fiber formation, attributable to compensation by the more widely expressed Fbln1 gene [73,74]. However, Fbln2 expression is required for angiotensin II-induced TGFβ signaling and cardiac fibrosis [75]. In humans, FBLN2 variants have been associated with atrioventricular septal defects [76] and intracranial aneurysm [77], providing additional support for a role in vascular remodeling. We hypothesize that, in the pulmonary vasculature, gain of function variants may lead to increased TGF-β signaling, increased proliferation and medial hypertrophy. The FBLN2 protein contains 10 consecutive EGF protein-protein interaction domains, nine of which are calcium-binding. All seven of the case variants are missense variants, two of which are predicted to alter the conformation of an EGF domain, and a recurrent variant carried by four cases is predicted to disrupt the Ca ++ binding site of another EGF domain. The carriers of FBLN2 variants have adult-onset disease with mean age-of-onset similar to the overall cohort or IPAH alone. Five of seven carriers also have a diagnosis of systemic hypertension (HTN), and it is possible that gene damaging variants in FBLN2 contribute to the development of HTN. However, given the frequency of HTN in the overall US/UK combined cohort (32% for adult-onset IPAH; similar to that reported in the REVEAL registry [65]), there may be other age-related genetic and environmental factors contributing to HTN. Finally, two of our cohort cases, 08-018 and 29-031, have additional diagnoses of renal or heart anomalies, and FLBN2 has been identified as a key regulator of development in those tissues [78][79][80].
PDGFD is a member of the PDGF family that functions in recruiting cells of mesenchymal origin during development or to sites of injury [81]. PDGFD is widely expressed including arterial endothelial cells, adventitial pericytes and smooth muscle cells, lung endothelial cells, and smooth muscle cell progenitors of distal pulmonary arterioles. Secreted PDGFD specifically binds PDGFRβ, a widely expressed protein that co-localizes with PDGFD in vascular smooth muscle cells. Pdgfd knockout mice are phenotypically normal with the exception of a modest increase in systemic blood pressure [82], However, cardiac-specific PDGFD transgenic mice, overexpressing the active core domain of human PDGFD in the heart, exhibit vascular smooth muscle cell proliferation, vascular remodeling with wall thickening, severe cardiac fibrosis, heart failure, and premature death [83]. While effects of Pdgfd overexpression on the pulmonary vasculature have not been investigated, the cardiac vasculature data are consistent with a gain of function mechanism. Further evidence for the role of PDGFD as an effector molecule in cardiovascular diseases and cancer has been reviewed [81,84,85]. Full-length PDGFD contains two conserved protein domains, an autoinhibitory CUB domain and an enzymatic PDGF/VEGF domain; the protein undergoes proteolytic cleavage at Arg247 or Arg249 to produce an active growth factor promoting angiogenesis and vascular muscularization [86]. All ten of the case variants are missense variants; four reside in the CUB domain and five reside in the active processed protein. Variant p.(Asp148Asn), carried by two patients, is predicted to disrupt the Ca ++ binding site of the CUB domain; variants p.(Arg295Cys) and p.(Ser309Cys), carried by one and two patients respectively, are predicted to alter the conformation of the PDGF/VEGF domain. All but one of the PDGFD variant carriers have adultonset disease with mean age-of-onset similar to the overall cohort or IPAH alone. Four out of ten of the PDGFD variant carriers have additional diagnoses of other pulmonary fibrotic and/or vascular fibrotic diseases including bronchopulmonary dysplasia, emphysema, asthma, and one patient (E010173) with both mixed pulmonary valve disease and peripheral vascular disease (Table 3). Targeting the PDGF pathway with small molecule inhibitors of tyrosine kinase is an active area of investigation and several inhibitors are FDAapproved [87]. Notably, imatinib reduced cardiac fibroblast proliferation and PDGFD expression 15-fold [88]; data regarding effects on pulmonary arterial smooth muscle cells are warranted. A limitation of tyrosine kinase inhibitors is that they target multiple tyrosine kinases. Sequestering PDGFD with neutralizing antibodies or DNA/RNA aptamers, or preventing PDGFD-PDGFRβ interaction via oligonucleotides, may provide more specific targeting.
To test the plausibility of the new candidate PAH genes identified by association analysis, we leveraged publicly available single-cell RNA-seq data. PDGFD, and recently identified KDR, have very similar expression patterns as BMPR2 and SOX17, two established PAH genes. PCA indicated that the PAH risk genes can largely be separated from non-risk genes based on PC2. The majority of known PAH risk genes rank in the top 5% of PC2 among 16,744 genes queried, and the new genes-FBLN2 and PDGFD-rank within the top 1.8%, providing support for their candidacy as PAH risk genes. Other risk genes, like KCNK3 and EIF2AK4, exert important PAH-related functions in cell types other than endothelial cells, and GDF2 is excreted from liver; thus, it will be important to consider expression patterns on a gene-specific basis. In addition, the dataset utilized in this study was based on adult-staged murine cells and is not well-suited for developmental genes such as TBX4 and other genes likely to contribute to pediatric-onset disease. Thus, additional datasets from different time points are needed.
Rare deleterious variants in established PAH genes are clearly pathogenic based on segregation analyses, enrichment of rare deleterious variants in PAH cases compared to controls with replication over time, and demonstrated loss of function or aberrant function in vitro, in vivo (model organisms), or ex vivo [6]. However, none of the PAH genes are fully dominant and many carriers are never diagnosed with PAH. BMPR2 variants exhibit variable penetrance with~14% penetrance in male carriers and 42% in females, suggesting sex as an important modifier of penetrance [6]. Other factors influencing penetrance likely include additional genetic factors, epigenetic factors [89], environmental factors, and gene × environment interactions. Explicit testing of oligogenicity for rare diseases, or gene-environment interactions, require much larger cohorts than those currently available for PAH. However, as more putative risk genes are identified and more PAH cases are studied [7,8,15,17], formal tests to assess the contributions of multiple genetic and environmental risks should be included. In the current study, five of the seventeen cases identified with rare deleterious variants in FBLN2 or PDGFD also carry variants in one or two established or recently reported candidate risk genes. For example, participant 12-207 carries variants in FBLN2 as well as ABCC8 and GGCX, and participant W000073 carries variants in PDGFD and TBX4. We acknowledge the possibility that at least some of the variants identified to date may not be causal and that the relative contribution of individual variants requires further investigation. How multiple rare variants interact to affect PAH pathogenesis, penetrance, endophenotypes, or clinical outcomes will require much larger cohorts and will be one of the major aims of future large international consortia.
Our pediatric data indicate that children present with slightly higher mean pulmonary arterial pressure, decreased cardiac output, and increased pulmonary vascular resistance compared to adults at diagnosis. The early age-of-onset and increased severity of clinical phenotypes suggest that there may be differences in the genetic underpinnings. De novo mutations have emerged as an important class of genetic factors underlying rare diseases, especially early-onset severe conditions [34,90], due to strong negative selection decreasing reproductive fitness [91]. Pediatric-onset PAH fits this category of diseases based on the high mortality during childhood [92][93][94][95][96]. Previously, we reported an enrichment of de novo variants in a cohort of 34 PAH probands with trio data [17]. We have now expanded this analysis to 124 trios with pediatric-onset PAH probands and confirmed the 2.45× enrichment of de novo variants in cases compared to the expected rate. Seven of the variant carriers have variants in known PAH risk genes (TBX4, BMPR2, ACVRL1), and three of the APAH-CHD variant carriers have variants in known CHD or CHD-associated risk genes (NOTCH1, PTPN11, PSMD12). We previously reported rare inherited LGD or D-Mis variants in CHD risk genes NOTCH1 (n = 5), PTPN11 (n = 1), and RAF1 (n = 2) carried by APAH-CHD cases [18]. Specific inhibition of the protein encoded by PTPN11 (SHP2) [97], and induction of mir-204 which negatively targets SHP2 [98], improved right ventricular function in the monocrotaline rat model of PAH, suggesting a more general role of PTPN11 in PAH.
At least eight of the other genes with case-derived de novo variants have plausible roles in lung/vascular development but have not been previously implicated in PAH: AMOT (angiomotin), CSNK2A2 (casein kinase 2 alpha 2), HNRNPF (heterogeneous nuclear ribonucleoprotein F), HSPA4 (heat shock protein family A member 4), KDM3B (lysine demethylase 3B), KEAP1 (kelch-like ECHassociated protein 1), MECOM (MDS1 and EVI1 complex locus), and ZMYM2 (zinc finger MIM-type containing 2). A common single-nucleotide polymorphism in MECOM has been implicated in systemic blood pressure [99]. KEAP1 encodes the principle negative regulator of transcription factor NF-E2 p45-related factor 2 (NRF2). The NRF2-KEAP1 partnership provides an evolutionarily conserved cytoprotective mechanism against oxidative stress. Under normal conditions, KEAP1 targets NRF2 for ubiquitin-dependent degradation and represses NRF2dependent gene expression. KEAP1 is ubiquitously expressed and aberrant oxidative stress response in the pulmonary vasculature is a recognized mechanism underlying PAH. Together, our analysis indicates that 15% of PAH cases are attributable to de novo variants. A larger pediatric cohort will be necessary to confirm some of these genes via replication and identify additional new genes and pathways that will likely be unique to children and not identifiable through studies of adults with PAH.

Conclusions
We have identified FBLN2 and PDGFD as new candidate risk genes for adult-onset IPAH, accounting for 0.26% and 0.35% of 2318 IPAH cases in the US/UK combined cohort, respectively. We note that five of seven FBLN2 variant carriers also have a diagnosis of systemic hypertension. A few cases carry rare variants in more than one PAH risk gene, consistent with oligogenic nature of PAH in some individuals. Analysis of single-cell RNAseq data shows that the new candidate genes have similar expression patterns to well-known PAH risk genes, providing orthogonal support for the new genes and providing more mechanistic insight. We estimate that1 5% of all pediatric cases are attributable to de novo variants and that many of these genes are likely to have important roles in developmental processes. Larger adult and pediatric cohorts are needed to better clinically characterize these rare genetic subtypes of PAH.

Supplementary Information
The online version contains supplementary material available at https://doi. org/10.1186/s13073-021-00891-1.   Table 5. Burden of de novo variants in pediatric-onset IPAH. Supplementary Table 6. Rare de novo risk variants identified in pediatriconset PAH. Supplementary Table 7. Clinical characteristics of pediatric PAH cases with rare de novo variants.
Authors' contributions WCN, WKC, MWP, and YS had full access to all of the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis.

Availability of data and materials
The datasets used and/or analyzed during the current study are available via contact with the senior authors. For PAH Biobank data, a Confidentiality Agreement with the collaborating Regeneron Sequencing Center grants to Dr. Nichols a nonexclusive, worldwide, irrevocable, perpetual, royalty-free sublicensable license to access and use the genomic data for any and all purposes. Therefore, while the PAH Biobank data are not uploaded to a publicly available database, direct access to the data are granted by the corresponding author on reasonable request who has full administrative access to all of the data. The data from the NIHRBR-RD study have been deposited in the European Genome-Phenome Archive [100]. Data from most of the affected participants in the US/UK combined cohort were included in previous publications from our group [7,8,17,18,27,61,100]. The following scripts are available: association test of rare variants with variable threshold of predicted functional scores [48] and principle component analysis of rare variants from Tabula Muris [58].

Declarations
Ethics approval and consent to participate This study was approved by the Institutional Review Boards (IRBs) of the Cincinnati Children's Hospital Medical Center, the East of England Cambridge South national research ethics committee, Columbia University Irving Medical Center as well as the individual IRBs at each of the Enrolling Centers' institutions. All patients have signed consent forms which are on file at the individual Enrolling Centers. No protected health information (PHI) on any patients enrolled in the PAH Biobank, UK NIHR BioResource, or CUIMC has been forwarded to the data analyzing group. Only the individual Enrolling Centers have the ability to re-contact any of the patients enrolled in the study. All research using these patient samples conformed to the principles of the Helsinki Declaration.

Consent for publication
Written informed consent for publication was obtained from the patients/ participants/legal guardians for minors at enrollment.