Rare-variant collapsing analyses identified risk genes for neonatal acute respiratory distress syndrome.

Background: Acute respiratory distress syndrome (ARDS) could account for a considerable proportion of neonatal death, while the genetic etiology and pathophysiology of neonatal ARDS remain elusive. In this case-control study, 515 neonates were enrolled in the China Neonatal Genomes Project (CNGP, NCT03931707) from August 2016 to June 2021, including 196 ARDS and 319 non-ARDS matched by sex, gestational age, birth weight, perinatal asphyxia, pneumonia, sepsis, and necrotizing enterocolitis. Clinical exome sequencing was used to detect genetic variants. Collapsing analyses together with permutation tests were used to identify ARDS risk genes enriched for rare variants in ARDS samples. In silico functional interaction analysis and expression pattern studies at different stages of lung development were used to investigate the biological functions of the risk genes.
Results: Collapsing analyses identified that rare variants were significantly abundant in the genes associated with the precursor of the lamellar body and there were eight predicted risk genes with strong confidence (P < 0.01). Among them, the expression of EDNRB increased significantly in lung development and was up-regulated in ARDS (P < 0.05). In addition, 151 predicted transcriptional target proteins of EDNRB were highly enriched in the lamellar body responsible for pulmonary surfactant storage and secretion. Conclusions: In our study, the genes associated with pulmonary surfactant storage and release were highly enriched with rare variants. A novel neonatal ARDS risk gene EDNRB may be a key gene for neonatal lung development and pulmonary surfactant homeostasis. Additional validation in independent patient populations and further exploration of underlying molecular mechanisms are warranted.

Introduction

Neonatal acute respiratory distress syndrome (ARDS), which has a prevalence of 1.5 % and overall mortality of 17 %–24 % in the neonatal intensive care unit [1], is characterized by extensive pulmonary inflammation and surfactant catabolism leading to pulmonary dysfunction [2]. It is considered to be triggered directly or indirectly by a range of etiologies [3], such as meconium aspiration syndrome, perinatal asphyxia, necrotizing enterocolitis, sepsis, and infectious pneumonia [4], [5]. Surfactant is liquid-protein complexes that line the alveolar air–liquid interface, contributing to alveolar patency and effective gas exchange [6]. However, the lungs and immune systems of neonates are still rapidly developing, and the pathophysiology of respiratory insufficiency and the role of surfactant homeostasis remain unclear.

In addition to clinical risk factors, genetic predisposition has recently been identified as a potential etiology of respiratory distress in neonates [7], [8]. In the Online Mendelian Inheritance in Man database, several genes have been reported to be responsible for surfactant catabolism and possibly cause clinical features of respiratory distress syndrome in term neonates, including SFTPB, SFTPC, and ABCA3. Moreover, transcription factor NKX2-1 was also responsible for the onset of respiratory distress as it participates in the expression, catabolism [9], and function [10] of surfactant by regulating the expression of the above genes [11]. Pathogenic variants of NKX2-1 can lead to choreoathetosis, hypothyroidism, and neonatal respiratory diseases [MIM:610978]. Researchers have found that single nucleotide polymorphisms having a higher allele frequency (AF) in the neonates with respiratory distress than in the control groups, such as SFTPD rs1923537 [12], LPCAT1 rs9728 [13], and NOS3 Glu298Asp[14], may also be another essential type of genetic risk factor for respiratory distress in neonates. However, only a few neonatal ARDS cases could be attributed to the genetic cause mentioned above and no studies have investigated their genetic risk of susceptibility.

To investigate potential genetic risk factors in extreme respiratory phenotypes, genome-wide association studies have been widely used for the common variants with a minor allele frequency (MAF) > 5 % [15], [16], [17]. Recently, empirical evidence suggested that rare variants (MAF < 1 %) may be a novel potential genetic risk factor associated with complex diseases [18]. Genes related to the onset of complex disease or trait variation could be identified by applying rare-variant collapsing analyses [19].

The variants of risk genes can be valuable markers for mapping different pathologies, allowing neonates to benefit from early diagnosis. We hypothesized that a high burden of rare variants may lead to a greater susceptibility of neonates to ARDS. Therefore, we conducted a case-control study to investigate the potential genetic etiologies in neonatal ARDS using rare-variant collapsing analyses and attempted to validate their biological plausibility through expression studies at different stages of lung development as well as in silico functional interaction analysis.

Methods

Study design and participants

This case-control study is based on the retrospective analysis of single-center data from the China Neonatal Genomes Project (CNGP, NCT03931707). In this project, clinical exome sequencing was performed on patients who have the suspicion of genetic diseases and underwent consultation with geneticists. We included 1,019 Han Chinese neonates admitted within 24h of birth at the Children’s Hospital of Fudan University from August 2016 to June 2021. The diagnosis criteria of neonatal ARDS refer to the definition of Montreux[2] and the following points need to be satisfied: 1) acute onset (ie, within one week) from a known or suspected clinical insult but respiratory distress syndrome, transient tachypnoea of the neonate, or congenital anomalies is not a primary current acute respiratory condition; 2) diffuse, bilateral, and irregular opacities or infiltrates, or complete opacification on chest radiographs, which are not fully explained by local effusions, atelectasis, respiratory distress syndrome, transient tachypnoea of the neonate, or congenital anomalies; 3) absence of congenital heart disease explaining the pulmonary edema (this includes ductus arteriosus with pulmonary overflow if no acute pulmonary haemorrhage exists). Echocardiography is needed to verify the origin of edema; 4) oxygenation deficit expressed as oxygenation index (FiO2*Paw*100/PaO2) ≥ 4. Hospitalized neonates who did not diagnose with ARDS and without any type of respiratory support were classified as non-ARDS. Neonates with known genetic syndromes or chromosomopathies were excluded from the study. The neonatal ARDS diagnoses were evaluated by two neonatologists back-to-back.

Patient information was acquired from electronic medical records. The samples used in this study were collected with appropriate informed consent and the approval of the ethics committee of the Children’s Hospital of Fudan University (2016–235).

Gene-based collapsing analyses

This study focused on predicted rare damaging variants (including potentially damaging variants) of 3,203 genes captured by the Agilent ClearSeq Inherited Disease Kit, which included 2,742 confirmed disease-causing genes [20]. The process of gene testing and damaging variant filtering are described in the eMethods in the Supplement. We performed the following steps to prepare the rare variant collapsing analyses. First, each variant was annotated to obtain its mutation type, which was further grouped into four types, including protein-truncating variants (PTVs), missense or non-synonymous variants (MISs), synonymous variants, and non-coding variants, according to the transfer table shown in eTable 1 in the Supplement. Next, the number of variants in each mutation type was calculated and summarized for each gene, and the genes with at least three PTV or MIS variants in the matched samples were included in the analysis. Fisher’s exact test was used to test whether the number of variants from each mutation type was significantly higher in the ARDS group than in the non-ARDS group. An empirical (permutation-based) expected probability distribution was used for quantile–quantile plots (eMethods in the Supplement) to estimate risk genes with more robust statistics. We will consider the gene to be a risk gene if its observed P value is significantly smaller than the expected P value. In particular, tests of PTV and MIS variants were used to identify ARDS-related burden genes. Synonymous variants and non-coding variants were treated as the near-neutral background, and genes with significant differences found at the synonymous or non-coding level were filtered out.

Gene expression profiles and in silico functional interaction analyses

We downloaded the RNA transcriptome data of mouse and human samples from the GEO database. The mouse dataset GSE74243 contains lung tissue samples at 26 time points across all five canonical stages of lung development. The integrated human dataset (GSE14334 and GSE43767) included lung tissue samples from 6 to 24 postovulatory weeks (PWs), which is the initial branching period of airways. GSE14334 contains 38 samples from 53 to 154 days (6–24 weeks) post-conception, representing 29 distinct time points. GSE43767 contains samples from 6 to 8 PW (n = 10) and 16 to 24 PW (n = 9). The pediatric dataset GSE152980 contains 48 ARDS and 25 non-ARDS nasal cytology brushings. Next, we normalized the gene expression profiles of mice and humans to visualize the changes across lung development and analyze the expression differences between ARDS and non-ARDS in pediatrics. Finally, we constructed a lung-specific functional interaction network in silico study using the gene expression profile of human lung tissue (eMethods in the Supplement).

Statistical analysis

The propensity score matching strategy was used to screen the original collected ARDS cases and the non-ARDS controls using the R package MatchIt. Continuous variables were evaluated by the Kruskal–Wallis H or T-tests. Categorical variables were tested by a Chi-Squared and Fisher Exact Test. Gene expression correlation was performed using the Pearson correlation coefficient. The threshold of significance was P < 0.05. All analyses were performed using R software (version 4.0.3, http://cran.r-project.org).

Results

Study population

Altogether, we recruited a total of 1,019 neonates (eFigure 1 in the Supplement), including 374 ARDS and 645 non-ARDS. Their demographic and clinical characteristics were described in eTable 2. To investigate the genetic risk factors for ARDS, we performed propensity score analysis to establish a matching cohort. The matching features included basic information (sex, gestational age, and birth weight) and ARDS-related clinical manifestations (perinatal asphyxia, pneumonia, sepsis, and necrotizing enterocolitis). Based on propensity score matching, we established a cohort of 319 non-ARDS and 196 ARDS neonates (Fig. 1). No significant difference in clinical features was detected between the two groups (Table 1). All of these neonates underwent clinical exome sequencing, and no diagnosable genetic findings were detected in either the non-ARDS or ARDS group after curation by two geneticists in a blinded fashion.

Fig. 1

Study design. Based on a case-control study with matched clinical information, we searched for strong confidence risk genes of neonatal ARDS through clinical exome sequencing and rare variant collapsing analyses. The specific process is divided into five steps. First, 515 samples were sequenced by clinical exome sequencing after propensity score matching by clinical characteristics (step a). Variant calling and quality control were performed by GATK and an in-home pipeline (step b). Ensembl Variant Effect Predictor (VEP) and ANNOVAR were used to annotate all variants, including protein-truncating variants (PTVs) and missense or non-synonymous variants (MISs) (step c). Rare-variant collapsing analyses were applied, and 55 ARDS potential risk genes were identified (P value < 0.05) (step d). Next, the results were evaluated using a quantile–quantile (QQ) plot (step e). Finally, we focused on eight strong confidence risk genes.

Table 1

Clinical information of both ARDS and non-ARDS groups.

Clinical characteristics	All samples			Matched samples
	ARDS	Non-ARDS	P value	ARDS	Non-ARDS	P value
No. (%) with data	374	645		196	319
Sex, Male	126 (33.7)	381 (59.1)	<0.001	107 (54.6)	186 (58.3)	0.462
Gestational age, mean (SD), weeks	34.10 (2.73)	36.53 (1.94)	<0.001	35.92 (1.94)	36.03 (1.88)	0.501
Birth weight, mean (SD), kg	2.40 (0.66)	2.81 (0.67)	<0.001	2.71 (0.63)	2.63 (0.74)	0.220
Perinatal asphyxia	90 (24.1)	68 (10.5)	<0.001	33 (16.8)	42 (13.2)	0.309
Pneumonia	76 (20.3)	29 (4.5)	<0.001	24 (12.2)	28 (8.8)	0.264
Sepsis	135 (36.1)	136 (21.1)	<0.001	69 (35.2)	91 (28.5)	0.136
Necrotizing enterocolitis	33 (8.8)	24 (3.7)	0.001	5 (2.6)	8 (2.8)	1.000

Detection of ARDS burden gene sets and genes

In ARDS cases versus non-ARDS controls, we found that both PTV and MIS variants were significantly enriched in the set of genes associated with the late endosome (ORPTV = 1.85, P = 0.03; ORMIS = 1.22, P = 0.01; Fig. 2a), which is the precursor of the lamellar body [21]. In addition, we identified 55 candidate burden genes with more PTV or MIS variants in the ARDS group compared with the non-ARDS group (eTable 3 in the Supplement) from 171 and 2,050 genes that were included in the PTV and MIS variant collapsing analyses. Five of the genes (MUTYH, ATXN3, HTT, TBP, and DCHS2) had a higher PTV allele frequency, and 50 genes had a higher MIS allele frequency (P < 0.05) in the ARDS group. There was no overlap between these two types of candidate burden genes.

Fig. 2

Evaluation of the gene sets and strong confidence risk genes (SCGs) for ARDS. (a) The significantly enriched gene sets in rare variants. “N Genes Hit by Variants” refer to the number of genes in the set and the number of genes that had at least one rare variant hit, respectively. CI, confidence interval; OR, odds ratio. (b) Differential expression of eight SCGs in the pseudoglandular and canalicular stages of human lung development. (c) Expression patterns of eight SCGs in the five canonical stages of mouse lung development (EMB, embryonic; PSG, pseudoglandular; CAN, canalicular; SAC, saccular; ALV, alveolar). (d) Differential expression of EDNRB in the nasal cytology brushings of pediatric ARDS and non-ARDS. Expression levels were normalized for processing. * P < 0.05; ** P < 0.01; *** P < 0.001.

The permutation-based expected probability distribution helped us identify the significant reference thresholds (P values) of the statistics. For 55 candidate burden genes, P < 0.01 indicates a strong level of confidence, and 0.01 < P < 0.05 indicates a moderate level of confidence (eFigure 2a-b in the Supplement). The PTV burden genes all exhibited moderate confidence levels. Eight of the MIS burden genes (LDLR, FANCA, NHS, SHROOM4, EDNRB, COL9A1, CNGB3, and ATP6V0A1) exhibited strong confidence levels (Table 2, eFigure 2b in the Supplement), in which the four variants (LDLR rs13306511, FANCA rs56369086, SHROOM4 rs189694750, and EDNRB rs2070591) with high allele frequencies (AFs > 5/515) in CNGP did not show genetic heterogeneity (P > 0.05, eFigure 3a-d). Finally, we considered the eight MIS burden genes with strong confidence as strong confidence risk genes (SCGs).

Table 2

The number of variants in eight SCGs and the results of rare-variant collapsing analyses.

Gene	ARDSPTV	ARDSMIS	Non-ARDSPTV	Non-ARDSMIS	PTV		MIS
					P	adjP	P	adjP
LDLR	0	14	0	2	–	–	5.21E-05***	0.107
NHS	0	13	0	4	–	–	1.25E-03**	1
SHROOM4	0	20	0	11	–	–	1.91E-03**	1
FANCA	2	29	2	21	4.91E-01	1	2.11E-03**	1
EDNRB	0	13	2	5	–	–	2.97E-03**	1
COL9A1	1	16	1	8	–	–	3.49E-03**	1
CNGB3	0	10	0	3	–	–	4.66E-03**	1
ATP6V0A2	0	7	1	1	–	–	5.83E-03**	1

SCGs, strong confidence risk genes; P value is the significance of the Fisher’s exact test (one-sided) and adjP is the adjusted p value with Bonferroni correction; * P < 0.05; ** P < 0.01; *** P < 0.001.

Expression signature of strong confidence risk genes

To detect the potential biological function of these SCGs, we detected their expression level in human and mouse lung tissue. Human lung development at the early embryonic stage can be divided into two periods: the pseudoglandular stage (branching morphogenesis, 6–16 weeks) and the canalicular stage (formation of epithelial sacs and appearance of capillaries, 16–24 weeks) [22]. We found that the expression levels of three SCGs (EDNRB, SHROOM4, and CNGB3) significantly increased from the pseudoglandular stage to the canalicular stage (P < 0.05, Fig. 2b), and the expression levels of the other five genes (LDLR, FANCA, NHS, COL9A1, and ATP6V0A2) progressively decreased (P < 0.01, Fig. 2b). Among them, EDNRB had the highest relative expression compared with other SCGs (Fig. 2b). Similar results were shown in the transcriptional profiles of mice at two identical periods of lung development and Ednrb gene expression increased rapidly from embryonic to saccular stages in the prenatal mice and peaked at birth (Fig. 2c). In addition, EDNRB expression was significantly upregulated in patients with pediatric ARDS compared to the non-ARDS group (P = 0.035, Fig. 2d). This finding demonstrates that EDNRB has a more important role in the whole process of lung development and in the disease onset.

Functional effects of variation on the in silico functional interaction network

In the in silico functional interaction network, EDNRB was estimated to positively regulate 101 target genes (Fig. 3a) and negatively regulate 50 target genes (Fig. 3b) at the transcriptional level. In addition, we found that EDNRB may positively regulate SFTPB and SFTPC expression (Fig. 3a), known as surfactant proteins. Moreover, EDNRB was highly correlated with the expression of SFTPB and SFTPC (R = 0.9, P < 2.2e-16, Fig. 3c). The enrichment analyses of the 151 target genes found that the lamellar body, the surfactant storage organelle in type II epithelial cells, was the top item in Gene Ontology analyses (P = 1.74e-7, Fig. 3d). These results may indicate the potential important function of EDNRB in the onset of ARDS.

Fig. 3

EDNRB targets genes and their functions. (a-b) The 101 positively regulated genes (a) and 50 negatively regulated genes (b) of EDNRB were predicted by in silico functional interaction analysis. The thickness of the line represents the Spearman correlation. (c) Correlation of EDNRB with SFTPB and SFTPC expression in two stages of human lung development. (d) The top ten items of Gene Ontology enrichment analyses of the 151 EDNRB target genes were sorted by P value. The rich factor is the number of genes belonging to this term out of 151 NDNRB target genes divided by the number of genes belonging to this term out of 16,648 genes included in the human lung expression profile data.

Discussion

Neonatal ARDS is a secondary deficiency of pulmonary surfactant induced by multiple inflammatory causes both inside and outside the lung, in which alveolar epithelial cells and vascular endothelial cells are disrupted by pro-inflammatory factors, resulting in a dysfunctional pulmonary capillary barrier and increased permeability[23]. The risk genes neonates inherit may put them at an increased risk for ARDS and play a role in how the lungs respond to damage. The study of ARDS risk genes helped us to understand the variability in susceptibility to the development of lung lesions. With the widespread application of next-generation sequencing, there are>40 candidate genes associated with the development or outcome of ARDS have been identified, including ACE, VEGF, and SFTPB[24]. However, we found only a marginal significant burden of rare missense variants in SFTPB among the candidate genes (PMIS = 0.056, eTable 4 in the Supplement), which may be due to the different pathogenesis of ARDS at different ages.

The essential attributes characterizing the biology of the lung were identified by analyzing the expression pattern at different stages of lung development[25]. Similarly, we highlighted the importance of EDNRB by outlining the expression patterns of the eight SCGs in the different periods of lung development in humans and mice. EDNRB is a non-selective type endothelin receptor of endothelin-1. Endothelin-1 signaling was found to mediate dynamic changes in fetal pulmonary blood flow [26]. The 151 predicted target genes of EDNRB were mainly enriched in the lamellar body-related gene set, and its precursors were also significantly enriched by rare variants. This organelle is found in type II alveolar cells in the lungs and contains densely packed concentric membrane layers for storing and secreting surfactants [27]. Dysfunction of the lamellar body can lead to neonatal respiratory distress syndrome and interstitial lung disease in children, such as a rare pathogenic variant of ABCA3 [28]. Recent studies have found that lung inflammatory cell infiltration can also cause abnormal lamellar bodies [29]. Therefore, we hypothesized that EDNRB may participate in the surfactant secretion process by regulating the function of the lamellar body and late endosome. Regarding the EDNRB variants identified in this study, we found a missense variant annotated as a disease-causing mutation by the Human Gene Mutation Database and MutationTaster, namely, rs2070591, that was highly enriched in the ARDS neonates (P < 0.05, eTable 5 and eTable 6 in the Supplement). As a risk variant for Hirschsprung's disease [30], it may also be the main cause of EDNRB dysfunction in ARDS. Based on these results, EDNRB may be a potential ARDS risk gene and is involved in surfactant homeostasis.

Because neonatal ARDS is a definition that has been around for a short time, there are few studies have been conducted to identify candidate risk genes. However, several other studies have been published on the investigation of risk genes affecting other similar extreme respiratory phenotypes, such as neonatal respiratory distress syndrome and bronchopulmonary dysplasia. Several investigators have done statistical association analyses of individual patient characteristics with non-common variants, albeit in small cohorts (<50 patients) [31], [32]. Another small cohort study using whole genome sequencing to identify rare variants in the genome identified 258 genes with rare nonsynonymous mutations in a cohort of 50 twin pairs [33]. Several other studies have begun to use whole-exome data from large cohorts of several hundred patients to identify novel genes and associated pathways that play an important role in the development of neonatal lung disease [34], [35]. The relative strength of our study is that we used a large patient cohort and a systematic combination of multi-omics approaches to identify a more accurate and informative set of risk genes. However, there were few overlaps between our result and candidate risk genes from previous studies (overlap = 1.3 %, 7/55, eFigure 2c in the Supplement). This finding may be related to the differences of samples in the regional distribution and phenotypes, including gestational ages, disease types, and races of the study populations.

A major strength of this study is that it was the largest cohort of ARDS neonates to date to conduct a risk genetic study. The approach of using only genetic variation information - does not provide an integrated view. In our study, network-based in silico functional interaction analysis can elucidate the genes and pathways involved, as well as biomarkers and potential drug targets [36]. In silico functional interactions and gene-set enrichment analysis have been used to integrate the vast amount of complex data obtained from genomics and transcriptomics studies [37].

Limitations

Some limitations need to be improved in the present study. Firstly, our results require additional validation in an independent patient population and further exploration of the underlying molecular mechanisms. Secondly, given the limitations of clinical exome sequencing, we acknowledged that the systematic analyses of whole-genome sequencing in patients from different regions or different races would provide more comprehensive results.

Conclusions

In this study, we applied clinical exome sequencing in a cohort of neonatal ARDS and used multiple analytical approaches to identify EDNRB as a possible ARDS risk gene. Although EDNRB was not the most statistically significant gene identified in the initial analysis, we clarified several biologically plausible proofs by gene expression at different periods of lung development and by in silico functional interaction analysis.

Ethics Statement

The study was conducted ethically following the World Medical Association Declaration of Helsinki and approved by the Research Ethics Committee of the Children’s Hospital of Fudan University (CHFudanU_NNICU11). Written informed consent was required from the legal guardian for genetic analysis.

Author Contributions

WZ and LY conceived and supervised the project. HC, XC, and XD designed the study and participated in data analysis and manuscript drafts. LH, XD, and GC participated in the design of the study and critical revision of the manuscript. CY, JZ, and YL contributed to data evaluation and data interpretation. All authors provided comments on the draft manuscript and approved the final manuscript.

Funding/Support

This study was supported by the Shanghai Municipal Science and Technology Project Shanghai (Grant No. 22Y11906800), the Shanghai Municipal Science and Technology Major Project (Grant No. 2017SHZDZX01), and the National Key R&D Program of China (Grant No. 2020YFC2006402).

Role of the Funder/Sponsor

Funders provide financial support for this study and do not participate in study design, data collection, data analyses, interpretation, or writing of the report.

Data Sharing Statement

The datasets generated during and analyzed during the current study are available from the corresponding author on reasonable request.

Additional Contributions

We are very grateful to the patient's families and thank them for their trust in our lab. Thank our colleagues for their invaluable contribution to this project. Thanks for the following foundations: the Shanghai Municipal Science and Technology Project Shanghai (Grant No. 22Y11906800), the Shanghai Municipal Science and Technology Major Project (Grant No. 2017SHZDZX01), and the National Key R&D Program of China (Grant No. 2020YFC2006402).

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.