This article is mentioned in:

Introduction

Cartilage-hair hypoplasia (CHH; OMIM #250250) (1) is an autosomal recessive disorder characterized by disproportionately short stature, sparse hair and immunodeficiency (2). Manifestations can be highly variable among individuals with CHH; even siblings of the same genotype can exhibit different phenotypes (3). Moreover, patients with CHH often have reduced life expectancies. Two prospective follow-up cohort studies revealed that patients with CHH may develop immunodeficiency or malignancy in adults without immune defects (4), or are prone to recurrent respiratory tract infections, which contribute significantly to mortality (5). CHH is caused by mutations in the RNA component of mitochondrial RNA processing endoribonuclease (RMRP; OMIM 157660) gene, which encodes the RNA subunit of the RNase MRP complex and is a long non-coding RNA (lncRNA) (6). Pathogenic variants of RMRP exhibit a high degree of heterogeneity and can be broadly classified into two categories. The first category encompasses insertions, duplications, or triplications, which are often located in the region between the TATA box and the transcription initiation site. The second category comprises single nucleotide substitutions or alterations involving no more than two nucleotides, typically found within the transcribed region (7–9). It has been conclusively established that compound heterozygous or homozygous variants in the RMRP gene can lead to functional impairment, such as inhibition of ribosome synthesis (10), cell cycle regulation (11,12), promoter efficiency and RNA transcript instability (13). However, the mechanisms underlying the relationship between the RMRP gene and CHH remain unclear.

Instead of focusing solely on the mutation itself, RNA-sequencing offers an alternative perspective for exploring the pathogenic mechanisms of RMRP variants in CHH. Two relevant studies utilizing this technology have been retrieved from the literature. In one study, primary fibroblasts were isolated from patients with CHH and healthy control donors by skin biopsy. The fibroblasts were then subjected to modified single-cell tagged reverse transcription sequencing (STRT-seq) after several passages. The results revealed that CHH fibroblast cells have a slower growth rate and are specifically delayed in the passage from the G2 phase to mitosis (14). In another study, a fibroblast-derived chondrocyte (FDC) model combined with transcriptome sequencing demonstrated that fibroblasts from patients with CHH exhibit a reduced commitment to terminal differentiation. Some key factors in the bone morphogenetic protein, fibroblast growth factor (FGF) and insulin-like growth factor-1 signaling axes are significantly upregulated in CHH fibroblasts during their transformation into chondrocytes (15). Despite encouraging results, the mechanistic link between RMRP mutations and changes in gene expression remains unclear. The sequencing data from these two previous studies still hold value for further analysis due to the good representativeness of their samples. The present study utilized Ingenuity Pathway Analysis (IPA), which incorporates a curated collection of information from the biomedical literature and various databases (16), to conduct an in-depth exploration of the mechanism underlying CHH using data from the two aforementioned studies.

Materials and methods

Participants

A family of four Chinese individuals was recruited for the present study at Shenzhen Maternity and Child Healthcare Hospital (Shenzhen, China). Peripheral venous blood was collected from the 5-year-old patient, who was conclusively diagnosed with CHH after genetic testing in the present study, as well as from three other healthy family members: A 31-year-old father, a 29-year-old mother and a 3-year-old sister. The blood samples were collected in tubes containing EDTA as an anticoagulant and stored at −80°C. Genomic DNA (gDNA) was extracted from the venous blood of each individual using the commercial DNeasy Blood and Tissue Kit (cat. no. 69506; Qiagen GmbH). The present study was approved by the Medical Ethics Committee of Shenzhen Health Development Research and Data Management Center (Shenzhen, China; approval no. 2019-016) and was performed in accordance with The Declaration of Helsinki. Written informed consent was obtained from each participant or, in the case of a minor, from their parents.

Whole-exome sequencing (WES)

The exome library was prepared using the Ion AmpliSeq Exome RDY Kit (cat. no. A38262; Thermo Fisher Scientific, Inc.), according to the manufacturer's protocol. Briefly, gDNA was quantified using the Qubit dsDNA HS Assay Kit (cat. no. Q32851; Thermo Fisher Scientific, USA) on the Qubit 2.0 Fluorometer (Thermo Fisher Scientific, Inc.) and subsequently utilized for exome library target amplification. The libraries were purified using AMPure XP (cat. no. A63881; Beckman Coulter, Inc.). Specific Ion Xpress Barcode Adapters (cat. no. 4471250; Thermo Fisher Scientific, Inc.) were ligated to the amplicons at 22°C for 30 min, followed by 72°C for 10 min. The libraries were purified with AMPure XP. Subsequently, the barcoded exome libraries, with a final concentration of 8 pM, were loaded onto the Ion PI Chip Kit V3 (cat. no. A26770; Thermo Fisher Scientific, Inc.) for WES sequencing using the Ion Torrent Proton platform (Thermo Fisher Scientific, Inc.). This platform performs single-end sequencing and has no fixed read length; the read length is related to the electrochemical signal attenuation caused by the characteristics of the sequence being tested.

WES data analysis

The Ion Torrent platform transformed the electronic sequencing signals into raw DNA sequence data, which were subsequently analyzed using the Ion Torrent Suite V4.4 (Thermo Fisher Scientific, Inc.). An embedded TMAP tool was used to automatically align them with the human genome assembly (GRCh37/hg19). After alignment, a series of criteria were created for variant annotation. Only samples with >200× mean depth, >98% coverage of the designed region and >90% uniformity passed quality control. All minor variants were called using the embedded plugin TVC with default parameters (germline, low stringency). All variants were annotated with information such as gene location (RefSeq; http://www.ncbi.nlm.nih.gov/refseq/), gene function, population frequency [1000 Genomes (https://www.internationalgenome.org/), ExAC and gnomAD (https://gnomad.broadinstitute.org/)], impact prediction [dbNSFP (https://www.dbnsfp.org/) and InterVar (https://wintervar.wglab.org/)] and disease databases [ClinVar (https://www.ncbi.nlm.nih.gov/clinvar/) and HGMD pro (https://www.hgmd.cf.ac.uk/ac/index.php/)] using the ANNOVAR (https://annovar.openbioinformatics.org/en/latest/) pipeline. An in-house script was used to evaluate the pathogenicity of each variant based on the guidelines established by the American College of Medical Genetics and Genomics and the Association for Molecular Pathology (17,18). Consequently, all variants were classified into five categories: Pathogenic, likely pathogenic, uncertain significance, likely benign or benign. Finally, only mutations classified as pathogenic or likely pathogenic were deemed harmful and confirmed using Sanger sequencing. Furthermore, a freely available tool, Ensembl Variant Effect Predictor (VEP; http://asia.ensembl.org/info/docs/tools/vep/index.html), was employed for the annotation and prediction of the effects of genomic variants. WebLogo 3.7.11 (https://weblogo.threeplusone.com/) (19) was used to analyze the variant conservation.

Sanger sequencing

Sanger sequencing was used to validate and evaluate the co-segregation of RMRP variants in all family members. The following primers were used for PCR amplification: Forward, 5′-CAACAGGTGAAAATCCGTCTC-3′ and reverse, 5′-TGCCTCTGAAAGCCTATAGTCT-3′. PCR was performed in a 25-µl volume, using ~25 ng gDNA as the template in a reaction with PCR Master Mix (cat. no. M7502; Promega Corporation). The thermal cycling amplification procedure consisted of 2 min at 95°C for pre-denaturation, followed by 35 cycles of amplification at 95°C for 15 sec, annealing at 58°C for 15 sec and extension at 70°C for 15 sec. The reaction was completed with a final extension step for 8 min at 70°C. The PCR products were subsequently purified using the QIAquick PCR Purification Kit (cat. no. 28104; Qiagen GmbH). Finally, Sanger sequencing was performed using the ABI 3730XL DNA analyzer (Thermo Fisher Scientific, Inc.).

Bioinformatics analysis

The first sequencing data profile was downloaded from EMBL-EBI ENA (accession no. PRJEB23608; http://www.ebi.ac.uk), including five adults with CHH (all were RMRP g.70A>G homozygous) and five matched controls. Fibroblasts from skin biopsies were passaged 2–5 times and underwent STRT-seq (14). The DEGs used for subsequent analyses were sourced from the article that originally published the PRJEB23608 sequencing data. These DEGs were identified by the author based on a differential-expression P-value <0.05 and q-value <0.05 (14). The second set of RNA-sequencing data was downloaded from EMBL-EBI ArrayExpress (accession no. E-MTAB-10996; http://www.ebi.ac.uk/biostudies/arrayexpress), containing four patients with CHH (one RMRP g.70A>G homozygous and three RMRP compound heterozygous variants) and four control donors. Fibroblasts from skin biopsies underwent FDC transdifferentiation and were used for RNA-sequencing (15). Data from day 3 were selected, and DEGs were filtered in the present study using the DESeq2 R package (https://bioconductor.org/packages/release/bioc/html/DESeq2.html) with the following criteria: Fold change ≥2 and false discovery rate <0.05. After common DEGs were identified from the two sequencing data profiles, a comprehensive analysis was conducted using the IPA knowledge base (content version: 111725566; release date: March 21, 2024; Qiagen GmbH). The ‘Core Analysis’ module was executed to select the most significant canonical pathways, upstream regulators, causal networks and biological functions based on a suite of implanted algorithms or score tools, such as an enrichment score (Fisher's exact test P-value) that measures the overlap of observed and predicted regulated gene sets, and a Z-score assessing the effects of molecular changes within a dataset on biological processes or functions (20). A detailed flowchart of the dataset analysis is shown in Fig. 1.

Figure 1. Detailed information of the two datasets and the data analysis flowchart. The PRJEB23608 dataset consisted of STRT-seq data from fibroblasts that were subcultured 2–5 times. Fibroblasts were derived from skin biopsies of 10 Finnish individuals, including 5 control donors and 5 patients with CHH carrying homozygous mutations. The E-MTAB-10096 dataset comprised RNA-seq data obtained from dermal fibroblasts cultured in fibroblast-derived chondrocyte medium for 3 days. Dermal fibroblasts were obtained from 8 German individuals, comprising 4 control donors and 4 patients with CHH and either homozygous or heterozygous mutations. CHH, cartilage-hair hypoplasia; RMRP, RNA component of mitochondrial RNA processing endoribonuclease; STRT-seq, single-cell tagged reverse transcription sequencing; RNA-seq, RNA-sequencing; IPA, ingenuity pathway analysis.

Results

Clinical manifestations of the patient with CHH

A patient with normal intelligence, a birth height of 47.0 cm and a birth weight of 2.9 kg was born at 38 weeks of gestation from an uneventful pregnancy. The patient was the first girl born to non-consanguineous healthy parents with a height below average (father, 163 cm; mother, 158 cm). The height and weight (56.0 cm and 6.2 kg) of the patient at 7 months were both below the 3rd percentile for age. There was a noticeable disproportionate short stature (height, 72.0 cm; weight, 9.5 kg) at 1.7 years of age, accompanied by sparse hair, short fingers and toes, and elbows that could not be straightened. Radiographs of the hand skeleton at 4.5 years old (height, 78.4 cm; weight, 11.1 kg) showed prominent thick and short metacarpal and phalangeal bones, and metaphyseal dominant chondrodysplasia. The left ulna was slightly bent (Fig. 2A). Growth hormone (GH) was injected 3 months later, with an initial dose of 2.5 U/day (6 days/week) and a maximum dose of 3.0 U/day maintained until 6.8 years old (height, 88.5 cm; weight, 13.8 kg). However, the radiographs were still unsatisfactory, showing lumbar scoliosis, upward warping of sacral vertebrae, uneven width of lumbar pedicle spacing, and abnormal bone in the metaphysis of the bilateral femoral neck (data not shown). As shown in the growth curves (Fig. 2B), the height and weight of the patient increasingly deviated from the normal range with age. Notably, treatment with GH for ~2 years resulted in a gradual improvement; however, the patient should be tracked for a longer period to observe the final effect.

Figure 2. Radiographs and clinical manifestations of the patient with cartilage-hair hypoplasia. (A) Main clinical features of the proband. (Aa) Radiograph of the left hand, taken at the age of 4.5 years showed short metacarpals and phalanges with bullet- and cone-shaped epiphyses. (Ab) Radiographs of the left arm, taken at the same age, revealed metaphyseal chondrodysplasia at the elbow epiphysis and shortening of the distal ulna. (Ac) The patient still had a disproportionately short stature at 6.8 years of age, (Ad) exhibited sparse hair, (Ae) had short and thick fingers, and (Af) short toes and feet. (B) Growth curve of the proband from birth to 6.8 years old. The blue lines represent the growth standard for Chinese girls, and the numbers to the right of the blue lines indicate the percentiles for height and weight, based on data from the Capital Institute of Pediatrics (54); the red lines depict the growth curve of the patient with CHH. Their height and weight deviation increased with age before the injection of growth hormone.

Genetic causal analysis

Before arriving at the Institute of Maternal and Child Medicine Research (Shenzhen, China) for genetic testing, the proband (II:1) was examined at different hospitals or institutions, with a summary of their test results provided as follows: Achondroplasia (ACH), which is an autosomal dominant disorder caused by variants in the FGF receptor 3 gene, was highly suspected based on the clinical phenotype. However, hotspot variants (NM_000142.5:c.1138G>A and NM_000142.5:c.1123G>T) in the coding region of exon 10 were not detected using Sanger sequencing. In addition, there were no significant abnormalities in chromosomal karyotype or urinary organic acid levels. At 5.3 years old, the patient was diagnosed with congenital heart disease (atrial septal defect), and a single heterozygous mutation (NR_003051.4: RMRP n.197C>T) was detected in the father using a commercial exon sequencing kit that simultaneously detects 4,811 genes but does not cover the 10-base range at the exon-intron junction. No clinically significant chromosomal aberrations were detected using a CytoScan HD chip. Based on the results previously obtained from other institutions, WES was performed on all family members. Finally, two evolutionarily conserved variants of RMRP, NR_003051.4: n.-21_-2dup (previously known as NR_003051.3: n.-22_-3dup) and n.197C>T (previously known as NR_003051.3: n.196C>T), which form a novel compound heterozygous mutation, were identified in the affected patient (Fig. 3A-C). Variants and zygosity were confirmed using targeted Sanger sequencing. Both parents and the sister of the patient, with normal phenotypes, were heterozygous carriers of either variant (Fig. 3D). In addition to fitting the recessive inheritance pattern, the main clinical manifestations of the patient matched the phenotypic spectrum of CHH. Furthermore, Ensembl VEP also predicted that both mutations were pathogenic.

Figure 3. RMRP gene variants. (A) Schematic diagram illustrating the localization of mutations on the RMRP gene, with one mutation located in the promoter region and another in the transcription region. (B) Pedigree of the family with CHH. The arrow indicates the proband; black shapes represent patients with CHH; half-filled shapes represent carriers of RMRP variants. (C) Conservation analysis of variants using WebLogo 3.7.11. Each logo consists of stacks of symbols, one stack for each position in the sequence. The overall height of the stack indicates the sequence conservation at that position, while the height of symbols within the stack indicates the relative frequency of each nucleic acid at that position. (D) Sanger sequencing validation results in family members. Except for the proband, the family numbers were all carriers of a single heterozygous mutation. RMRP is located on the reverse strand. The arrow in the n.197C>T sequence indicates G instead of C because the reference sequence provided by Sanger sequencing is the complementary reverse strand. CHH, cartilage-hair hypoplasia; RMRP, RNA component of mitochondrial RNA processing endoribonuclease.

IPA reveals the mechanistic network of RMRP

A total of 31 common DEGs were intersected from the PRJEB23608 dataset (165 DEGs) and E-MTAB-10996 day-3 dataset (774 DEGs). As shown in Fig. 4A, ‘cell cycle checkpoints’ was the most prominently enriched canonical pathway based on a -log(P-value) and was revealed to be suppressed with a negative Z-score in the IPA core analysis. Notably, enrichment analysis of the two independent DEGs showed that this pathway also ranked first and was inhibited. The other pathways enriched by both were highly similar, such as the ‘kinetochore metaphase signaling pathway’, ‘mitotic prometaphase’ and ‘RHO GTPases Activate Formins’. Fig. 4B shows the top enriched diseases and functions, from which cancer, immunological diseases, respiratory diseases and developmental disorders were all markedly enriched. These results reaffirmed that CHH is a disease involving multiple organs. IPA comprehensive analysis was executed based on E-MTAB-10996 day-3 data because this dataset contained more DEGs and downregulated RMRP. The key node genes associated with ‘cell cycle checkpoints’, such as CCNE2, CDC25A, CHEK1 (CHK1), CCNB2, CDC25C and CCNA2, were all dysregulated in these data (Fig. 4C). These findings demonstrated that RMRP, as the major upstream regulator, may regulate the aberrant expression of downstream target genes, mainly by inhibiting the transcription factor TP53. The imbalanced expression of these genes could ultimately lead to the inhibition of ‘cell cycle checkpoints’ and consequently activate the function of ‘growth failure or short stature’, which are closely related to the occurrence of the CHH (Fig. 4C).

Figure 4. Enrichment analysis results of DEGs from different datasets and the RMRP regulatory network. The cells used for analysis were primary fibroblasts that were cultured following skin biopsy. (A) Dot plots showing the enrichment of canonical pathways. (B) Dot plots showing the enrichment of diseases and functions. (C) A schematic network of the RMRP regulatory mechanism. RMRP and the downstream DEGs, such as PLK1, CHEK1 and HK2, were colored according to their fold change in the E-MTAB-10996 (day 3) dataset, with green indicating downregulation and red indicating upregulation. Others were colored according to ingenuity pathway analysis prediction, with orange for activation and blue for inhibition. DEGs, differentially expressed genes; RMRP, RNA component of mitochondrial RNA processing endoribonuclease.

Discussion

The present study detected a new compound heterozygous variant in a Chinese girl with typical features of CHH. Two previous reports were retrieved related to the n.-21_-2dup variant. One was detected in a Japanese boy with CHH, along with g.218A>G (21); the other was detected in a German boy with CHH, along with 193G>A (8). Similarly, it has been reported that g.-19_-3dup (previously known as NR_003051.4: n.-18_-2dup), which is only 3 nt shorter than n.-21_-2dup variant, forms a compound heterozygous mutation with g.193G>A (9,22) or g.4G>T (23), leading to the occurrence of CHH. The n.197C>T variant of the paternal allele is a known pathogenic mutation evolutionarily conserved in the RMRP stem-loop pairing region (11,24). Moreover, Gomes et al (25) reported that >50% of Brazilian patients with CHH carry the n.196C>T variant, suggesting a possible founder effect. However, to the best of our knowledge, the novel compound heterozygous model formed by these two variants, as discovered in the present study, has not been previously reported. In addition to intrafamilial co-segregation validation, the Ensembl VEP (26) was used to predict the variation effects, and the results indicated that they were pathogenic in CHH. The parents and sister of the patient were asymptomatic; therefore, it may be concluded that the compound heterozygous variant was the genetic cause of the condition. Notably, immunodeficiency was not observed in this patient, although most patients with CHH usually have variable degrees of immune dysfunction (27–29). Collectively, these results reaffirmed the diversity of CHH and RMRP variants.

In the E-MTAB-10996 data profile, the mRNA levels of the RMRP gene were markedly reduced, which has also been confirmed by previous studies in patients with CHH with diverse variants (13,25,30). Whether the disease progresses to CHH or cancer may depend on RMRP expression (31). Decreased RMRP levels, caused by mutations, can lead to CHH. The possible underlying mechanisms include decreased rRNA processing, cell proliferation, changes in the transcriptome and increased Wnt/β-catenin signaling (8,10,13,14,31). Elevated RMRP levels have been observed in various types of cancer and are thought to contribute to malignancy by sequestering various microRNAs (31). In the present study, it was revealed that the downregulation of RMRP may regulate the aberrant expression of downstream targets by activating the transcription factor TP53. This could eventually lead to the inhibition of ‘cell cycle checkpoints’ and the manifestation of phenotypes such as ‘growth failure of short stature’. The tumor protein TP53, which encodes the cellular tumor antigen p53, is a well-known transcription factor that serves a central role in maintaining genome stability under various stress signals and determines the outcome of the DNA damage checkpoint response (32,33). The increased expression of p53 may be related to elevated β-catenin in synovial tissues, which blocks p53 proteolysis (34–36). In a zebrafish model of CHH, a mutation in RMRP has been reported to disrupt chondrogenesis and bone ossification through enhanced Wnt/β-catenin signaling (37). Excessive activation of the canonical Wnt/β-catenin signaling pathway may contribute to phenotypic instability in chondrocytes and the loss of cartilage homeostasis (38–40). Cartilage homeostasis is essential for sustaining the proper phenotype and metabolism of chondrocytes. Disruption of extracellular matrix formation and breakdown processes can lead to the loss or dysfunction of cartilage homeostasis, promoting the development of degenerative cartilage disorders (41–43). Previous studies have established that long non-coding RNAs serve important roles in cartilage development, degeneration and regeneration (44–46). In the final stage of cartilage development, the lncRNA RMRP promotes the differentiation of chondrocytes into hypertrophic chondrocytes. Conversely, interference with RMRP leads to the deregulation of chondrogenic differentiation (46–48). Cell cycle checkpoints are critical for enabling an orderly cell cycle, responding to irreparable DNA damage and maintaining genomic stability during cell division. Based on their distinct functions, cell cycle checkpoints are classified into two groups: DNA damage checkpoints (ATM/CHK2/p53) and DNA replication stress checkpoints (ATR/CHK1/WEE1), which are involved in the surveillance of the G1/S and G2/M checkpoints (49–51). The present study revealed that most of the key node genes associated with ‘cell cycle checkpoints’ were dysregulated in patients with CHH. In summary, the present study indicated that mutated RMRP may act as a stress signal to activate TP53 by elevating β-catenin or Wnt/β-catenin signaling, consequently resulting in cell cycle arrest and chondrodysplasia. Moreover, multiple systemic diseases were enriched in both transcriptome datasets assessed in the present study and were associated with previously reported CHH follow-up outcomes, such as malignancy, immunodeficiency and respiratory disease (4,5), which may explain the variable phenotypes of CHH.

There are currently no standard treatments for CHH. A previous study exhibited that a total of 7 years of GH treatment markedly improved bone growth and had a positive effect on growth rate; however, the height velocity decreased with the interruption of GH treatment (21). In another case, the height SD score of the patient changed from −4.00 to −2.98 after 4 years and 7 months of treatment (52). GH is recommended in cases of persistent short stature and small for gestational age. Gonadotropin-releasing hormone agonist (GnRHa) treatment may be considered when a short adult height is expected at pubertal onset (53). Albeit with only 2 years and 4 months of follow-up, the present study also observed a positive effect of GH on height, which was the third case retrieved regarding GH treatment for CHH. Therefore, it may be concluded that GH can be utilized as a treatment for CHH and, if possible, GnRHa can be used as an adjunct to treatment.

In conclusion, the present study identified a novel compound heterozygous pathogenic variant in a Chinese girl with CHH who exhibited a notable improvement in height following GH treatment. RMRP variations should be considered when a patient has typical short stature and the possibility of ACH has been excluded. In addition, two transcriptome datasets containing nine patients with CHH and homozygous or compound heterozygous mutations were selected for analysis. As the sampling and sequencing methods were similar, the results were highly comparable and representative. The findings of the present study provide valuable insights into the mechanism of action of RMRP in CHH and its clinical treatment. Next, we aim to investigate the function of RMRP variants through gene-edited cartilage lineage cells and/or mouse models. These experiments will verify the discovered mechanisms, and elucidate the relationship between mutated RMRP and cartilage development, encompassing degeneration, regeneration and homeostasis. In addition, we aim to assess a cohort of patients with CHH of multiple ethnicities to enhance the robustness of the present findings.

Acknowledgements

Not applicable.

Funding

The present study was supported by the Shenzhen Key Laboratory of Maternal and Child Health and Diseases (grant no. ZDSYS20230626091559006).

Availability of data and materials

The data generated in the present study may be requested from the corresponding author. The raw WES data generated in the present study may be found in the NCBI SRA under accession number PRJNA1165399 or at the following URL: https://www.ncbi.nlm.nih.gov/sra/?term=PRJNA1165399. The variants NR_003051.4:n.-21_-2dup and NR_003051.4:n.197C>T were submitted to the NCBI ClinVar under accession numbers SCV005324796 and SCV005328392 or at the following URLs: http://www.ncbi.nlm.nih.gov/clinvar/variation/550387/?oq=SCV005324796&m=NR_003051.4(RMRP):n.-21_-2dup and http://www.ncbi.nlm.nih.gov/clinvar/variation/633393/?oq=SCV005328392&m=NR_003051.4(RMRP):n.197C>T.

Authors' contributions

JG and SD contributed to laboratory investigations, data analysis, manuscript preparation, and manuscript drafting and revision. JG and SC performed experiments and managed the program. JG and JZ completed data acquisition, analysis and interpretation. JZ and SL assisted in the bioinformatics analysis and interpreted the data. JG, JZ and SD oversaw the research program and reviewed the manuscript. JG, JZ and SD confirm the authenticity of all the raw data. All authors read and approved the final version of the manuscript.

Ethics approval and consent to participate

The present study was reviewed and approved by the Medical Ethics Committee of the Shenzhen Health Development Research Center (approval no. 2019-016). All procedures involving human participants were performed following the ethical standards of the institutional or national research committee and in accordance with The 1964 Helsinki Declaration and its later amendments or comparable ethical standards. Written informed consent was obtained from all of the participants or the parents of minor participants.

Patient consent for publication

The parents of minor participants provided written informed consent for their personal or clinical details, along with any identifying images, to be published in this study.

Competing interests

The authors declare that they have no competing interests.

References