Children below 60 months with community-acquired blood-culture positive P. aeruginosa sepsis were eligible. Children with any comorbidities (prematurity, congenital malformations, previous surgery, immunosuppression, known immunodeficiency, and chronic diseases) and children who had been exposed to broad-spectrum intravenous antibiotics were excluded. This study was carried out in accordance with the recommendations of the ethics committees of participating centers (Kantonale Ethikkommission Bern, KEK Ref Nr 029/11, Bern, Switzerland; Mater Health Services HREC13/MHS/4, Brisbane, QLD, Australia). Parents/guardians of all included patients gave written informed consent in accordance with the Declaration of Helsinki. Patients were recruited prospectively and retrospectively in hospital databases and infectious diseases networks. Whenever possible, DNA was also collected from the parents.

Genomic DNA was extracted from whole blood (N = 8), frozen skin biopsy or lymphatic tissue (N = 2), or paraffin-fixed histology slides (N = 1). Exome sequencing libraries were prepared using Agilent SureSelect (V5, 50.4 Mb). Cluster generation was performed using Illumina TruSeq PE Cluster Kit v5 reagents. Libraries were sequenced as 100 bp long, paired-end reads on Illumina HiSeq 2500 using TruSeq SBS Kit v5 reagents.

Sequencing reads were processed using CASAVA v1.82. Reads were aligned to the human reference genome hg19 using BWA (11, 12) v0.6.2. PCR duplicates were removed using Picard 1.27-1 (http://picard.sourceforge.net/).

We used genome analysis toolkit (GATK) (13, 14) version 3.1-1 to call single nucleotide variants (SNVs) and small insertion and deletions (indels) from duplicate-marked bam files. We used HaplotypeCaller for multi-sample variant calling on all samples followed by GATK best practice to call the variants and included only variants that were flagged as PASS by GATK in subsequent analysis.

We used SnpEff (15) version 4.1B to predict the functional impact of variants. As variants can have several predicted effects, we only considered the most severe effect for each variant. The effects in decreasing order of severity are frameshift and in-frame for indels and stop-gain, splice site-disrupting, non-synonymous, synonymous, intronic, UTR, non-coding exon, and intergenic for SNVs. Variants were annotated as putative loss-of-function (LoF) if they were stop-gain or splice site-disrupting SNVs or frameshift indels mapping to the first 95% of coding region, or larger deletions removing either the first exon or more than 50% of the protein-coding sequence of the affected transcript.

We checked the minor allele frequency (MAF) of all previously described variants in the following datasets: the NHLBI exome sequencing project (ESP, N = 6503) (16), the 1000 genomes project phase 2 (1KG, N = 876) (17) and the UK10K project (UK10K, N = 3621) (18), the Exome Aggregation Consortium database (ExAC, N = 60,706) (19), and a set of 533 in-house control exomes.

We restricted analyses to non-synonymous and LoF exonic variants with MAF <1% in ESP, 1KG, UK10K, ExAC, and in-house control exomes. We analyzed each family separately assuming autosomal recessive, autosomal dominant, and X-linked recessive inheritance models. Only variants with >10× coverage in both the parents and the offspring were included in the de novo analysis. We assumed full penetrance for the potentially causal variants and given the fatality of the phenotype did not consider a mutation as potentially causal if it was present in the parents, any of the abovementioned databases or in the in-house control exomes in the same zygosity form (i.e., heterozygous or homozygous) as in our patients.

We restricted analyses to non-synonymous and LoF exonic variants with MAF <1% in ESP, 1KG, UK10K, ExAC, and in-house control exomes. We analyzed each individual separately assuming autosomal recessive and X-linked recessive inheritance models. We assumed full penetrance for the potentially causal variants and given the fatality of the phenotype did not consider a mutation as potentially causal if it was present in any of the abovementioned databases or in the in-house control exomes in the same zygosity form (i.e., heterozygous or homozygous) as in our patients.

We ran a targeted search for rare (MAF <1%), non-synonymous, and LoF variants in a list of 252 known PID genes. This list includes known genes in 50 PID syndromes compiled by The International Union of Immunological Societies (IUIS) Expert Committee in 2015 (20) and 3 newly published PID genes since the latest release of IUIS till April 2016 (Table 1).

We used the database for annotation, visualization, and integrated discovery (DAVID) (24, 25) v6.7 with default options and highest classification stringency for functional annotation clustering of PID genes carrying rare, non-synonymous variants. We used hypergeometric test for assessing the significance of gene enrichment results.

Serum of a patient heterozygous for a complement C9 variant was diluted with GVB⁺⁺ in a series starting at 1:10 to 1:640 in 1:2 steps. The same was done with a serum pool as a control. This pool consisted of serum from over 300 healthy donors. Forty microliters of these dilutions were pipetted in duplo in a round-bottom 96-wells plate. To these wells, 40 μl of GVB⁺⁺ was added to the samples and the blanc, to the positive control 40 μl of 1.7% saponin was added (100% lysis). Ten microliters of purified C9 (Quidel) or GVB⁺⁺ were added (end conc 50 μg/ml). Next, in all wells, 150 μl of EA’s were added. EA’s are sheep erythrocytes (Hatunalab, Sweden) coated with an optimal dose of Amboceptor (Rabbit-anti-Sheep erythrocyte; Dade Behring) brought to a concentration of 0.25 × 10⁸ cells/ml. This was incubated for 1 h at 37°C while agitating. After incubation, the plate was centrifuged for 5 min at 2000 rpm. Fifty microliters of each supernatant were pipetted to a flat-bottom 96-wells plate and diluted with 150 μl GVB⁺⁺. Extinction was measured at 415 nm using a spectrophotometer. The percentage of lysis was calculated as follows: ((average of sample duplo−average blanc)/(average 100%−average blanc))×100.

GVB⁺⁺ buffer consisted of 2 mM 5,5 di-ethylbarbituurzuur (Genfarma), 1.15 mM Na 5,5 di-ethylbarbituraat (Bufa), 96 mM NaCl (Merck), 0.5 mM CaCl (Merck), and 0.18 mM MgCl₂ (Merck).

We identified 11 previously healthy children with community-acquired P. aeruginosa bacteremia from whom DNA could be obtained. Child–parent DNA trios were available in seven affected families. The presenting age ranged from 6 months to 4 years, and 7/11 (64%) of patients were male (Table 2). All patients had never been admitted to hospital or exposed to intravenous antibiotics prior to presenting with P. aeruginosa sepsis. Nine patients had developed respiratory symptoms within 72 h of leading to hospital admission, one presented with Ecthyma gangrenosum (Figure 1), and eight (73%) died. Three children died in the emergency department shortly after presentation, and all deaths occurred within 48 h of hospital admission. Postmortem examination reports were available in seven cases. In four, P. aeruginosa grew in high quantities from nasal and oral swabs, tracheal secretions, and lungs. Viral co-infections were found in five deceased patients, including human herpes virus-6, parainfluenza virus-3, human metapneumovirus, respiratory syncytial virus, and varicella virus. One surviving patient presented with acute abdomen, and P. aeruginosa bacteremia was thought to result from intra-abdominal perforation. Another surviving patient was diagnosed with urosepsis and P. aeruginosa grew in urine and blood. The third surviving patient presented with recurrent parainfectious neutropenia. Six patients received aminoglycosides and/or anti-Pseudomonas beta-lactam antibiotics during sepsis, including all three survivors.

For each sample, 95% of reads passing filtering criteria were unique (not marked as duplicate); 97% of unique reads could be aligned to the human reference genome hg19. The mean on-bait coverage was 66×, with 95% of target bases achieving at least 10× coverage and 76% achieving at least 30× coverage (Table S1 in Supplementary Material). In 118 of the PID genes, at least one exon had an average coverage of <2×. Further investigation of these low-coverage intervals showed that the majority of them are overlapping with untranslated regions (UTRs).

A total of 115,604 SNVs and 10,814 indels passed GATK quality control, including 28,795 synonymous variants, 28,284 non-synonymous variants, 713 in-frame indels, 660 frameshift indels, 283 stop-gained, and 151 splice-site variants (Tables S2 and S3 in Supplementary Material).

Thirty-eight rare, non-synonymous coding variants (MAF <1%) were observed in the 11 patients after testing for different inheritance models (Table 3), among which 12 were never seen in the same zygosity form (heterozygous or homozygous) in publicly available databases (Table 4). In addition to the above 38 variants, we performed a targeted search for rare, non-synonymous coding variants in 252 previously known PID genes and found 76 variants in 61 genes (Table 5). Functional classification of these 61 genes using DAVID highlighted the complement pathway as the most enriched cluster (DAVID enrichment score: 15.73). In total, 14 of the 28 complement pathway genes present in the list of 252 known PID genes carried at least one rare, non-synonymous coding variant in our cohort (hypergeometric probability: p(x > 14) = 0.003).

In one patient (patient S4) who died of fulminant septic shock, we identified a novel single-base insertion on the X chromosome (C > CT, chromosome 20, position 100617191) leading to a frameshift in the Bruton agammaglobulinemia tyrosine kinase (BTK) gene. Other than a middle ear infection during infancy, the previous history, growth, development, and vaccination history of the patient had been unremarkable. Taqman genotyping and clinical genetic testing confirmed presence of mutation in the deceased patient and revealed a de novo occurrence in the mother, while other family members were healthy. BTK plays a crucial role in B-cell development. Postmortem serological testing confirmed the absence of immunoglobulins, consistent with BTK LoF: IgG < 1.3 g/l (4.22–11.9); IgA < 0.2 g/l (0.2–1.58); IgM < 0.2g/l (0.48–1.9).

A novel homozygous deletion spanning 6 bp on chromosome 20 (ACTCGAG > A, X chromosome, position 31393171) was observed in a 2-year-old boy with consanguineous parents (patient S13), leading to an in-frame deletion in a conserved region of the catalytic domain of the DNMT3B gene. Previously described missense mutations in the same protein domain are known to cause immunodeficiency-centromeric instability-facial anomalies (ICF) syndrome (26), a rare disease characterized by variable immunodeficiency and recurrent infections with mild facial abnormalities. The fatal septicemia was preceded by uncomplicated varicella, and the patient had a history of recurrent mild respiratory and upper airway infections that had been attributed to poor health conditions. The patient had normal T and B cell counts but reduced immunoglobulin levels, consistent with ICF syndrome: IgG 0.16 g/l (4.22–11.9); IgA <0.06 g/l (0.2–1.58); IgM 0.05 g/l (0.48–1.9).

We also observed a rare, known (MAF = 0.001 in ExAC) heterozygous stop-gain variant in the complement C9, on chromosome 5 (G > T, chromosome 5, position 39342214) in a male patient (patient S15) with severe Pseudomonas septicemia and recurrent parainfectious neutropenia. The patient survived without sequelae. Complement reconstitution assay with classical or terminal pathway proteins (C1–C9) showed normal lytic activity (Figure 2), suggesting that this variant is unlikely to be causal to increased susceptibility to P. aeruginosa.