Multi-omics analysis of a case of congenital microtia reveals aldob and oxidative stress associated with microtia etiology

Abstract

Background

Microtia is reported to be one of the most common congenital craniofacial malformations. Due to the complex etiology and the ethical barrier of embryonic study, the precise mechanisms of microtia remain unclear. Here we report a rare case of microtia with costal chondrodysplasia based on bioinformatics analysis and further verifications on other sporadic microtia patients.

Results

One hundred fourteen deleterious insert and deletion (InDel) and 646 deleterious SNPs were screened out by WES, candidate genes were ranked in descending order according to their relative impact with microtia. Label-free proteomic analysis showed that proteins significantly different between the groups were related with oxidative stress and energy metabolism. By real-time PCR and immunohistochemistry, we further verified the candidate genes between other sporadic microtia and normal ear chondrocytes, which showed threonine aspartase, cadherin-13, aldolase B and adiponectin were significantly upregulated in mRNA levels but were significantly lower in protein levels. ROS detection and mitochondrial membrane potential (∆ Ψ m) detection proved that oxidative stress exists in microtia chondrocytes.

Conclusions

Our results not only spot new candidate genes by WES and label-free proteomics, but also speculate for the first time that metabolism and oxidative stress may disturb cartilage development and this might become therapeutic targets and potential biomarkers with clinical usefulness in the future.

Backgrounds

Microtia-anotia (MIM: 600674) is a congenital malformation of the external and middle ear caused by the abnormal development of the first and second pharyngeal arch and the first sulcus at the embryonic stage. It may be accompanied by hearing loss, and also mandibular and facial soft tissue dysplasia [1]. Estimates of the incidence of microtia vary widely and range from 0.83 to 17.4 per 10,000 births [2]. Existing data indicate that inheritance is more likely in syndromic or familial microtia whereas multifactorial or polygenic causes are more probable in sporadic cases [3]. Since the fact that the ear is not only a functional organ, but also a crucial part of appearance, therefore, patients with microtia and their caretakers have suffered both physical and psychosocial burdens. Hence, to discovering the etiology of microtia has been a major topic in the plastic surgery field and public health.

About 40% of microtia patients present with other structural abnormalities as a syndrome [4], which involves chromosomal abnormalities and single-gene disorders. As goosecoid gene (GSC) has been identified as the most interesting candidate gene within the 14q32 chromosomal region, which was found to be connected with hemifacial microsomia and Goldenhar syndrome [5]. On the other hand, genetic variations that disrupt neural crest cells (NCCs) development can also be considered as an aetiology factor for sporadic cases of microtia. It is widely believed that HOXA1 and HOXA2 deficiency affecting the branchial arch development can cause microtia and other associated deformities. Meanwhile, some important pathways are also likely to affect the occurrence of microtia, such as bone morphogenetic proteins (BMPs), Wingless/INT (WNTs), fibroblast growth factors (FGFs) [6], and retinoic acid et al. [3]. Recently, epigenetic modifications such as DNA methylation and noncoding RNAs have been verified to affect the occurrence of microtia [7]. However, up till now, no mechanism mentioned oxidative stress and energy metabolism as the microenvironment’s imperative role in the pathogenesis of congenital microtia.

Here in this report, we study a rare case of congenital microtia with costal-chondrodysplasia based on the whole-exome sequencing, several candidate genes and related pathways are identified and new CNVs are additionally spotted. Secondly, label-free proteomics is used to assess microtia-associated protein changes. Altogether, we find that genes and the significant pathways are significantly related to the oxidative stress and energy metabolism process. Then real-time PCR and IHC show that ALDOB, ADIPOQ, CDH13 and TASP1 are probable pathogenic genes and ROS analysis and JC-1 assay both prove that oxidative stress exists in microtia. The expectation is that more insights into the molecular circuitry of microtia could be gained and more potential biomarkers and amelioration of microenvironment during embryonic development should benefit disease treatment and prevention.

Results

Clinical Report

Here we present a sporadical microtia case of non-consanguineous parents. The male patient A, who was initially assessed at 20 years of age presenting a weight of 93 kg, height of 166 cm, head circumference of 59 cm, and showed left ear with third degree microtia [8]. Mother (C), father (D), and the sister (B) show no abnormality (Fig. 1D-F). None of the family members show psychomotor retardation. A is diagnosed as left-side congenital microtia. During our three-dimensional (3D) chest computed tomography (CT) (Fig. 1A, B), the patient has abnormal rib cartilage which led to the failure of the reconstruction of the outer ear because the lack of a well-sculptured cartilage framework. Meanwhile, the patient does not show abnormality in the vertebral column (Fig. 1 C). The patient showed no gestational exposure to specific medications according to the mother’s recall, so they did not perform prenatal diagnosis by karyotyping or chromosomal microarray analysis.

Fig. 1

The phenotype of the patients and family members A Three-dimensional (3D) Computed Tomography (CT) reconstructions of the patient head showed left-side microtia B 3D-CT reconstruction of the patient costal area showed costal chondrodysplasia with sternal dysplasia. The costal cartilage is missing on the 6–10 and not attached to the rib. C X-ray thorax shows no deformity of the vertebral column (D-F) 3D-CT results of the patient's father (D), mother E and sister’s F normal costal area

WES

ACMG Classification of single nucleotide polymorphisms (SNPs)

One hundred fourteen deleterious InDel and 646 deleterious SNPs are screened out (Table 1 & Supplementary Table 1 and 2). SNPs are also classified by American Society for Medical Genetics and Genomics (ACMG) system [9] (Supplementary Table 3), among which 2 are pathogenic: Potassium voltage-gated channel, KQT-like subfamily, member 4 (KCNQ4) and Lipopolysaccharide-responsive, beige-like anchor protein (LRBA). 8 are likely pathogenic: Leptin receptor (LEPR), G protein-coupled receptor 161 (GPR161), Superoxide dismutase-3 (SOD3), WD repeat-containing protein 36 (WDR36), Plasminogen (PLG), Protein kinase, DNA-activated, catalytic polypeptide (PRKDC), unc-45 myosin chaperone B (UNC45B), Threonine aspartase (TASP1).

Table 1 Variants in 4 samples identified by WES

These variants include 1 frameshift variant: LEPR (NM_001198687: exon17: c.2534delT:p.I845fs) and 6 missense variants: KCNQ4 (NM_004700: exon4:c.C546G:p.F182L),GPR161(NM_001267613:exon2:.C802G:p.R268G),SOD3(NM_003102:exon2:c.C691G:p.R231G),WDR36(NM_139281:exon9:c.G1100A:p.G367D), PLG (NM_000301:exon15:c.G1858A:p.A620T), UNC45B (NM_001308281: exon15:c.G2071A:p.E691K). 1 stop gain: LRBA (NM_001199282: exon35:c.G5586A:p.W1862X) and 2 splice variants: PRKDC and TASP1 (NM_001323603: exon7:c.98–1 > TA) (supplementary Table 3). Then, the following candidate genes ranked in descending order according to their relative impact on the development or severity of microtia: PRKDC > LEPR > PLG > MUC6 > DGKK > GPR161 > WDR36 > BPTF > LRBA > SOD3 > KCNQ4 (Fig. 2A).

Fig. 2

Candidate genes classified by the American Society for Medical Genetics and Genomics (ACMG) system and GO enrichments A Candidate genes are ranked according to how strongly they were associated with the disease by the ACMG system; The maximum correlation score is 1. B GO enrichment distribution map by molecular function (MF); C biological process (BP); D cell component (CC). The abscissa represents the proportion of enriched genes in the pathway to the total enriched genes, and the ordinate represents the name of the enriched GO term. The size of the dot represents the number of genes enriched in the term of the gene, and the color represents the p-value

Novel Variant or CNV Identification

The patient lacks a family history of genetic disease, that is, the pathogenic variants in the patient may be de novo mutations. Combining exome analysis, we identified de novo CNVs in espin pseudogene (ESPNP), which is located in chromosomal 1p36.13, starting from 17,017,581 and ending at 17,030,650, the size is 13,070 bp, the CNV type is duplication and function as ncRNA exonic. Meanwhile, CNVs defined as possibly deleterious locate on chromosomal X, start from 1,571,538 to 2,139,310, with the size of 567,772, are acetyl serotonin methyltransferase (ASMT), dehydrogenase/reductase X-linked (DHRSX), acetyl serotonin O-methyltransferase like (ASMTL), A-kinase anchoring protein 17A (AKAP17A), P2Y receptor family member 8 (P2RY8). The CNV type is duplication (Supplementary Table 4). None of them have been reported before in microtia. Worth notice is that, the patient shares possibly deleterious CNVs on ASMT, AKAP17A and DHRSX with his father but not his mother (Supplementary Table 4) while his parents show no abnormality and this may be due to the incomplete penetrance.

GO Analysis of the Mutant Genes

The GO analysis showed that the metabolism-, embryonic development-, communication-, and membrane-related genes were significantly enriched in the microtia. 11 pathways were significantly involved in molecular functions (MF) (Table 2), they are mainly about lipid metabolism and endopeptidase activity, some may involve copper ion binding and G-protein coupled amine receptor activity, while extracellular matrix structural constituent may also play a role in the formation of cartilage.

Table 2 GO enrichment of WES DEGs

One hundred thirteen pathways of biological process (BP) were potentially related to the pathogenesis of congenital microtia since most are associated with embryo development related process (Table 2), top 10 among which are single-multicellular organism process (10/13), multicellular organismal process (10/13), multicellular organismal development (8/13), single-organism developmental process (8/13), developmental process (8/13), system development (7/13), embryo development (5/13), anatomical structure homeostasis (4/13), chordate embryonic development (4/13).

Twenty-two pathways of cell component (CC) reflect that the endocrine system may play crucial roles in auricle development (Table 2), among which are plasma membrane (7/13), cell periphery (7/13), organelle lumen (6/13), cytoplasmic membrane-bounded vesicle (3/13), cytoplasmic vesicle (3/13) etc. (Table 2).

Label-free proteomics

Identification of proteins involved in microtia

Forty-nine were significantly differentially expressed between the patient and sister (Table 3). Among those are 28 upregulated. Haptoglobin (HP) is increased by 10 times in group of A vs. B while in group of A vs. D the fold change is more than 6(Table 5). Protein S100-A9 (S100A9) in this comparison is 10 times higher, and repeat in group A vs. D by 4 times higher. Protein S100-A8 (S100A8) in the firs A vs. B and the A vs. D group (Table 5) is also upregulated with about 4 times and 2 times respectively. 21 are downregulated including Cadherin-13 (CDH13) and Adiponectin (ADIPOQ), Sex hormone-binding globulin (SHBG), Coagulation factor XIII A chain (F13A1), Cell surface glycoprotein MUC18 (MCAM), this trend remains in the group A vs. D and the latter two remain the same in group A vs. C (Table 4).

Table 3 Significantly differently expressed proteins (DEPs) between the patient(A) and the sister (B) / father (C) /mother (D)

Table 4 Top 5 over-represented GO terms mapped to each category of enriched proteins (A vs B)

Among the 23 proteins that were significantly different between patient and mother (Table 4), 18 are upregulated. Immunoglobulin delta heavy chain (IGHD) repeat in three groups, but the highest fold change is in this group by more than 7 times. Apolipoprotein A-II (APOA2) increases in A vs. B by 6 times and in A vs. C by 2 times. Pantetheinase (VNN1) in the A vs. B and A vs. C is both about 4 times higher. Peroxiredoxin-2 (PRDX2) is upregulated in the first two groups by 2 times. Five proteins are downregulated, including EGF-containing fibulin-like extracellular matrix protein 1 (EFEMP1) repeats in the first two groups.

There are 27 significantly different proteins between the patient and the father (D), 20 are upregulated including and Serum amyloid A-1 protein (SAA1) by 21-times higher, Alpha-1-antitrypsin (SERPINA1) repeat in A vs. C and A vs D, the fold changes are 13 and 11 respectively. The same trend goes with Serotransferrin (TF), in A vs. C the fold change is about 2 while in A vs. D fold change is more than 10 times. Immunoglobulin alpha-2 heavy chain (IGA2) in the first and this group is increased by the same 2-times (Table 5)0.7 proteins are downregulated in this group.

Table 5 Top 5 over-represented GO terms mapped to each category of enriched proteins (A vs C)

Ten proteins that overlap between three groups are Alpha-1-acid glycoprotein 2 (ORM2), SAA1, IGHD, Polymeric immunoglobulin receptor (PIGR), Hemoglobin subunit beta (HBB), Fructose-bisphosphate aldolase B (ALDOB), SOD3, Plasminogen (PLG), F13A1, MCAM. Especially for ORM2 is even for more than 52 times up-expressed in A vs. B and 22 times higher in A vs. D group. SOD3 is halved in the comparison between the patient and the sister, while upregulated more than 10 times in the group between the patient and the mother but remains unchanged in comparison of patient to father. PLG appears the same in three groups, of which fold-change is about 1. F13A1 and MCAM are downregulated in three groups.

When compare the proteomics results to the WES SNPs genes, there are overlapped results: Mannan-binding lectin serine protease 2, MASP2 (chromosome 1 11,105,527 C > T exonic, missense SNVs) and Ficolin (collagen/fibrinogen domain-containing lectin) 2, (FCN2) (NM_004108: exon5: c.429 + 1G > A; NM_015837: exon4: c.315 + 1G > A).

GO enrichment analysis

The top 5 over-represented GO terms, including MF, BP, and CC categories of the common up and downregulated proteins were summarized in Tables 6, 7 and 8. GO level 2 analysis shows that proteins in three groups in the BP were principally enriched by biological process, biological regulation, single-organism process, and cellular process. In the MF, proteins were mainly enriched in molecular function, protein binding, ion binding, receptor binding and so on. In the CC, proteins were mainly enriched in cellular component, cell part, organelle, extracellular region part, and extracellular region.

Table 6 Top 5 over-represented GO terms mapped to each category of enriched proteins (A vs D)

Table 7 Proteins GO terms Enrichment level 2

Table 8 Top 5 over-represented GO terms mapped to each category of enriched proteins (A vs D)

When comparing the proteins of the patient to the database at level 2, GO enrichment showed that the proteins were mainly enriched in lipid metabolism-related GO terms and oxidative stress related areas (Table 9 and Fig. 3).

Table 9 Proteins GO terms Enrichment level 2

Fig. 3

GO annotation of proteins level 2 Analysis A. Significantly different proteins when the patient A compared with sister B assigned by GO annotation; B patient A compared with mother C; C patient A compared with father D. Dark red means the biological process proteins contribute to, blue columns represent proteins’ molecular function, and orange means in which cellular component the proteins are active

KEGG pathway analysis

To gain an initial understanding of the role and function of the identified protein differences, we used KEGG pathway analysis to identify the biological pathways of the proteins that were significantly differentially expressed (> 1.5 fold upregulated or downregulated) between microtia and control samples. In Fig. 4, the KEGG classification results showed that genes were mainly enriched in PPAR signaling pathway, Cholesterol metabolism, Fox O signaling pathway, HIF-1 signaling pathway, IL-17 signaling pathway and so on (Fig. 4). In the group of A Vs. C (mother), the top 30 signaling pathway results showed that proteins are mainly enriched in Fox O signaling pathway, PPAR signaling pathway, HIF-1 signaling pathway, Amyotrophic lateral sclerosis (ALS), EGFR tyrosine kinase inhibitor resistance etc. (Fig. 4). When compared to his father, the IL-17 pathway is the only pathway that is significantly affected.

Fig. 4

KEGG pathways enrichments of proteins. A The significantly different proteins in group patients A vs. sister B enrich in PPAR signaling and cholesterol metabolism pathway based on the KEGG database; B. The significantly different proteins in group A vs. father C enrich in many pathways, including FoxO, MAPK, PPAR, etc. C. The significantly different proteins in group A vs. mother D enrich in IL-17 signaling pathway

PPI network analysis.

Proteins with similar expression patterns may have similar functions or participate in the same biological pathways, or be in adjacent regulatory positions in the pathways. To this end, to categorize the data based on functional similarity, we show the PPI network of differently expressed proteins. Some differently expressed proteins between the patient and the sister (Fig. 5 A V B) are correlated with each other, S100A8 is highly related to S100A9, APOA2, APOC1-3, APOE, and APOD are in the center of the STRING while SAA1, SAA4, HP, are closely related to the proteins mentioned above. SOD3 and PRDX2 repeat in the A vs. C.

Fig. 5

Protein–protein interaction (PPI) network. Network representation of signed PPIs shows significantly different proteins complexes and signaling pathways. Circles represent proteins while lines indicate molecular interactions

Further verifications on Other Microtia Patients

As the variations of proteins between the patient and the controls don't necessary cause anomaly, further evaluation is necessary to investigate the potential roles of candidate genes identified in the multi-omic analysis. To accomplish this, we collected auricular cartilage tissues from individuals with microtia as well as from individuals who underwent rhinoplastic surgery or experienced accidental ear injuries. Firstly, we screened the differentially expressed genes. Real-time PCR results in Fig. 6 showed that TASP1, CDH13, ALDOB and ADIPOQ were significantly upregulated than the normal control group while the absolute expression level of ADIPOQ was not as high as other three, and the variance of ALDOB among different patients is higher. Secondly, we performed IHC of the four genes aforementioned on microtia patients and normal people as the control to verify the expression on protein level. The results showed that ALDOB, TASP1, and CDH13 were all significantly lower in the microtia group, while ADIPOQ showed no significance (Fig. 6).

Fig. 6

Candidate genes verifications between microtia patients and normal people. A Real-time PCR results mRNA level comparison. TASP1, CDH13, ALDOB and ADIPOQ are upregulated significantly in microtia patients’ ear chondrocytes when compared to the normal people as control. B IHC images and protein expression level comparison IHCs involved in the revealed markers of TASP1, CDH13, ALDOB and ADIPOQ are shown for microtia and normal ear cartilage. TASP1, CDH13 and ADIPOQ levels are significantly lower in the microtia group. Scale bars:10X: 200 μm, 20X: 100 μm.*p < 0.05, **p < 0.01. Error bars represent the mean ± SD of three independent experiments

ROS Detection and mitochondrial membrane potential (∆Ψm) detection proved that oxidative stress exists in microtia chondrocytes

As we mentioned above, both WES and proteomic results indicate that cellular redox homeostasis may play a role in microtia, therefore, we test this hypothesis by ROS experiment and JC-1 experiment. Chondrocytes from different microtia patients were used as the experimental group, normal people’s ear chondrocytes were the control group and the positive control group was stimulated with 1ug/ well Rosup for 20 minutes. After detecting fluorescence intensity, ROS levels in the microtia group were significantly higher than those in the control group (p < 0.01) (Fig. 7A).

Fig. 7

ROS Detection and Mitochondrial Membrane Potential (∆Ψm) Detection in Microtia Chondrocytes. A ROS levels of chondrocytes in patients with microtia in control cells. B: mitochondrial membrane potential (MMP/ΔΨm) based on JC-1 staining. C MMP fluoresence intensity. MMP is significantly lower in microtia chondrocytes compared to the ctrl group while maintaining a balance in healthy chondrocytes. D The monomers/aggregates fluorescence ratio of normal and microtia chondrocytes. The mitochondrial membrane potential in microtia chondrocytes was significantly lower than that of normal ear chondrocytes. Error bars represent the mean ± SD of three independent experiments. *p < 0.05, **p < 0.01

As mitochondria are the predominant source of ROS production, we measured mitochondrial membrane potential (MMP/ΔΨm) based on JC-1 staining (Fig. 7B). Red fluorescence represented a potential-dependent aggregation in the mitochondria, while green fluorescence represented the monomeric form of JC-1, appearing in the cytosol after mitochondrial membrane depolarization. We found that MMP was significantly lower in microtia chondrocytes compared to the ctrl group, moreover, ΔΨm maintained a balance in healthy chondrocytes (Fig. 7C). Our results showed that the green/red fluorescence ratio increased in the microtia group. In normal cells, JC-1 aggregated in mitochondria and the ratio was 0.91 ± 0.16. Microtia chondrocytes showed significantly a higher ratio (420.30 ± 154.60, P < 0.01) because the monomeric form JC-1 appeared in the cytosol, which indicated the dissipation of ΔΨm (Fig. 7D). These results collectively support the hypothesis that oxidative stress exists in microtia patients and may affect the mitochondrial function in chondrocytes.

Discussion

Microtia is one of the leading congenital malformations in the plastic surgery field, it brings harm to children's psychology, and puts pressure on their families and increases society's burden, experts are attaching increasing importance to the exploration of microtia. One major obstacle in treating microtia is that the underlying mechanism is complicated and not clearly understood. Owing to the fact that auricular cartilage originates from NCCs [10], therefore, previous reports about microtia mainly focused on the development of NCCs. In this study, we further investigated the etiology of microtia accompanied by costal chondrodysplasia and sternal bone deformity by bioinformatic tools. WES and label-free analysis were conducted between the family members. Interestingly enough, both WES and proteomic results indicated that metabolic program, associated immune responses, and inflammatory and cellular redox homeostasis may play a pathological or at least a modulating role in microtia, and this has never been reported before.

The candidate genes SOD3 and PLG by ACMG classification overlap with the proteomic results. PLG not only regulates coagulation and complement pathway, it is also responsible for neutrophil apoptosis and efferocytosis in inflammation, and functions as apolipoprotein binding [11], whereas in our case, it is downregulated in microtia patients when compared to his mother but upregulated than his father. It is widely believed that excessive reactive oxygen species (ROS) can lead to oxidative stress and DNA damage which can decrease MSCs self-renewing and multidirectional differentiation [12]. SOD3 protects cells from the toxic effect of RO intermediates by converting superoxide radicals into hydrogen peroxide and oxygen. Nightingale et al. found that decreased SOD3 level following BMSCs chondrogenesis [13]. Gavriilidis et al. knocked down SOD2 in chondrocytes, resulting in increased ROS levels, mitochondrial DNA strand breaks and decreased antioxidant capacity, and this may be a potential contributor to osteoarthritis [14].

The variant identified in the PRKDC gene is predicted to be pathogenic according to ACMG, it can recognize and repair DNA double-strand breaks [15], and interacts with autoimmune regulator (AIRE) to regulate B cells and natural killer (NK) cells, causing inflammation and immune dysregulation, and downstream S100A8 expression [16]. Plus, PRKDC is also involved in nervous system development [17]. Although S100A8 does not show significance in WES, it is upregulated in proteomics results. S100A8 is a calcium- and zinc-binding protein that plays a prominent role in the regulation of inflammatory processes and immune response. S100A8/A9 promotes cell death via autophagy and apoptosis and this occurs through the crosstalk of mitochondria and lysosomes via ROS [18]. The S100 family is closely related to chondrogenesis [19], and S100A8 and S100A9 up-regulate inflammatory cytokines through TLR-4 And MMP-1, -3, -9, and -13 to destroy bone [20].

TASP1 encodes taspase 1(threonine aspartase 1, TASP1). TASP1 proenzyme intramolecularly proteolyzes and generates an active N-terminal 28 kDa α subunit and a 22 kDa C-terminal β subunit heterodimer whereas the β subunit is the active site for cleavage activity [21]. TASP1 cleaves various nuclear factors after an aspartate with substrates including the histone-methyltransferase mixed lineage leukemia 1 (MLL1/KMT2B), transcription factor IIA (TFIIAα-β/GTF2A1), TFIIA-like factor (ALFα-β) and so on [22]. TASP1 activateing MML1 is crucial for HOX and cyclin genes expression that participate in body segmentation and cell proliferation [23] and TASP1 cleaving TFIIAα-β make it more susceptible to degradation impedes its transcriptional activity in embryonal cell proliferation and morphogenesis. Tasp1-/- mice show craniofacial malformations [24] and facial features containing low-set ears [25] and TASP1 mutation has been reported in complex syndromes manifesting with facial and skeletal abnormalities [26]. As mitochondrial redox status can affect SAM level to modify histone methylation [27], and studies have identified SET/MLL as redox sensitive units of the H3K4me3 [28]. The relationship between Taspases and oxidative stress-mediated alterations in DNA/histone methylation and other post-translational modifications is worth to be explored in future studies.

The enrichment analysis of differently expressed genes (DEGs) revealed that mutated genes are majorly enriched in the cellular metabolic process, embryonic development process, oxidative stress and cytoplasmic vesicles, which are concordant with the results in proteomics analysis. It may reflect that metabolism may play an important regulating role in the pathogenesis of microtia. Expression differences could be a sign of the biochemical or physiological changes that might lead to microtia. Many researches have proved that immunometabolism is closely related to the development of diseases [29]. Oxidized low-density lipoprotein can trigger strong proinflammatory responses that can potentially contribute to the development of atherosclerosis [30].

What gives us implications in our study is the upregulated proteins in three groups are closely related to the immune system and oxidative stress. The highest change of protein is ORM2, which is reported to modulate the activity of the immune response and the complement cascades pathway [31]. ORM2 is at the same time a suppressor of Fe-related ROS [32]. SAA1 is also a major acute-phase inflammatory protein [33], it contributes to bone and cartilage destruction [34]. IGHD as its name implies, mediate the effector phase of humoral immunity [35]. Polymeric immunoglobulin receptor (PIGR) specifically binds to J-chain-containing secretory IgM and IgA, therefore, it is a critical protein in the mucosal immunity [36]. Vertebrate haemoglobin, a heterotetramer of paralogous α- and β-subunits that mediate respiratory oxygen transport and exchange [37], and can be cleared by binding to the haptoglobin (HP), in this way to scavenging HBB-mediated toxic effects caused by free radical and participates in the complementary pathways [38].

When comparing the proteomics results to the WES SNPs genes, there are overlapped results: MASP2 and FCN2, which are all belong to key proteins of the lectin pathway of complement. Ficolin 2 is the initiator and their N-terminal collagen-like domain interacts and forms complexes with the MASPs, which leads to the C4 and C2 cleavage and thus activation of the downstream complement cascade [39]. In our results, FCN2 is at least fivefold higher, and the variation is heterogenous in his mother. MASP2 is also at least fivefold higher and the variation is heterogenous in the patient’s sister and father, both are in the same trend with the proteomics results. It is reported that a higher level of ficolin is related to autoimmune disorders [40]. There have been reports about complement-associated and prostaglandin-dependent bone metabolism and physiologic remodeling [41], such as axial spondyloarthritis [42]. It is also believed that complement activation can release cytokines, proteases, and ROS, thereby contributing to inflammation and generates deleterious effects [41]. But so far, none of these molecules have been reported in the pathology of microtia.

Cadherin-13 (CDH13) that repeat in the groups A vs. B and A vs. D is an important adhesion molecule that mediates cell–cell adhesion in neural tube development [43], while there is a report about CDH13 and CRTAC expression in the limb mesenchyme control contrarily proprioceptor identity [44]. Coincidently, CRTAC is downregulated in our WES result. Cdh13 protects against cardiac stress through its association with adiponectin through the AMPK signaling pathway in mice [45]. What is more, research has proved that, Cdh13 serves as an anchor for tethering of adiponectin to M2 macrophages and can promote the cell proliferation by activation of Akt [46]. In our proteomic results, CDH13, adiponectin and CRTAC are all downregulated, which may imply the disruption of the cell proliferation and neural tube development, and inflammation may participate in the malfunction.

Besides, the commonly dysregulated genes were significantly enriched in GO terms in Fig. 3 and Table 9 shows that most proteins are significantly enriched in lipid metabolic process, peroxidase activity, and the hemoglobin related functions. LEPR⁺ MSCs give rise to most bone and adipocytes in adult bone marrow and can form bony ossicles that support hematopoiesis in vivo [47]. ApoE is involved in lipoprotein metabolism and ApoE -/- mice showed intervertebral disc degeneration partially due to increased matrix metalloproteinase (MMP) -3, -9, and -13 levels [48]. In addition, Farnaghi et al. found that hypercholesterolemia caused the degeneration of articular disc cartilage in osteoarthritis by the production of mitochondrion-derived oxidative stress [49]. Evidence shows that there is a link between oxidative stress and cartilage degradation.

What is more, emerging evidence has shown that cell metabolism influences gene expression by modifying the epigenome, which can regulate stem cell pluripotency, differentiation and somatic cell reprogramming [50]. ALDOB participates in fructose and other central carbon metabolism by reversibly cleaving fructose 1,6-diphosphate into glyceraldehyde 3-phosphate (G3P) and dihydroxyacetone phosphate (DHAP) [51]. Lipogenesis requires oxidative phosphorylation to provide a large amount of energy for lipid storage, while chondrocytes mainly rely on glycolysis for energy supply [52]. We hypothesize that reduced ALDOB protein levels in microtia patients may lead to decreased glycolysis and decreased MSC’s chondrogenic ability in some way which would be studied in future studies.

In label-free proteomics analysis, ALDOB protein is 5.41239 times higher in the patient compared to his sister(B) and 3.87038 times higher when compared to the mother in the A/C group, CDH13 and ADIPOQ are 0.35024 and 0.26241 times in A vs. B group respectively and show the same trend in A vs. D group by 0.31112 and 0.29045 time respectively. This is coherent with the result in IHC that CDH13 and ADIPOQ are lower in microtia patients. But the real-time PCR of other microtia patients shows that mRNA levels of ALDOB, ADIPOQ, CDH13, and TASP1 are all higher in the microtia group while IHC results show all of the four proteins are significantly lower when compared to the control group as shown in Figs. 7 and 8. This is probably due to the translational modifications of the mRNAs and shall be further explored as a candidate disease mechanism.

Fig. 8

Graphical Abstract. A graphical overview of this study

The collective findings provide novel insights into the pathogenesis of congenital microtia-atresia. It can be seen that the screened proteins are mainly involved in lipid metabolism, immunology, and inflammation process, and our in vitro verifications have proved that ROS do exist in microtia patients, as some fault programmed gene expression and misregulated microenvironment during the embryonic stage could be kept in the entirety and outer ear keeps on developing after birth [53], therefore these faults may affect ear length, width, prominence and degree of auricular development. The relationships between oxidative stress and these proteins in cancers or inflammatory processes have been widely studied [29], but the relationship between oxidative stress and chondrogenic defect of microtia has not been thoroughly realized. Alterations of the redox state in microtia may provide new insight into pathogenic mechanisms and refine therapeutic strategy, with the ultimate goal of improving the quality of life for the patients.

Potential limitations

There are also some limits of proteomic analyses as protein may not fully reflect the genetic variants and can only show the pathology of micro-environmental effect in this case. To overcome this limitation, we use the whole-exome sequencing to find out the possible variations in microtia at the genome level. However, the exome region constitutes only 1–2% of the entire genome, and WES inevitably leaves out the non-coding mRNA elements and/or intronic or intergenic regulatory regions that may be responsible for the disease phenotypes [2]. What is more, due to individual variances, more microtia patients are needed to further verify our findings. Last but not least, more explorations of how oxidative stress and malfunction of genes related to embryonic development in vivo shall be done in our future studies.

Conclusion

In summary, based on comprehensive analyses of WES from the patient with microtia and costal chondrocyte dysplasia, the present study found out SNPs including PRKDC, LEPR, PLG, MUC6, DGKK, GPR161, WDR36, BPTF, LRBA, SOD3, KCNQ4 and de novo CNVs in ESPNP may be the candidate genes that may function at the onset of the disease. Besides, there are likely deleterious CNVs including ASMT, DHRSX, ASMTL, AKAP17A, and P2RY8. Multi-omics data find out proteins with lipid metabolism, immune response, the complement cascades pathway and oxidative stress may take part in the pathogenesis of this rare case of congenital microtia. The verification of the multi-omics results with other microtia patients identified that ALDOB, ADIPOQ, CDH13 and TASP1, and oxidative stress may play a role in chondrogenic anomalies. In view of this, we speculate the importance of the interrelationship between gene network and environmental factors on the disease. Up till now, there is no further study examining the ROS pathogenic mechanism involved in microtia. A deeper understanding of microtia etiology will not only shed light on the occurrence of this defect, but also facilitate better prevention and therapeutic strategies for microtia.

Methods

Patients and ethics

A patient showing congenital microtia-atresia was hospitalized at The First Affiliated Hospital of Xi'an Jiao Tong University. Family history was obtained and clinical features of all the study subjects and evaluated via a review of medical records. This study was approved by the ethics committee of the First Affiliated Hospital of Xi'an Jiao Tong University Medical College and a written informed consent was provided. Authorization was also obtained for the disclosure of recognizable individuals in photographs and the collection of blood for further analyses.

WES

Whole blood (3 ml) was obtained from each individual. Genomic DNA extracted from peripheral blood for each sample was fragmented to an average size of 180 ~ 280bp and subjected to DNA library creation using established Illumina paired-end protocols. Whole-exome capture was performed using an Agilent SureSelect Human All Exon V6 enrichment capture kit (Agilent Technologies, Santa Clara, CA, USA). The Illumina Novaseq 6000 platform (Illumina Inc., San Diego, CA, USA) was utilized for genomic DNA sequencing in Genechem Bioinformatics Technology Co., Ltd (Beijing, China) to generate 150-bp paired-end reads with a minimum coverage of 10 × for ~ 99% of the genome (mean coverage of 100 ×). After sequencing, basecall files conversion and demultiplexing were performed with bcl2fastq software (Illumina). The resulting FASTQ data were submitted to in-house quality control software for removing low quality reads, and then were aligned to the reference human genome (hs37d5) using the Burrows-Wheeler Aligner (bwa), and duplicate reads were marked using SAMtools. SNP/INDEL calling: Single nucleotide variants (SNVs) and indels were called with samtools to generate gVCF. The raw calls of SNVs and INDELs were further filtered with the following inclusion thresholds: 1) read depth > 4; 2) Root-Mean-Square mapping quality of covering reads > 30; 3) the variant quality score > 20. CoNIFER (V0.2.2) was used to detect CNVs. Annotation was performed using ANNOVAR (2017 June 8). Annotations included minor allele frequencies from public control data sets as well as deleteriousness and conservation scores enabling further filtering and assessment of the likely pathogenicity of variants. The classification system of the American College of Medical Genetics and Genomics (ACMG) was used to predict the harmfulness of variations.

Label-free proteomics

Serum pools were depleted of the most abundant proteins from each individual using the Agilent Human 14 Multiple Affinity Removal System Column following the manufacturer’s protocol (Agilent Technologies). The 10 kDa ultrafiltration tube (Sartorius) was used for desalination and concentration of low-abundance components. One volume of SDT buffer (4%SDS, 100 mM Tris–HCl, pH 7.6) was added, boiled for 15 min and centrifuged at 14000 g for 20 min. The supernatant was quantified with the BCA Protein Assay Kit (Bio-Rad, USA). The sample was stored at -20 °C. 20 µg of proteins for each sample were mixed with 6X loading buffer respectively and boiled for 5 min. The proteins were then separated on 12.5% SDS-PAGE gel. Protein bands were visualized by Coomassie Blue R-250 staining. 200 μg of proteins for each sample went filter-aided sample preparation (FASP digestion) and LC–MS/MS analysis was performed on a Q Exactive Plus mass spectrometer (Thermo Fisher Scientific) that was coupled to Easy nLC (Thermo Fisher Scientific). 2μg peptide was loaded onto the C18-reversed phase analytical column (Thermo Fisher Scientific, Acclaim PepMap RSLC 50um X 15 cm, nano viper, P/N164943) in buffer A (0.1% formic acid) and separated with a linear gradient of buffer B (80% acetonitrile and 0.1% Formic acid) at a flow rate of 300 nl/min. The linear gradient was as follows: 5% buffer B for 5 min, 5–28% buffer B for 90 min, 28–38% buffer B for 15 min,38–100% buffer B for 5 min, hold in 100% buffer B for 5 min. MS data was acquired using a data-dependent top10 method dynamically choosing the most abundant precursor ions from the survey scan (350–1800 m/z) for HCD fragmentation. MS1 scans were acquired at a resolution of 70,000 at m/z 200 with an AGC target of 3e6 and a maxIT of 50 ms. MS2 scans were acquired at a resolution of 17,500 at m/z 200 with an AGC target of 2e5 and a maxIT of 45 ms, and the isolation width was 2 m/z. Only ions with a charge state between 2–6 and a minimum intensity of 2e3 were selected for fragmentation. Dynamic exclusion for selected ions was 30 s. The normalized collision energy was 27 eV. The MS data were analyzed using MaxQuant software version 1.6.14.0. MS data were searched against the Uniprot_HomoSapiens_20386_20180905 database, downloaded on http://www.uniprot.org/. An initial search was set at a precursor mass window of 6 ppm. The search followed an enzymatic cleavage rule of Trypsin/P and allowed maximal two missed cleavage sites and a mass tolerance of 20 ppm for fragment ions. Carbamidomethylation of cysteines was defined as fixed modification, while protein N-terminal acetylation and methionine oxidation were defined as variable modifications for database searching. The cutoff of the global false discovery rate (FDR) for peptide and protein identification was set to 0.01. Protein abundance was calculated on the basis of the normalized spectral protein intensity (LFQ intensity). |log2fold change (FC) |≥ 1 and p value < 0.05(Student’s test) were selected as the threshold for screening.

Real-time polymerase chain reaction (PCR)

Total RNA was extracted using TRIzol reagent (Invitrogen, USA). The total RNA was reversely transcribed using oligo dT (Promega, USA), M-MLV reverse transcriptase (Promega, USA) and dNTPs (Sigma–Aldrich, Germany) according to the manufacturer’s instructions. Real-time PCR was performed with a Light Cycler 480 SYBR Green Master and Light Cycler 96 system (Roche, Switzerland) according to the manufacturer’s instructions. Primers were designed by Primer Premier 5 software (Canada) and are listed in the Supplementary Table S5. The relative expression of each gene was normalized by comparison with the reference GAPDH and measured with the 2-△△CT method.

Immunohistochemical staining

Samples from the patients’ and normal people’s cartilage pieces were fixed in 4% paraformaldehyde for 24 h prior to embedding in paraffin. The samples were sectioned into 5-μm slices, mounted on glass slides, and stained with hematoxylin and eosin (HE). The primary antibodies used were mouse anti-human monoclonal antibody (ADIPOQ 1:200, 21,613–1-AP; ALDOB 1:200,18,065–1-AP; CDH13 1:50, 12,618–1-AP; TASP11:50, 16,739–1-AP; all from Proteintech, China), followed by incubation of biotin-labeled second antibody(Cat# kit-9730; MXB, Fuzhou, China) at RT for 50 min followed by enzyme-conjugated streptavidin (Cat# kit-9730; MXB, Fuzhou, China) incubate at RT for 50 min. Develop color with 3,3’-diaminobenzidine (Cat# DAB-1031; MXB, Fuzhou, China) and wash in water. Finally, the sections are stained with hematoxylin staining. The slides were then examined and photographed using the Leica SCN400 slide scanner and the integrated density was analyzed by Image-Pro Plus 6.0 software.

Cell culture

Auricular cartilage of microtia patients or normal people who had done rhinoplastic surgery was washed in PBS + 100 U/mL penicillin G, 100 μg/mL streptomycin (Gibco, Gaithersburg, MD, USA) and fragmented into 1mm³ pieces in 0.2% type IV collagenase (Gibco, Gaithersburg, MD, USA) and digested overnight at room temperature. The cell suspension was centrifuged at 1000 rpm for 5 min and filtered by 70 μm mesh (BD, USA) and cultured in the complete medium consisted of DMEM (Gibco, Gaithersburg, MD, USA) containing 10% heat-inactivated fetal bovine serum (BIOEXPLORER, USA). Chondrocytes were maintained in a humidified incubator containing 5% CO² at 37 °C and the first or second passages were used for the following experiments.

ROS analysis

Intracellular ROS levels were measured using 2’,7’-dichlorofluorescin diacetate (DCFH-DA; Cat#S3300S, Beyotime, China) according to the manufacturer’s instructions. Briefly, ear chondrocytes of microtia patients, normal people, ALDOB^OX and ALDOB^pEN groups were seeded to 48-well plates. Rosup was added as positive control, and one well was left as the blank control. DCFH-DA/Rosup was diluted with the serum-free medium at 1:1000 to achieve a final concentration of 10 µM and were incubated with cells for 30 min at 37 °C in the dark. The cells were washed 3 times with PBS to fully remove the DCFH-DA that did not enter the cells. The fluorescence intensity was detected with the excitation wavelength of 488 nm and the emission wavelength of 525 nm using the microplate reader (Cytation 5 M, Bio-Tek, USA).

Mitochondrial membrane potential detected by JC-1

Mitochondrial transmembrane potential (Δ Ψ m) was detected using the JC-1 mitochondrial membrane potential assay kit (Cat#C2005, Beyotime Biotech, China), following the manufacturer’s protocol. Cells cultured in 6-well plates were incubated with 1 ml JC-1 staining solution (10 µg/mL) at 37℃ for 20 min, and then the supernatant was removed and washed twice with JC-1 dyeing buffer (1X). Loss of mitochondrial membrane potential (Δ Ψ m) was assessed by fluorescence microscope (Axio Observer 7, Zeiss, German).

Statistics analysis

Results are expressed as the mean ± SD (standard deviation) with at least 3 independent biological repeats. The unpaired, two-tailed Student’s t-test was used to analyze the difference between groups using Prism 6 (GraphPad) software. P < 0.05 was considered as significant.