Patients were diagnosed and recruited based on typical clinical symptoms of mandibular, ear and facial soft tissue hypoplasia, physical examination results, and imaging reports including x-ray and CT reports. All samples were collected from Department of maxillofacial surgery, Plastic Surgery Hospital, Chinese Academy of Medical Sciences, Peking Union Medical College. All participants enrolled in our present study were Han Chinese, Asia population.

Six patients with HFM and four healthy controls were enrolled in our study. All patients' adipose tissue samples in our manuscript were collected through scar excision when removing the mandibular distractor surgery. The adipose tissue from controls was also collected through the buccal fat pad resection procedure from patients with facial fat accumulation. After taking adipose tissue away from operating table, we dissected it into small pieces with surgical scissors. Samples were maintained in centrifuge tubes with RNAlater (Ambion Inc) and stored at −80°C as the first part, preserved in cryotub.

According to the Declaration of Helsinki Principles, the protocols for the study were approved by the Ethics Review Committee of Plastic Surgery Hospital. Written informed consent was obtained from participants.

According to the manufacturer's instructions, total RNA was isolated from adipose tissues using the RNeasy Lipid Tissue Mini Kit (QIAGEN). In accordance with the manufacturer's instructions, stranded total RNA LT sample prep kit (Illumina) was used to prepare the compatible library. A Bioanalyzer (Agilent) was used to analyze the quality and concentration of all libraries. Sequencing of mRNA was performed on an Illumina Hiseq 2,500 sequencing system (Illumina), and 150-bp paired-end FASTQ files were generated.

Total RNA was isolated with the TRIzol reagent (Thermo Fisher Scientific) according to the manufacturer's protocol. The RNA concentration and purity were evaluated with NanoDrop ND-2000 spectrophotometer (Thermo Fisher Scientific Inc). With 1 µg of RNA in the 20 µl reaction system, reverse transcription reactions were performed with PrimeScript Reverse Transcriptase (Takara Bio) according to the manufacturer's instructions.

To validate the confidence of RNA-Seq, several differentially expressed genes were selected and analysed by quantitative real-time PCR (qPCR) utilizing the SYBR Premix EX Taq reagent (Takara) in a QuantStudio 7 Flex Real-Time PCR System (Applied Biosystems). Primers (Supplementary Table S1) were designed for the coding sequences of the candidate genes in Primer 3 software (http://frodo.wi.mit.edu/cgibin/primer3). GAPDH was used as the internal control. qPCR replicates were performed in a final volume of 10 µl containing primers, SYBR Premix EX Taq reagent (Takara) and cDNA templates. All quantitative PCRs were performed for three biological replicates. The relative expression levels of the candidate genes were calculated as the averaged normalized Ct value of each sample compared with the GAPDH Ct value of the corresponding sample based on the 2−ΔΔCt method.

Our RNA expression levels were estimated using HISAT2 (version 2.1.0). Genome sequence data (hg19, Genome Reference Consortium GRCh37) and annotation data were obtained from UCSC website (https://genome.uscs.edu). With featureCounts (v1.5.0-p3), we calculated transcript counts at the gene level and relative abundances in FPKM (fragments per kilobase of exon per million fragments mapped). We build Principal Component Analysis (PCA) and identified 1,234 differential expressed genes (DEGs) by using DEseq2 (P-value <0.05, log2Foldchange >2) (8). Heatmap show the gene expression pattern in different groups, all DEGs expression level were scale to −2 and 2 based on Z-score calculation. Diagrams were visualized by ggplot2 in R version 4.1.0.

The KEGG (Kyoto Encyclopedia of Genes and Genomes database) and gene ontology (GO) pathway enrichment analysis was performed using ClusterProfiler package from BioConductor (http://www.bioconductor.org/). A hypergeometric exact test was used to classify enrichment GO categories and KEGG pathways. Only terms with a P-value 0.05 were considered enriched. According to the rank of P-value, we display the top 20 GO enriched terms for both upregulated and downregulated genes.

The buccal fat pad from healthy controls were collected in D-Hanks solution with 1,000 U/L penicillin streptomycin (Solarbio, Beijing, China). An equal volume of 0.25% of EDTA trypsin supplemented with 0.1% of type I collagenase solution was added, and the samples were gently vibrated in a 37°C incubator shaker for 45 min. Cells were collected and centrifuged at 1,000 rpm for 10 min; the supernatant was removed, and the cells were resuspended with DMEM supplemented with 10% FBS. After filtering using a 100-mesh sieve, cells were transferred into a culture dish and incubated in 37°C with 5% CO 2 for 24 h. Cells were observed using a phase-contrast microscope daily. Cell morphology and proliferating rates were recorded, and when cells reached 90%–95% confluence, cells were digested with 0.25% EDTA trypsin and subcultured into several culture dishes. Lentiviral interference vectors, pCMV.GFP&PR.U6 (sh-, interference vector), and lentiviral overexpression vector, VP001-CMV-MSC-3flag-EF1-ZsGreen-T2A-PURO (oe-, overexpression vector), were purchased from Sangon Biotech (Shanghai, China) and General Biosystems (Anhui,China). A lentiviral-based HOXB2 interference vector, a HOXB2 overexpression vector were constructed. Lentiviral infections were performed using protocols available online (www.broadinstitute.org). Next, cells were inoculated into a six-well plate at a density of 5 × 104 cells/ml and infected with NC, si-HOXB2, vector, HOXB2-OE, respectively.

For the cell proliferation test, a CCK8 Kit (Solarbio, China) was used. In 96-well plates, cells were seeded in 96-well flat-bottomed plates with each well containing 2,000 cells in 200 µl of culture medium and cultured in ambient environment described above. CCK8 reagent was added into wells, after incubation at 37°C for 24 h in a humidified incubator with 5% CO2, the proliferative ability of the cells was measured at 450 nm. For transwell migration assay, different cells were suspended in 200 µl of DMEM without FBS and seeded on the top chamber of 24-well plate-sized Transwell inserts (Corning Falcon). The lower chambers contained DMEM with 20% FBS. After incubation for 8 h, the inserts were fixed and stained with crystal violet. Cells in the upper chamber were removed with cotton swabs. The average confluence of migrated cells was analyzed by ImageJ according to three random fields captured by 100× microscope. Each experiment was conducted in triplicate.

The 6 patients with HFM enrolled in our study who were aged from 7 to 26. All patients were classified with Pruzansky II and III mandibular hypoplasia, including facial soft tissue deficiencies, shortened mandibular ramus, small glenoid fossa, malformed condyle, preauricular tags, microtia, excluding other concurrent clinical symptoms (epibulbar dermoids and vertebral anomalies). The 4 matched controls were all adults who aged from 20 to 25. The adipose tissue was collected when matched controls underwent buccal fat pad resection. The clinical details of the participants are provided in Supplementary Table S2.

A transcriptomic analysis of adipose tissue from six patients with HFM and four healthy control was conducted to better understand the pathogenesis of the disease. A total of 28 million read pairs were obtained from each sample and were compared between the HFM and control groups in adipose tissue. The total number of annotated mRNAs identified was 333,341 and 1,244 mRNAs were significantly deregulated (q-value <.05; absolute value Log2 Ratio ≥1) (Figure 1A). A comparison of the adipose tissues from patients with HFM and those from controls demonstrated that 710 genes were up-regulated, and 534 genes were down-regulated. The differentially expressed mRNA data sets were analyzed using principal component analysis and hierarchical clustering. HFM adipose tissue exhibited distinct gene expression profiles compared with control adipose tissue based on PCA (Figure 1B). According to a hierarchical clustering analysis of the differentially expressed genes in adipose tissue from patients with HFM and control subjects, each gene expression pattern clustered separately (Figure 1C). Subsequently, we draw a correlation heatmap between various gene expression (Figure 1D). The results indicated that HFM has a profound impact on the expression of mRNA in the facial adipose tissue.

To confirm the accuracy of RNA-Seq, five genes including up-regulated (HOXB2, HAND2, COL1A1, MACH1) and down-regulated (SIX2) between adipose tissue from patient with HFM and controls were selected for further validation in additional set of clinical samples. Consistent with our expectation, the validation results of HOXB2 and HAND2 were consistent with the RNASeq data. While no significant difference (P > 0.05) was detected in the levels of COL1A1, MACH1, SIX2 mRNA, according to the results of quantitative real-time PCR (Figure 2).

Gene ontology and pathway analysis of genes were performed to find possible biological alterations related to HFM. A wide range of biological functions and signaling pathways are involved with DEGs. The immune system and the skeletal system exhibited the most significant differences, out of many involved pathways.

Based on the enrichment score, we selected the top ten items from GO (Figure. 2). the most significant biological processes related to craniofacial morphogenesis included skeletal system development, ossification, regulation of vasculature development, embryonic skeletal system development and embryonic skeletal system morphogenesis. According to KEGG analysis (Figure. 3), signal pathways including PI3K−Akt signaling pathway and Human papillomavirus infection showed significant differences. The PI3K−Akt signaling pathway was previously shown to be associated with NCCs development (9). The relationship between human papillomavirus infection and HFM has not been reported yet. Additional studies are required to further investigate these possibilities.

To examine the effects of HOXB2 on ADSC proliferation and metastasis, we conducted a series of cell function experiments. HOXB2 overexpression significantly decreased the migration and invasion of ADSC (Figures 4A–D), as well as reduced their proliferation (Supplementary Additional file S1, Figure S2B). On the contrary, Loss of HAND2 further enhanced cell migration and invasion of these ADSC (Supplementary Additional file S1, Figure S2F, S3B–D).