{'fraction_non_alphanumeric': 0.042425904230126876, 'fraction_numerical': 0.027095744946032854, 'mean_word_length': 5.145526960784314, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 46, 'https://': 11, 'lorem ipsum': 0, 'www.': 6, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Candidate genes and sequence variants for susceptibility to mycobacterial infection identified by whole-exome sequencing

Inborn errors of immunity are known to influence susceptibility to mycobacterial infections. The aim of this study was to characterize the genetic profile of nine patients with mycobacterial infections (eight with BCGitis and one with disseminated tuberculosis) from the Republic of Moldova using whole-exome sequencing. In total, 12 variants in eight genes known to be associated with Mendelian Susceptibility to Mycobacterial Disease (MSMD) were detected in six out of nine patients examined. In particular, a novel splice site mutation c.373–2A>C in STAT1 gene was found and functionally confirmed in a patient with disseminated tuberculosis. Trio analysis was possible for seven out of nine patients, and resulted in 23 candidate variants in 15 novel genes. Four of these genes - GBP2, HEATR3, PPP1R9B and KDM6A were further prioritized, considering their elevated expression in immune-related tissues. Compound heterozygosity was found in GBP2 in a single patient, comprising a maternally inherited missense variant c.412G>A/p.(Ala138Thr) predicted to be deleterious and a paternally inherited intronic mutation c.1149+14T>C. Functional studies demonstrated that the intronic mutation affects splicing and the level of transcript. Finally, we analyzed pathogenicity of variant combinations in gene pairs and identified five patients with putative oligogenic inheritance. In summary, our study expands the spectrum of genetic variation contributing to susceptibility to mycobacterial infections in children and provides insight into the complex/oligogenic disease-causing mode.


Introduction
Tuberculosis (TB) is a serious infectious disease usually caused by Mycobacterium tuberculosis (M. tuberculosis), but also by M. africanum, M. bovis and M. canettii species (Lawn and Zumla, 2011). About one-quarter of the world's population has a TB infection (WHO, 2020), but most people who are infected with TB bacteria do not develop symptoms. Immunocompromised adults and young children are particularly at risk of progressing to active disease. They also develop more severe forms of TB (e.g., disseminated/milliary TB, TB meningitis). Therefore, the vaccination with attenuated M. bovis Bacille Calmette-Guérin (BCG) substrains is recommended to prevent severe TB disease (Roy et al., 2014). The vaccine is given routinely at birth or during childhood in countries with a high burden of TB and it is generally recognized as a safe and effective measure. Serious complications (or BCG infection) are, however, known to occur in children receiving BCG (Kourime et al., 2016).
The use of next-generation sequencing (NGS) offers a rapid and efficient approach to find disease-causing mutations in affected individuals and to discover new disease genes (King et al., 2018). In this study, we implemented whole-exome sequencing (WES) in an attempt to identify risk genes and rare variants in a cohort of seven case-parent trios and two single pediatric cases with mycobacterial infections (BCGitis and disseminated/generalised TB) from the Republic of Moldova. The functional effect of identified candidate variants was investigated through experimental assessment of splicing.

Patient description
Nine children (six girls and three boys aged 2-24 months) with M. tuberculosis complex infection were recruited through the department of Paediatrics at Municipal Hospital in Chisinau, Moldova, during [2016][2017][2018]. Among these patients, eight were diagnosed with BCGitis and one with generalised M. tuberculosis infection. The diagnosis relied on the assessment of clinical features and for six cases was confirmed bacteriologically and/ or histologically. The clinical, demographic and immunologic patient data presented in this study are based on medical records and summarized in Table 1. Immunologic values were within the reference range for most parameters studied in most patients. We note only decreased levels of IgA in patients P2 and P7 and diminished counts of B-lymphocytes (CD20 + ) in P5. These may indicate an immune system defect. However, overall clinical/ immunologic phenotypes in these cases are inconsistent with those reported for BCG-effected patients with SCID or CGD (Fekrvand et al., 2020). Also, other known PIDs associated with BCG complications, such as GATA2 deficiency, which are even rarer than SCID and CGD, were excluded for the same reasons (i.e., inconsistency with the provided clinical/immunologic phenotypes), and thus a diagnosis of MSMD was suspected.
Both parents were present for the interviews and all of them, except one, denied mycobacterial disease during his/her lifetime. For the only father (P2) who announced TB during childhood, the diagnosis of primary (Ghon) complex of pulmonary TB at the age 12 years has been witnessed/confirmed by the medical record. No parents were consanguineous. DNA from both parents was available for seven probands (P1-P7). Therefore, a total of seven trios and two singletons (P8 and P9) were included in the study.
Three probands were selected based on sequencing results to experimentally validate the effects of the identified variants on RNA splicing. Two of them (P6 and P9) were available for direct  Demultiplexing/conversion of the sequencing data was done using bcl2fastq, version 2.18. The resulting fastq files were aligned to the human reference genome (GRCh37) by bwa mem, version 0.7.17, the local realignment was done using abra, version 2.05. After the variant calling by freebayes, version 1.1.0, the annotation of the resulting vcf file was done by SnpSift/SnpEff, version 4.3t, and vcfannotate, version 1.0.0-rc1.

Variant filtering
The filtering of the resulting variants was performed using GSvar, Version 2018_03-12-gbe02570, following steps and criteria outlined in Figure 1. First, only variants with acceptable quality scores (MQM >20 and QUAL >10) were put through the filtering process. Variants were further annotated with population databases (1000G, ESP, ExAC, plus in-house database of variants, n = 900) and only rare variants (MAF cutoff of ≤0.01) with in-house occurrence ≤20 in heterozygous state and ≤2 in homozygous state were considered. We set up a MAF threshold of 0.01 (in the overall world population), because this value consistently explains the incidence of BCG-itis in the Republic of Moldova-15-20 cases per 30000-35000 newborns, assuming autosomal recessive inheritance of the disease. Of note, we primarily relied on the predicted pathogenicity of the variants rather than their frequencies, although we realize that the lower is the MAF the higher is the expected clinical effect.
Two filtering approaches were used for the analysis of WES data in this study (Figure 1). In seven of nine cases in which both parents were available (family cases), WES analysis was triobased (an affected proband with unaffected parents). In this approach, the variants were filtered according to compound heterozygous, homozygous recessive, de novo dominant and X-linked inheritance models. All detected variants were then manually checked using the Integrative Genomics Viewer (IGV) to confirm their quality. Subsequently, genes were prioritized based on predicted pathogenicity of mutations and/or their biological proximity to the known MSMD-causing genes as described in next subsections.

In-silico prediction
Detailed in silico predictions for identified coding nonsynonymous (missense) mutations were made with tools Align GVGD (Tavtigian et al., 2008), SIFT (Sim et al., 2012), MutationTaster (Schwarz et al., 2010), Provean (Choi and Chan, 2015) and Polyphen2 (Adzhubei et al., 2010) implemented in ALAMUT visual (interactive biosoftware, Version 2.7 rev.1, Rouen, France). The missense variants which were predicted as deleterious by at least three out of five methods were assumed likely as pathogenic. For predicting splicing effects, MutationTaster and the Human Splicing Finder (HSF) software (http://www.umd.be/HSF3/) were used (Desmet et al., 2009), and the overall effect was considered deleterious if at least one tool predicts it to be deleterious. Information from the variant genetics database ClinVar (https://www.ncbi.nlm.nih.gov/clinvar/), which adhere to the guidance of the American College of Medical Genetics (ACMG), was also employed for all types of variants.
The effects of mutations on protein structure were assessed using ExPASy translate tool (https://www.swissmodel.expasy.org/) (Waterhouse et al., 2018). In addition, I-Mutant 2.0 software (http:// folding.biofold.org/i-mutant/i-mutant2.0.html) was used to predict the impact of mutations on protein structure stability based on the change in Gibbs free energy (ΔΔG) (Capriotti et al., 2006); the predictions were performed starting from the protein 3D structure with the temperature of 36°C and pH of 7.0 as default parameters. Finally, the Oligogenic Resource for Variant AnaLysis (ORVAL) platform was used to predict clinical consequences of intergenic interactions (https://orval.ibsquare.be/). Both heterozygous and homozygous rare variants identified through candidate genes and trio-based approaches were analyzed in ORVAL in accordance with recommendations of the developers (Renaux et al., 2019).

Functional annotation and categorization of genes
We used DAVID classification system (https://david.ncifcrf. gov/) and PubMed (https://pubmed.ncbi.nlm.nih.gov/) for information about gene function and their involvement in diseases.
Biological affinity/similarity between candidate genes identified by the trio-based prioritization analysis and 19 known MSMD-causing genes was evaluated using the Human Gene Connectome Server (HGCS) (https://hgc. rockefeller.edu/) (Itan et al., 2014). The level of proteinprotein/gene-gene interaction was also analyzed and graphically visualized as nodes (genes/proteins) and edges (the relationship between proteins) using STRING software (https:// string-db.org/). In this analysis, the novel and known MSMDcausing genes were considered as the seed molecules from which direct and indirect protein-protein/gene-gene interactions were obtained in accordance with the developers' recommendations (Szklarczyk et al., 2019).
The genes were also categorized based on their expression pattern as reported in the GTEx database (GenotypeTissue Expression Project database, https://gtexportal.org/home/). Genes with the higher expression levels in whole blood or immune-related cells (i.e., spleen or EBV-transformed lymphocytes) were given the highest priority in this study. Additionally, RD-Match (http://genetics.bwh.harvard.edu/ rdmatch/) was used to estimate probability of finding mutation by chance in a given gene (Akle et al., 2015).

Whole blood RNA extraction and RT-PCR
Peripheral blood samples were collected via PAXgene RNA tubes (Qiagen, Valencia, CA) and stored at -80°C until RNA extraction. RNA was extracted via PreAnalytiX PAXgene blood RNA kit from PreAnalytiX/Qiagen. DNase I treatment was performed at the same time as the RNA isolation using the reagents from the PreAnalytiX PAXgene blood RNA kit. Quantity and quality of the extracted RNA was verified by the NanoDrop ® ND-1000 spectrophotometer (NanoDrop Frontiers in Genetics frontiersin.org 05 Technologies, Wilmington, DE. United States) and the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA. United States).

GBP2 minigene reporter assay
A minigene splicing assay was performed to verify/validate whether the candidate mutation c.1149+14T>C in GBP2 affected splicing products. Briefly, a 1,157 bp fragment harboring exon 7 (281 bp) and part of the flanking introns (246 bp of intron 6 and 581 bp of intron 7) was amplified by standard PCR with patient genomic DNA as template employing a proofreading polymerase and primers: GBP2S-NotI (5′-AATCGGGCGGCCGCTTGACCTAGTATCTCA GCCCCA) and GBP2AS-BamH (5′-CGATTCGGATCCCCA GAACCTCTGACTGTAGCA). The fragment was cloned into the pSPL3 reporter vector between the NotI and BamHI sites (Life Technologies).
Sanger sequencing was used to evaluate whether the wildtype (WT) and mutant-type (MT) expression vectors had been successfully constructed without errors. The minigenes were transfected in duplicates into HepG2 cells using X-tremeGENE transfection reagent (Roche Life Science) according to the manufacturer's instructions and our previous publications (Palhais et al., 2015;Palhais et al., 2016). Forty 8 hours after transfection cells were harvested and RNA was isolated with QIAzol (Qiagen) and chloroform extraction. cDNA synthesis was performed using SuperScript ® VILO Master mix (Invitrogen). Splicing of GBP2 exon seven was examined by PCR using Tempase Hot Start DNA Polymerase (Ampliqon) and vector specific primers SD6 (5′-TCTGAGTCACCTGGACAACC-3′) and SA2 (5′-ATCTCAGTGGTATTTGTGAGC-3′). The PCR products were separated and verified by electrophoresis in an agarose gel and characterized by direct Sanger sequencing. The relative quantity of the bands corresponding to exon inclusion/exon skipping was estimated by capillary electrophoresis using the Fragment Analyzer ™ (Agilent Technologies, Santa Clara, United States) instrument and software and reported as percent exon skipping (relative to the total intensity of the signal in each lane).

Real-time qPCR
The expression of GBP2 was analyzed by quantitative real time PCR in five RNA samples: two independently isolated RNA samples from the carrier of the mutation c.1149+14T>C (i.e., the father of the patient P2) and three matched control samples. All reactions were run in duplicate on the Roche LightCycler 480 with SYBR-green premix FastStart Essential DNA Green Master (Roche Diagnostics, CH) according to manufacturer instructions employing the following primer pairs: GBP2 forward 5′-CAATTACGCAGCCTGTGGTG-3′, GBP2 reverse 5′-GAGCCTAGAGAGAAGCCGTTT and RPL13A forward 5′-AAGGTCGTGCGTCTGAAG, RPL13A reverse 5′-GAGTCCGTGGGTCTTGAG. RPL13A was used as a normalization control because it is stably expressed across different tissue types, including blood cells (Ledderose et al., 2011;Sabatino et al., 2013;Xiao et al., 2017;Dembic et al., 2019).

Statistical analysis
Data from RT-qPCR are expressed as mean ± standard error (s.e.m). Samples were normalized to RPL13A. Means of different groups were compared and analyzed using the Student's t-test.

Summary of exome sequencing data
Overall, 100 M raw reads per sample (range 76M-229 M) were generated, of which 98.7% (range 97.8%-99.5%) and 81.1% (range 74.11%-85.93%) were aligned to the human reference genome (hg19) with a mapping qualities of Q20 and Q30, respectively (Supplementary Table S1). On average, 97.9% of the targeted regions (range 95.3%-98.8%) were covered at least 10 times, and 97.0% of the targeted regions (range 90.6%-98.5%) were covered at least 20 times. The average read depth in the targeted region was 129x (range 96.7x-244.9x) reads per base. Above all, these data fully reflect the reliability of our sequencing, thus providing a solid basis for the follow-up analysis.
Two approaches (trio-based exonic screening and single cases candidate genes analysis) were employed to identify disease-associated variants (Figure 1). The filtering query results for each sample are given in Supplementary Tables S2, S3.

Variants in MSMD and PRR gene panels
Twelve qualifying low frequency variants within eight known MSMD-causing genes were identified in six of the nine patients investigated, including six missense, three splice site/region (±20 bp), two deep intronic and one synonymous mutations (Table 2 and Supplementary Table S3). All detected mutations were in a single heterozygous state. Three of them were found in genes STAT1 (c.373-2A>C in P9 and c.541+50T>C in P7) and IFNG (c.40G>A/p.(Val14Met) in P3), responsible for both autosomal dominant and autosomal recessive forms of mycobacterial disease. Three variants were found in the IL12B gene in two patients, with variants c.877A>G/p.(Lys293Glu) and c.89-14T>C present in the same patient (P6) in a compound heterozygous state (according to the trio-based analysis). Noteworthy, five variants, namely TYK2 c.2441C>T, TYK2 c.157G>A, IFNGR1 c.40G>A, IL23R c.257G>A, IL12B c.961G>A and IL12RB1 c.16196C>T, have been previously reported in patients with MSMD, but it is unclear whether they are actually connected to the disease (https://www.ncbi. nlm.nih.gov/clinvar/).
The study also identified twelve variants in nine PRR genes in six patients (Supplementary Table S4). All were single heterozygous changes, inherited (as seen from the caseparents trio analysis), reported in databases with little or no functional consequences (i.e., intronic, synonymous, benign missense variants).

Trio-based exonic screening of rare variants and phenotype-based gene prioritization
Using trio-based WES, a total of 455 high-confidence rare variants in 201 genes that fit the proposed inheritance criteria (i.e., de novo, homozygous, hemizygous, compound heterozygous variants) were detected in seven trios. These variants included 389 compound heterozygous mutations in 135 genes, 20 homozygous (6 missense, four splice region intronic, one UTR, three intronic, five synonymous, one in non-coding transcript), 32 X-linked hemizygous (12 missense, one in-frame insertion, four splice region intronic, seven UTRs, two intronic, six synonymous), and 14 de novo (4 missense, one frame shift, one in-frame deletion, one splice region, two intronic, five synonymous) variants. After in-silico assessment of pathogenicity, 20 most promising variants in 13 genes were finally selected as top candidates for susceptibility to mycobacterial disease (Tables 3 and 4; Supplementary  Table S5).
By analysis of biological distances between the identified 201 genes and MSMD core genes using HGCS software, two genes-GBP2 and IFNW1 were ranked as the closest to the MSMD genes. Two heterozygous variants in GBP2 (c.412G>A/ p.(Ala138Thr) and c.1149+14T>C) and one homozygous variant in IFNW1 (c.58G>A/p.(Gly20Arg)) were therefore prioritized (Table 3). Similar results were obtained by analysis of gene network using String software (Figure 2). In this analysis, genes GBP2 and IFNW1 are closely interconnected with other known MSMD-associated genes, forming a densely connected module, which implies their direct functional involvement in the disease-associated (i.e., IFN-γ signaling) pathway. Also, KDM6A showed some affinity with the MSMDassociated gene TBX21. Other genes found in this study locate outside the network, suggesting either their affiliation to different pathways or indirect interactions via yet unknown elements.
The significance of new candidate genes was also assessed by their expression in immune related tissues using the GTEx database. Of 15 prioritized in this study genes, four genes, namely GBP2, HEATR3, PPP1R9B and KDM6A, exhibited an elevated transcriptional level in whole blood, spleen or EBVtransformed lymphocytes (Supplementary Table S6), thus lending additional support to their role in increasing susceptibility to mycobacterial infections. Remarkably, genes HEATR3 and PPP1R9B harbor de novo variants, both with predicted damaging effects-a frameshift (Lys132Argfs*3) and a disruptive in frame deletion (Pro252_Pro254del), respectively. Assessment with RD-Match did not reveal true-positive matches between immunodeficiency and either our or known MSMDassociated genes (Supplementary Table S7), thus indicating its insufficient power for mitigating false-positive findings in this study.

Analysis of mRNA transcripts
Two patients were available for the analysis of mRNA, including one with BCG-tis (P6) and another with disseminated TB infection (P9). For the four candidate mutations identified in a patient with BCGitis, effects on splicing were predicted for two compound heterozygous variants c.89-14T>C and c.877A>G in the known MSMDassociated gene IL12B and one de-novo frameshift mutation c.395_396delAA in HEATR3. To examine the functional impact of these variants on mRNA splicing, we performed RT-PCR experiments on total RNA extracted from peripheral leukocytes. Agarose gel electrophoresis showed only single wildtype bands across four loci in all samples (the carriers and controls), indicating no visible effect of these mutations on the pre-mRNA splicing (Supplementary Figure S1).
A novel heterozygous STAT1 intron five splice site mutation, c.373-2A>C, abolishing the splice acceptor site of exon 6, was found in a patient with disseminated TB infection (P9). The mutation causes a dramatic reduction in the MaxEnt score for the mutant 3' splice site (http://hollywood.mit.edu/burgelab/ maxent/Xmaxentscan_scoreseq_acc.html).
RT-PCR with primers located in exons five and seven of STAT1 revealed three products in the mutation carrier sample, while only one middle-length band was observed in all samples used as controls ( Figure 3A). Sanger sequencing of the respective PCR products confirmed that the middle band (205 bp) corresponds to a wild type variant, the upper band (385 bp) corresponds to aberrant splicing where an upstream splice acceptor site is used, resulting in retention of 174 bp of intronic sequence, and the lower band (115 bp) corresponds to the product in which exon six is skipped ( Figure 3B). Thus, we conclude that the mutation STAT1 c.373-2A>C affects the pre-mRNA splicing in vivo ( Figure 3C).

Functional characterization of GBP2 c.1149+14T>C variant
Patient P2 was found to be compound heterozygous for two mutations in the GBP2 gene, namely, a missense mutation c.412G>A/p. (Ala138Thr) and an intronic mutation c.1149+14T>C immediately downstream of the 5′ splice site of exon 7. As RNA was not available from this patient, in vitro functional experiments (i.e., mini-gene assay and SPRi Workflow of WES analysis in 7 case-parental trios (A) and nine single cases (B). The filtered candidate rare variants in single-case and trio-based approaches are shown in Table 2 and Table 3, respectively. See also supplementary material for filtering query results in each trio/single case (Supplementary Tables S2, S3). MQM = mean mapping quality; QUAL = variant quality.
Frontiers in Genetics frontiersin.org 08 analysis) were designed and used for testing the effect of c.1149+14T>C on splicing. In the mini-gene assay, two constructed minigenes, covering exon seven of the GBP2 gene and carrying either a T (wild type) or a C (mutant) at position c.1149+14T>C, were transfected into HepG2 cells and splice products were subsequently analyzed by RT-PCR ( Figure 4A). The analysis showed a slightly decreased exon seven inclusion in HepG2 cells transfected with the MUT minigene as compared to the WT minigene. Quantitative analysis of exon seven inclusion/ exclusion showed a statistically significant difference (p = 0.0028) between WT and MUT lanes, suggesting that c.1149+14T>C has an effect on splicing ( Figure 4A). In particular, the mutation may increase the ability of negative splicing regulatory proteins, such as heterogeneous nuclear ribonucleoproteins (hnRNPs), to bind the pre-mRNA. To investigate this, the kinetics of RNA-protein interactions was analyzed by SPRi using recombinant hnRNP A1 protein and biotin-labeled wild-type or mutant (i.e., harboring c.1149+14T>C) oligonucleotides. The analysis revealed binding of the splicing inhibitory protein hnRNP A1 to oligonucleotides with wild type sequence and the binding was increased for the oligonucleotides carrying mutant sequence ( Figure 4B). Thus, the results from the SPRi analysis supports the effect of the c.1149+14T>C variant on the exon seven skipping via increased binding of hnRNP A1. The father of the child (P2) who is a heterozygous carrier of the c.1149+14T>C variant was available for in vivo mRNA analysis of GBP2. Reverse transcriptase PCR (RT-PCR) analysis of splice variants showed no structural difference in the cDNA between the mutant and control samples for c.1149+14T>C, with only one wild-type band detected ( Figure 4C). The absence of alternative transcript variants and inconsistency with the results from in-vitro mini-gene assay may be due to nonsense mediated decay (NMD), which recognizes and eliminates abnormal mRNAs making them virtually undetectable. To test this hypothesis, we compared the levels of GBP2 expression in the carrier of the c.1149+14T>C mutation (i.e., the father of the affected child) with the GBP2 expression levels in three matched healthy controls. As is shown in Figure 4D, the mutation carrier had significantly lower GBP2 expression compared to controls (p = 0.001). This finding is consistent with the hypothesis of mRNA degradation by NMD and supports the impact of the intronic mutation c.1149+14T>C in GBP2 on alternative splicing.

Protein structure prediction
Swiss ExPASy (Expert Protein Analysis System) Molecular Biology Server was used to predict the effects of identified mutations in HEATR3, STAT1 and GBP2 genes on the protein structure. The de novo deletion (-AA) in the third exon of HEATR3 is a frameshift mutation predicted to cause a severely truncated non-functional HEATR3 protein  For ALIGN GVGD, Class C0 indicates that a change is unlikely to be pathogenic, whereas class C65 represents the highest likelihood of a change to be pathogenic.
Frontiers in Genetics frontiersin.org  (Supplementary Figure S2). The skipping of exon six from the STAT1 gene transcript leads to an in-frame deletion of 90 nucleotides and the removal of 30 amino acids from a stretch linking the N-Terminal to the Coiled-Coil domain of STAT1, while the retention of a 174 bp intronic fragment leads to an out-of-frame change and generates a premature termination codon at position 191 ( Figure 3D). The missense mutation c.412G>A/p.(Ala138Thr) in GBP2 gene is located in an essential functional domain and was predicted to be deleterious by four out of five in silico tools used (Table 2). Noteworthy, the alanine (Ala) at position 138 is one of four hydrophobic residues located along one side of the α-helix ( Figure 5). While the structural analysis revealed no essential differences between the wildtype and mutant models, the substitution of Ala to the more hydrophilic threonine (Thr) leads to the reduction in hydrophobicity of the whole helix from -0.68 to -0.92 (QQAMDQLHYV vs. QQTMDQLHYV, http://www. gravy-calculator.de), that may affect the stability of the protein. Indeed, the predicted folding free energy change determined using the I-Mutant2.0 software was negative (DDG = -1.64) with a sufficient reliability (RI) of 7.0, indicating a destabilizing effect of the above substitution on the protein structure and possibly function.

Identification of an oligogenic etiology in patients with mycobacteriosis
In seven out of nine patients investigated (P1, P2, P3, P5, P6, P7 and P9) the identified variants occurred concomitantly in two or more different candidate genes (Tables 2 and 3). We next used the ORVAL platform to explore their potential digenic or oligogenic impact (i.e., epistatic interactions) in increased predisposition to mycobacterial disease. The analysis identified five patients with at least one pathogenic combination (cutoff >0.532 for CS and ≥50% for SS) between genes ( Figure 6). In particular, the variant combinations in the gene pairs H2BW2-LRP1B, H2BW2-KDM6A (both in P1), GBP2-TYK2 (in P2), SRPX2-ZCWPW1 and ZCWPW1-RBMXL3 (both in P5) were predicted as disease-causing candidates with 95% confidence, while the variant combinations in the gene pairs GAL3ST2-HMCN1, GAL3ST2-TTN and HMCN1-TTN (all in P7) were within the 99% confidence zone for pathogenicity.

Discussion
In this study, the genetic etiology of susceptibility to mycobacterial infections was explored in nine pediatric cases by whole-exome sequencing, bioinformatic analysis and functional assessment. As result, several potentially  For ALIGN GVGD, Class C0 indicates that a change is unlikely to be pathogenic, whereas class C65 represents the highest likelihood of a change to be pathogenic. Frontiers in Genetics frontiersin.org 13 causative rare variants (mostly SNVs) in the MSMD-known and newly discovered genes were identified via candidate gene (singleton WES) and child-parents trio approaches.
Using singleton WES, 12 heterozygous (monoallelic) mutations in the eight MSMD-causing genes, namely TYK2, ZNFX1, IFNGR1, IL23R, JAK1, IL12B, STAT1 and IL12RB, were   STRING interaction (network) analysis (https://string-db.org/) between 34 genes/proteins. Each circle represents a protein/gene: red nodes are known MSMD-causing genes and the blue nodes are novel candidate genes identified in this study via trio-based approach. Numbers inside the nodes denote codes for patients with mutations in the corresponding genes. Interactions between proteins/genes are indicated by lines with variable thickness showing the strength of data support. The confidence cutoff for showing interaction links has been set to 'medium' (0.400).
Frontiers in Genetics frontiersin.org 15 detected in six patients. As mutations in STAT1 and IFNGR1 are known to be responsible for both AD and AR forms of mycobacteriosis, the heterozygous variants identified in these genes  are the most promising causal candidates in this study.
The novel splice-site mutation (c.373-2A>C) reported here in the gene STAT1 in a patient with disseminated TB leads to both partial intron five retention and exon six skipping. The retention of the intronic fragment leads to a frameshift, resulting in early termination of translation and a nonfunctional protein. The in-frame exon skipping event causes a deletion in the protein which would lead to the loss of connection between the N-Terminal and Coiled-Coil domains, which are known to play a role in the dimerization and DNA binding (Mertens et al., 2006). Variants with similar molecular outcomes have been reported previously in several cases with mycobacterial infections (Kong et al., 2010;Kristensen et al., 2011;Vairo et al., 2011;Sakata et al., 2020). These however were in biallelic homozygous or monoallelic compound heterozygous states and were inherited according to an AR mode, while in our case the finding of only a single monoallelic variant is more consistent with an AD inheritance pattern. A simple haploinsufficiency at the STAT1 locus may provide an explanation for our case, as the dosage of normal STAT1 protein generated by the single wild-type allele could be insufficient for complete protection against the highly virulent M. tuberculosis infection in this unprotected (i.e., BCG-unvaccinated) infant. Alternatively, a compound combination with an additional monoallelic variant located in a deep intronic or regulatory region of STAT1 not captured by our analysis may be casual. This second scenario may also be a cause for BCGitis in five patients with mutations in the TYK2, IL23R, JAK1, IL12B and IL12RB genes responsible for AR MSMD. Further genetic and functional studies should be carried out to evaluate the exact biological effects and inheritance pattern of the identified mutants. Regarding IFNGR1, although the identified variant IFNGR1 c.40G>A/ p.(Val14Met) has a relatively high allele population frequency (0.002-0.003) in Europe, it was proven in previous studies to be damaging and was enriched in patients with eczema herpeticum (ADEH+)-a rare and serious skin infection caused by the herpes simplex virus (Gao et al., 2015).
In a trio approach, 25 variants in 15 unknown genes were bioinformatically prioritized. We then critically revised these genes by examining their expressions in immune related tissues. As a result, genes GBP2, HEATR3, PPP1R9B and KDM6A, showing elevated expressions in whole blood, spleen or EBV-transformed lymphocytes, were selected as top-ranked candidates for susceptibility to BCG infection in this study. Of these new candidate genes, GBP2 displays the closest biological relationship with the known MSMD-causing genes. The encoded by GBP2 protein (Guanylate Binding Protein 2) is a member of interferon-gamma induced guanylate-binding proteins (Ramsauer et al., 2007;Yan et al., 2016). Following induction, GBPs confer protective immunity against a variety of microbial pathogens (Kim et al., 2011;Praefcke, 2018;Tretina et al., 2019). Previous experimental studies using animal models and also functional gene-expression studies demonstrated its important role in the protection against mycobacteria (Berry et al., 2010;Kim et al., 2011;Marquis et al., 2011;Tretina et al., 2019;Perumal et al., 2021). In this study, two compound heterozygous mutations in GBP2 gene (c.412G>A/p.(Ala138Thr) and c.1149+14T>C) in a patient with BCG infection (P2) were identified. The maternally inherited missense mutation p.Ala138Thr was predicted to be deleterious by the majority of in silico tools, possibly having an impact on the protein folding and stability. Functional analysis of the second, paternally inherited intronic mutation c.1149+14T>C, demonstrated that this variant is likely to affect splicing of the GBP2 mRNA transcript, leading to decreased levels of functional transcript.
A recently discovered MSMD-causing gene SPPL2A was an outlier compared to the other MSMD-causing genes in our gene networking and connectome server analyses (Figure 2). Likewise, prioritized novel genes in the present study, although not directly Close-up view of residues surrounding Ala138Thr in the 3D structure of the GBP2 protein. Substitution of Ala with more hydrophilic Thr reduces the stretch of hydrophobicity along one side of the α-helix.
Frontiers in Genetics frontiersin.org 18 linked to the core MSMD-causing genes, may also be involved in the disease pathways. Indeed, HEATR3 is a component of the NOD2-mediated NF-kappaB signaling, which plays a crucial role in the control of Mycobacterium within host macrophages (Deng and Xie, 2012;Zhang et al., 2013). Here, we identified a deletion frameshift mutation c.395_396delAA/p.(Lys132Argfs*3) in HEATR3 in a patient with BCG-related complications. The variant introduces an early premature termination codon, which could either be associated with no protein due to nonsense-mediated mRNA decay or a severely C-terminally truncated protein missing all functional domains. Moreover, the de novo status of this variant in our patient increases it is likelihood for being a real cause of BCGitis. Also, undescribed de novo disruptive in-frame deletion was found in PPP1R9B, which encodes a scaffolding protein involved in multiple signaling pathways and cytoskeletal rearrangement in NK cells (Meng et al., 2009).
Other genes in this study, including TTN and LRP1B, are secondary candidates, because of their insignificant expression in immune cells, although their impact cannot be completely excluded.
With the increasing use of next generation sequencing, it has become clear that oligogenic inheritance is more common than previously thought for hereditary disorders, including PIDs Gene networks showing the epistatic interactions between the genes in which rare variants have been identified. Nodes are individual genes and edges represent pair-wise interactions between genes/variants. Numbers inside the nodes denote number of variants. Yellow, blue and light-green nodes represent heterozygous, homozygous and hemizygous variants, respectively. Edge color intensity is proportional to the pathogenicity classification score (CS), whereas edge thickness is proportional to the pathogenicity support score (SS) for a combination of variants. In case two variants within one gene (i.e., compound heterozygotes), the highest pathogenicity scores computed for variant combinations are depicted. Only disease-causing variant combinations/networks (cutoff >0.532 for CS and ≥50% for SS) are illustrated.
Frontiers in Genetics frontiersin.org (Rigaud et al., 2011;Chong et al., 2015;de Valles-Ibáñez et al., 2018). Here, we identified five patients with potential oligogenic inheritance. Noteworthy, two of these cases carried trans heterozygous (monoallelic) variants in two genes belonging to the IL-12/IL-23/IFN-γ axis, namely, IL12RB1 and STAT1 in the patient P9, and TYK2 and GBP2 (a novel candidate-gene) in the patient P2. Since IL12RB1, STAT1 and TYK2 separately are known to cause autosomal recessive MSMD and GBP2 is a novel plausible candidate for MSMD, reduced activity of both proteins in the same close interacting network is a likely disease model. This disease mechanism, termed 'synergistic heterozygosity' (Vockley et al., 2000), could also explain differences in the penetrance of mycobacterial disease in some families and the abovementioned inconsistency with previously published results in the inheritance of STAT1 c.373-2A>C. Also, H2BW2 (a histone protein) and KDM6A (a histone H3 demethylase protein) belong to the same chromatin remodeling pathway regulating gene expression. Although experimental validation of the oligogenic model still needs to be conducted, our data provide insight into the complex disease-causing mode of mycobacterial disease and suggest the combined effect of mutations in IL-12/IL-23/IFN-γ mediated and other pathways genes as an important mechanism in the susceptibility to mycobacterial infections.

Conclusion
The analysis presented here expand the spectrum of genetic variation underlying mycobacterial infections in children. Indeed, we identified 12 heterozygous variants, two of which were novel, in eight known MSMD-causing genes and found pathogenic or likely-pathogenic variants in 15 new candidate genes. Besides this, our study provides insights into a potential oligogenic disease-causing mode for conferring susceptibility to mycobacterial infections. The obtained results rise several topics for further research, including methodological application of WES to identify novel genes for susceptibility to mycobacterial and other infections, experimental and population-based validations of the identified variants and genes, role of oligogenic inheritance in susceptibility to infectious diseases, haploinsufficiency of STAT1 and other MSMD causing genes, relevance of intronic variants, compound heterozygous inheritance of MSMD, etc. Also, CNV and epigenetic modifications that were not tested are assumed for causality in this study and in general. Revealing novel genes and molecular mechanisms underlying mycobacterial infections in children might facilitate the development of therapeutic strategies to prevent and treat mycobacterial disease.

Data availability statement
Variants with functional assessments in this study are deposited in the ClinVar repository (https://www.ncbi.nlm.nih. gov/clinvar/; submission SUB12127021, accession numbers SCV002581932 -SCV002581934). All other relevant data supporting conclusions are included in the article/ Supplementary Material or available by request.

Ethics statement
The studies involving human participants were reviewed and approved by The Research Ethic Committee of the Institute of Phthisiopneumology (Republic of Moldova). Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

Author contributions
TI acquired funding and supervised the study. AV, TI and BA conceived and designed the research. VC and SG collected patient samples and performed clinical assessments. AV, GB and MD performed the experiments. ID and WH provided bioinformatical expertise and support. AV drafted the manuscript. All authors provided critical feedback and helped shape the research, analysis and manuscript.

Funding
This study was supported by the Hannover Unified Biobank. Alexander Varzari was sponsored by the Alexander von Humboldt Foundation (Georg Forster Research Fellowship) and by the Moldovan State Program Particularităţile recidivei tuberculozei pulmonare (theme No. 20.80009.8007.23). The research of Igor Deyneko was carried out within the assignment of Ministry of Science and Higher Education (theme No. 122042700043-9) and development of bioinformatics methods by RSF grant .