Identification of potential pathogenic mutations in Chinese children with first branchial cleft anomalies detected by whole-exome sequencing.

IMPORTANCE: First branchial cleft anomalies (FBCAs) are rare congenital malformations, accounting for < 8% of all branchial cleft anomalies. However, little is currently known about the cause of FBCAs at the molecular level.
OBJECTIVE: To identify genomic alterations related to the genetic etiology of FBCAs in Chinese children.
METHODS: We performed whole-exome sequencing of samples from 10 pediatric patients with FBCAs. Data analysis was carried out using the Burrow-Wheeler Alignment software package, and the dbSNP database for comparisons. Rare variants were further validated by Sanger sequencing. Insertion/deletions (indels) were examined using the Genome Analysis Toolkit.
RESULTS: We identified 14 non-synonymous mutations in seven potential FBCA-susceptibility genes (TRAPPC12, NRP2, NPNT, SH3RF2, RHPN1, TENM4, and ARMCX4). We also detected 133 shared small indels in 125 genes. Gene Ontology analysis indicated that most of the identified genes played critical roles in development and differentiation pathways involved in regulating organ development.
INTERPRETATION: We characterized the mutational landscape in pathways involved in development and differentiation in Chinese children with FBCA. The results identified potential pathogenic genes and mutations related to FBCA, and provide molecular-level support for the branchial theory of FBCA pathogenesis.

INTRODUCTION

First branchial cleft anomalies (FBCAs) are rare congenital malformations accounting for <8% of all branchial cleft anomalies. FBCAs can occur anywhere from the external auditory canal (EAC) to the mandibular angle. The most frequent clinical symptom of FBCA is recurrent parotid abscesses. The gold standard for the diagnosis of FBCA is preoperative magnetic resonance imaging combined with pathologic diagnosis after surgery. FBCAs are commonly divided into two types according to Work’s classification, based on anatomical and histological features: Type I lesions are of ectodermal origin, while Type II lesions arise from both ectodermal and mesodermal tissues. The recurrent parotid abscesses are insensitive to antibiotics and have a poor outcome following incision and drainage, making surgery the only effective option for treating FBCAs. However, FBCA excision is associated with a high frequency of postoperative EAC stenosis and facioplegia because of the close relationship between the EAC and the facial nerve.

Various theories have been proposed regarding the pathogenesis of FBCA, involving congenital, lymph node, branchial, and precervical sinus origins. In addition, Raguse et al recently suggested high levels of estrogen in pregnancy as a possible pathogenic mechanism for FBCAs. The widely accepted branchial theory is closely related to the development and differentiation of the branchial arches, involving incomplete obliteration of the branchial clefts during embryogenesis. Although some case reports have indicated a possible genetic etiology, the etiology of FBCAs has rarely been examined at the molecular level.

Whole‐exome sequencing (WES) is a next‐generation sequencing (NGS) technique for sequencing all genomic exon regions with high throughput. Although exons only account for about 1%–2% of the human genome, they play the most important role in coding for proteins with subsequent biological functions. Continuing reductions in the cost of WES and the development of NGS instruments have led to the wide application of WES in disease research, especially in the study of cancers and rare diseases. We therefore aimed to identify potential pathogenic genes in 10 Chinese children with FBCAs using WES.

METHODS

Ethical approval

This study was approved by the ethics committee of Beijing Children’s Hospital (No. 2014‐111). All patients’ parents signed informed consent.

Patients and samples

Ten patients with FBCAs treated in Beijing Children’s Hospital (BCH) from 2013 to 2015 were recruited for this study. They were diagnosed by MRI and pathologic analysis. The blood DNA was extracted using the QIAamp DNA Micro kit (Qiagen, Heidelberg, Germany) following the manufacturer’s protocol. Extracted DNA samples were quantified using a Qubit uorometer (Invitrogen, Carlsbad, CA, USA) and were qualified by the 0.8% agarose electrophoresis.

Whole‐exome sequencing

Agilent SureSelect Human All ExonV6 kit (Agilent Technologies, USA) was used for the Illumina Hiseq paired‐end exome sequencing library construction. The adapter‐modified sample was purified using Agencourt AMPure XP beads (Beckman Coulter, Brea, CA, USA) and analyzed on The Agilent 2100 Bioanalyzer system. The Agilent protocol was used to add index tags by post‐hybridization amplification. Finally, all samples were sequenced on an Illumina Hiseq 4000 PE150 system (Illumina, Inc., USA).

WES reads mapping and variants calling

The low‐quality short sequences in the original sequencing data were removed firstly. Paired‐end reads of 150 bp at each end were mapped to the reference genome (UCSC GRCh37/hg19) by Burrows‐Wheeler Aligner (BWA) software. After being sorted and indexed by samtools, PCR duplicate reads were removed using Picard (http://broadinstitute.github.io/picard). Finally, variant calls were performed by Samtools mpileup and bcftools. Meanwhile, Genome Analysis Toolkit was applied for InDel exploration.

Variant annotation and prioritization

Rare variants (minor allele frequency (MAF) < 0.01 in Asian cohort) were filtered by 1000 Genomes Project data (2012 April). The remaining rare variants were annotated by ANNOVAR software. Variants excepted synonymous, noncoding or non‐frameshift ones were submitted to SIFT and PolyPhen‐2 (http://genetics.bwh.harvard.edu/pph2/) for functional prediction.

Verification by Sanger sequencing

The in predispose genes’ variants identified through whole‐exome sequencing were verified by PCR and Sanger sequencing. Specific PCR primers were designed by Primer 5 software (Table S1).

Three‐dimensional protein modeling

Three‐dimensional structure of NRP2 (PDB No. 4QDQ) was obtained from NCBI (https://www.ncbi.nlm.nih.gov/). According to the domains’ information about NRP2 (O60462) from UniProt (https://www.uniprot.org/), visualization of NRP2 structure and mutants were performed by PyMOL (www.pymol.org).

RESULTS

Clinical characteristics of patients with FBCA

The patient characteristics are summarized in Table 1. There were ten patients (four male, and six female; median age 67 months, interquartile range 25–90 months), with no history of genetic or infectious diseases. Seven patients had left‐sided lesions, and the remainder were located on the right. All 10 patients had undergone incision and drainage. All four type II lesions were close to the facial nerve and adhered to the inferior wall of the EAC cartilage, potentially increasing the risk of postoperative facial palsy and EAC stenosis. However, the facial nerve was unaffected in all patients in the current study, with no recurrences to date.

TABLE 1

Features of 10 patients with FBCA

Characteristic	Number of patient
Sex
Male	6
Female	4
Incision History
Yes	10
No	0
Work’s Classification
I	6
II	4
Close to the facial nerve
Yes	4
No	6
Relationship with EAC
Adhere with the inferior wall of EAC cartilage	4
Adhere with the posterior wall of EAC cartilage	6

FBCA, first branchial cleft anomaly; EAC, external auditory canal.

WES and coverage

About 10 billion bases of effective sequence data were yielded using the Illumina HiSeq 4000 (Illumina, Inc., USA). The average read length was 120 bases. After mapping to the human reference genome (UCSC GRCh37/hg19) using the Bowtie 2 alignment tool, the average sequencing depth on target and the coverage for the target regions of 10× for each sample were obtained (Table S2).

Mutation detection

Fourteen rare mutations, including five reported variants related to seven genes, were discovered by WES of the 10 cases (Table 2). All the variants were verified by Sanger sequencing, and the results were in accordance with those of WES (Figure 1). Fourteen variant sites were verified: TRAPPC12 (c.346G>T, c.476T>C), NRP2 (c.1412G>A, c.1519C>T), NPNT (c.289C>T, c.1606C>T), SH3RF2 (c.37C>G, c.2092G>A), RHPN1 (c.295C>T, c.1093G>A), TENM4 (c.1132G>T, c.6666C>G), and ARMCX4 (c.1121C>G, c.3235G>A). ARMCX4 (c.1121C>G) was a homozygous mutant, and all the others were heterozygous mutants (Figure 1). According to the American College of Medical Genetics Standards and Guidelines, NRP2 (c .1412G>A, c .1519C>T), NPNT (c .289C>T), RHPN1 (c.1093G>A), and TENM4 (c.6666C>G) were considered as like lypathogenic variants, and the remainder were variants of uncertain significance. Detailed information on the affected amino acid sites is presented in Table 2. We also detected 133 shared small indels in 125 genes (Table S3).

TABLE 2

Rare variants identified by whole‐exome sequencing and possible effects in patients with FBCAs

Symbol	Chr	Position	ID	Nucleic acid	Amino acid	Domain	SIFT	Polyphen2_HVAR	Polyphen2_HDIV	Patient ID
TRAPPC12	2	3391740	–	G > T	G > C	–	Deleterious (0.04)	Probably damaging	Probably damaging	2
		3391870	–	T > C	V > A	–	Deleterious (0.01)	Possibly damaging	Probably damaging	3
NRP2	2	206608047	rs149849497	G > A	R > H	F5/8 type C 2	Deleterious (0)	Probably damaging	Probably damaging	4
		206608154	–	C > T	R > C	F5/8 type C 2	Deleterious (0)	Probably damaging	Probably damaging	3
NPNT	4	106858189	–	C > T	P > S	EGF‐like 2;	Deleterious (0.04)	Probably damaging	Probably damaging	7
		106890142	rs147786422	C > T	R > C	–	Deleterious (0)	Benign	Probably damaging	5
SH3RF2	5	145317528	–	C > G	P > A	–	Deleterious (0.01)	Probably damaging	Probably damaging	4
		145442166	rs199894169	G > A	G > R	–	Deleterious (0.03)	Probably damaging	Probably damaging	3
RHPN1	8	144458813	–	C > T	R > W	–	Deleterious (0)	Possibly damaging	Probably damaging	10
		144462146	rs576390102	G > A	D > N	BRO1	Deleterious (0.04)	Possibly damaging	Probably damaging	6
TENM4	11	78574130	rs1025280520	G > T	Q > K	–	–	Probably damaging	Probably damaging	8
		78380724	–	C > G	N > K	YD Repeat	–	Probably damaging	Probably damaging	3
ARMCX4	X	100744697	–	C > G	S > C	–	–	–	–	2
		100746811	–	G > A	G > R	–	–	–	–	1

FBCAs, first branchial cleft anomalies; Chr, chromosome; SIFT, sorting intolerant from tolerance; – , not available.

FIGURE 1

Confirmatory analysis of FBCA‐related genes detected via whole‐exome sequencing using Sanger sequencing. (A) (B) TRAPPC12, (C) (D) NRP2, (E) (F) NPNT, (G) (H) SH3RF2, (I) (J) RHPN1, (K) (L) TENM4 and (M) (N) ARMCX4.

Variant analysis

Alignment of multiple protein sequences of TRAPPC12, NRP2, NPNT, SH3RF2, RHPN1, TENM4, and ARMCX4 revealed that 11 of 14 amino acid sites (p.V159A of TRAPPC12, p.R471H and p.R507C of NRP2, p.R565C of NPNT, p.P13A and p.G698R of SH3RF2, p.R99W of RHPN1, p.Q378K and p.N2222K of TENM4, and p.S374C and p.G1079R of ARMCX4) were located in amino acid sequences that were highly evolutionarily conserved among different species, including Homo sapiens, Gorilla, Pan troglodytes, Pongo abelii, Canis lupus familiaris, Oryctolagus cuniculus, and Bos taurus (Figure S1). The other three mutations were not associated with highly conserved protein sequences (p.G116C of TRAPPC12, p.P97S of NPNT, and p.D365N of RHPN1).

Protein structure analysis

Out of the seven genes, only the crystal structure of NRP2 has been analyzed in humans ((PDB No. 4QDQ). Three‐dimensional protein modeling showed that R471 and R507 of NRP2 were located in the F5/8 type C2 domain (PROSITE‐ProRule: PRU00081), right next to the F5/8 type C1 domains, indicating that these mutations might cause functional changes in the NRP2 protein (Figure S2).

Gene ontology (GO) analysis

We used the DAVID Bioinformatics Resources 6.8 (NIAID/NIH) and PANTHER database to analyze GO terms for the seven genes (Figure S3). There were no significantly enriched GO terms according to DAVID, but NRP2, NPNT, and TENM4 all play key roles in differentiation and development. The PANTHER database indicated that NRP2, RHPN1, and TENM4 were all involved in developmental process (GO: 0032502). Four of the seven genes were therefore closely related to cell differentiation and development, thus supporting the branchial theory.

Regarding the 133 shared small indels from 125 genes, the PANTHER database showed that 13 different GO terms were enriched (Figure S4). Most genes were involved in developmental process (GO: 0032502), cellular process (GO: 0009987), cell proliferation (GO: 0008283), biological regulation (GO: 0065007), metabolic process (GO: 0008152), and metabolic process (GO: 0008152). These results suggest that cell differentiation and development may be disturbed at the embryonic stage in FBCA.

DISCUSSION

As a rare congenital malformation of the head and neck, FBCAs account for < 8% of all branchial cleft anomalies. Surgical repair is the only effective treatment for FBCAs; however, relapse may occur because the surgical procedures are relatively complex and because of wide heterogeneity among patients. More information is therefore needed regarding the pathogenesis of FBCAs. We accordingly examined 10 children with FBCA using WES, and revealed 14 differential rare mutations located in seven potential susceptibility genes. Database screening identified five out of the 14 variants. Variant analysis showed that 11 out of the 14 mutations were located in regions that were highly evolutionarily conserved among different species, suggesting that these seven genes may play key roles in the etiology of FBCAs.

The current clinical observations and previous study indicated that FBCAs may be caused by the abnormal development and differentiation of the branchial arches. The current GO analysis accordingly revealed that NRP2, RHPN1, TENM4, and NPNT were involved in developmental processes.

NRP2 is a member of the neuropilin family of receptor proteins with essential roles in endothelial cell migration, which in turn contributes to angiogenesis and vascular development. R471 is a necessary binding site on NRP2 for vascular endothelial growth factor (VEGF) proteins, which modulate the development of organs derived from all three germ layers (endoderm, mesoderm, and ectoderm). The R471H mutation may thus influence the interaction between NRP2 and VEGF, and thus disturb further functions of these two proteins. R507 of NRP2 is not a binding site, but R507C mutation may change the side chain of amino acids. Changes in the spatial structure of NRP2 may then disrupt the interaction between NRP2 and VEGF, further influencing organ development, and potentially inducing FBCAs.

RHPN1 is a Rho GTPase‐interacting protein that interacts with cytoskeletal components upon Rho binding. Rho GTPases are major determinants of actin cytoskeletal dynamics and thus control a wide variety of morphogenetic events. D365 of RHPN1 is located in the BRO1 domain, which binds with multivesicular body components such as yeast Snf7 and mammalian CHMP4b, and can function to target BRO1 domain‐containing proteins to endosomes. Although the biologic function of RHPN1 is largely unknown, recent findings suggest that it acts as a key regulator of the podocyte cytoskeleton and is necessary for glomerular filtration.

TENM4 is a transmembrane protein expressed in the nervous systems and mesenchymal tissues, including cartilage. TENM4 is required for mesoderm induction in mice and has been found to suppress chondrogenic differentiation in humans. TENM4 can also regulate differentiation of muscle satellite cells during muscle regeneration. N2222 is located in the YD repeat in TENM4, which appears to be involved in binding carbohydrates.

NPNT encodes an extracellular matrix protein and is considered to play critical roles in the development and functions of various organs and tissues, including the kidneys, bones, teeth, and atrioventricular canal. Moreover, the mutation of P97S of NPNT belongs to the epidermal growth factor (EGF)‐like calcium‐binding domain, which includes six cysteine residues shown to be involved in disulfide bonds. A calcium‐binding site has also been found at the N terminus of some EGF‐like domains, and calcium‐binding may be crucial to numerous protein–protein interactions.

Three other identified gene products, TRAPPC12, SH3RF2, and ARMCX4, also play critical roles in the development and differentiation pathways regulating organ development. TRAPPC12 is involved in protein transport from the endoplasmic reticulum to the Golgi apparatus. TRAPPC12 mutants have been associate d with the induction of progressive childhood encephalopathy, which involves progressive central nervous system dysfunction. SH3RF2 is an E3 ubiquitin‐protein ligase (SH3 domain‐containing ring) that acts as an anti‐apoptotic regulator of the JNK pathway and functions as an oncogene. The functions of ARMCX4 have been largely unexplored, except for a study by Chang et al, which revealed that ARMCX4 may be a potential passenger gene inducing endometrial cancer. Notably, ARMCX4 (c.1121c>G) was the only homozygous mutant in the current study, suggesting that a change in this locus could severely affect the function of ARMCX4. However the specific mechanisms require further study.

The current study also detected 133 shared small indels in 125 genes, enriched for 13 different GO terms, mostly related to cell differentiation and development. Although further studies are needed to determine the mechanisms of these genes in FBCA, these shared small indels may disturb cell differentiation and development at the embryonic stage, leading to FBCAs.

This study had some limitations. First, we could not confirm if the potential pathogenic mutations arose de novo or were inherited, because of the lack of parents’ blood samples. Second, the small sample size meant that it was not possible to validate the pathogenic mutations in a different cohort. In addition, although the 14 potential pathogenic mutations and their seven related genes were detected by WES and confirmed by Sanger sequencing, and their biological functions were analyzed by bioinformatics methods, further studies of the functions and mechanisms of these genes and their variants are required to clarify the pathogenesis of FBCA. We aim to continue to follow‐up the current patients and expand the sample size, and conduct prospective studies to verify these pathogenic mutations.

In summary, we identified 14 rare variants in seven genes (TRAPPC12, NRP2, NPNT, SH3RF2, RHPN1, TENM, and ARMCX4) and 133 common small indels from 125 genes by WES screening of 10 children with FBCAs. All these genes were involved in development and differentiation pathways. This study therefore identified potential pathogenic genes of FABC, further supporting the branchial theory of pathogenesis at the molecular level.

CONFLICT OF INTEREST

No potential conflict of interest relevant to this article was reported.

Supplementary Material

Click here for additional data file.