This systematic review followed the publishing guidelines and checklist established by PRISMA (Preferred Reporting Items for Systematic Review and Meta-Analysis). Articles were included and excluded based on the following eligibility criteria: 1) Inclusion criteria: described genes causing or potentially associated with amelogenesis imperfecta and enamel hypoplasia in species other than humans; published as original articles (not as review articles, editorials, dissertations, conference proceedings, or comments); and published in the English language; 2) Exclusion criteria: gene mutations were not described in the original articles; enamel defects resulting from exposure to environmental risk factors; cell-based experiments, molecular and biochemical analyses, structural and component analyses, and evolutional researches; and the articles failed to fit in any of the above criteria but did not include amelogenesis imperfecta candidate genes or related information.

The search for articles was conducted through three central literature databases: Medline (Ovid), PubMed (National Library of Medicine), and Embase (Ovid). In addition, relevant articles were searched in Scopus (Elsevier) to retrieve any studies missed in the database searches. Concepts included in the search to identify studies were amelogenesis imperfecta and genetics (gene mutation). No specific species was included in the keywords since our review included all species. A combination of Medical Subject Headings (MeSH) terms and titles, abstracts, and keywords was developed to obtain the initial Medline search string, and then adapted to the searches of the other databases. The Mouse Genome Informatics (MGI) database was searched using keywords “amelogenesis imperfecta,” “enamel hypoplasia,” “tooth enamel,” “tooth mineralization,” and “enamel mineralization” in order to provide a means of comparison and validation for the systematic review and identify genes that were potentially missed in the database searches.

The citations searched were stored in Rayyan (https://rayyan.qcri.org/welcome), an online application for systematic reviews that stores the citations/results, automatically processes the removal of duplicates obtained through various database searches, and tracks the decisions made during the systematic review. The primary Excel workbook designed for the systematic review (http://libguides.sph.uth.tmc.edu/excel_SR_workbook) was also used for tracking search strategies and results. A Cohen’s kappa test was conducted by two screeners to check the reliability of study selection during title and abstract screening. After achieving a >90% score for the Cohen’s Kappa test, all the titles and abstracts found through the database search were full-text reviewed by the two screeners independently. All the screening results were recorded in the Primary Excel workbook, and a codebook for data collection from eligible articles was developed as previously described (Sangani et al., 2015).

The Database for Annotation, Visualization, and Integrated Discovery (DAVID) (http://david.abcc.ncifcrf.gov/) was used for the gene set enrichment analysis. Gene Ontology (GO), including its Biological Process (BP), Molecular Function (MF), and Cellular Component (CC), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were used as reference gene sets (Sun H. et al., 2019). The top five most significant pathways or GO terms were selected for further analysis. k-means was used to cluster the gene functional enrichment results and the square error to extract the closest clusters. The highly-expressed mouse tooth miRNAs were retrieved from the publications (Cao et al., 2010). The miRNA-mAIGene regulations were integrated using the data from four databases: TargetScan (version 7.1) (Agarwal et al., 2015), miRanda (August 2010 Release) (John et al., 2004), miRTarBase (Release 7.0) (Huang et al., 2020), and PITA (version 6) (Kertesz et al., 2007). Considering the possibility of false results and multiple targets for each miRNA in these databases, the intersection of miRanda and PITA was merged with the intersection of TargetScan and miRTarBase to obtain reliable miRNA-mAIGene pairs. This conservative approach was demonstrated to effectively reduce the prediction of false-positive miRNA-mAIGene pairs (Jiang et al., 2016; Bonnet et al., 2020). Each gene set (GO term or KEGG pathway) containing at least two genes was used in the core miRNA family-based regulatory network. A Fisher’s exact test was applied to assess the enrichment significance of the miRNAs. All networks were visualized using Cytoscape (Shannon et al., 2003).

The mHAT9d mouse dental epithelial cell line originated from the apical bud of the incisors was a gift from Dr. Hidemitsu Harada (Iwate Medical University, Iwate, Japan). mHAT9d cells were cultured in Dulbecco’s Modified Eagle Medium: Nutrient Mixture F-12 (DMEM/F12; Thermo Fisher Scientific) supplemented with B-27 (Thermo Fisher Scientific), 25 ng/ml basic FGF (233-FB; R&D Systems), 20 ng/ml EGF (2028-EG; R&D Systems), and penicillin/streptomycin (Otsu et al., 2016). The LS8 cell line (Chen et al., 1992) was provided by Dr. Malcolm Snead (University of Southern California). Cells were plated at a density of 60,000 cells onto a 12-well cell culture plate and maintained until 80% confluence. The cells were treated with mimic for a negative control, miR-16-5p, miR-23a-3p, miR-23b-3p, miR-27b-3p, or miR-214-3p (mirVana miRNA mimic, Thermo Fischer Scientific) using Lipofectamine RNAiMAX transfection reagent (Thermo Fisher Scientific), according to the manufacturer’s protocol (24 pmol of mimic and 3 μL of transfection reagent in 1 ml of medium per well). After 24 h of treatment, the cells at 100% confluence were cultured with differentiation medium [including 15 μg/ml retinoic acid (R2625, Sigma Aldrich) and 0.1 μM dexamethasone (D4902, Sigma Aldrich)] in order to induce ameloblast differentiation.

mHAT9d cells were plated onto ibiTreat 8-well μ-slides (ibidi GmbH, Munich district, Germany) at a density of 10,000/chamber and cultured until 80% confluence. Cells were then treated with a mimic for miR-16-5p, miR-27b-3p, or control using Lipofectamine RNAiMAX transfection reagent (4.8 pmol of mimic with 0.48 µL of transfection reagent in 200 µL of proliferation medium). After 24 h of transfection, the cells were cultured under differentiation medium for 48 h. In addition, cells were treated with 100 μg/ml BrdU (Sigma Aldrich) for 1 h at day 2 of differentiation (n = 6 per group) and visualized with a rat monoclonal antibody against BrdU (ab6326; Abcam, 1:1,000), as previously described (Yoshioka et al., 2021a). BrdU-positive cells were quantified using images from six independent experiments.

Total RNAs were isolated from cells treated with mimics for the target miRNAs or negative control (n = 6 per group) using the QIAshredder and RNeasy mini extraction kit or the miRNeasy mini kit (QIAGEN), as previously described (Suzuki et al., 2019; Yan et al., 2020). In addition, total RNAs were isolated from ameloblasts at each stage of differentiation (pre-secretion, secretion, and maturation) in the lower incisors of 8-week old males C57BL/6J mice (n = 3). Briefly, the lower incisors were extracted, and ameloblasts were manually dissected and separated into three parts [apical 1/3 (pre-secretion), middle 1/3 (secretion), and incisal 1/3 (maturation) between the cervical loop and bony ridge of the incisor] under a dissection microscope. cDNA was reverse-transcribed with the iScript Reverse Transcription Super Mix (BioRad) and amplified with the iTaq Universal SYBER Green Super Mix (BioRad) using a CFX96 Touch Real-Time PCR Detection System (BioRad). The expression of genes was normalized with Gapdh. miRNA expression during ameloblast differentiation was detected with Taqman Fast Advanced Master Mix and Taqman Advanced miR cDNA Synthesis Kit (Thermo Fisher Scientific), according to the manufacturer’s instructions. The PCR primers used are listed in Supplementary Table S1.

The cells were plated onto ibiTreat 8-well μ-slides (ibidi GmbH, Munich district, Germany) at a density of 10,000/chamber and maintained until 80% confluency. The cells were then treated with mimics for miR-16-5p, miR-27b-3p, or a negative control, using Lipofectamine RNAiMAX transfection reagent (4.8 pmol of mimic with 0.48 µL of transfection reagent in 200 µL of differentiation medium) (n = 4 per group). After 24 h of treatment, the medium was replaced with differentiation medium for 2 days. AMELX expression was detected with anti-AMELX rabbit polyclonal antibody (ab153915, Abcam, 1:250), as previously described (Yoshioka et al., 2021b). Immunofluorescent images were captured with a confocal microscope (Ti-E, Nikon United States).

The cells were plated onto 12-well plates at a density of 60,000 per well, maintained until 80% confluence, and treated with either miR-16-5p, miR-27b-3p, or a negative control mimic, for 24 h (n = 3 per group). The cells were then cultured in ameloblast differentiation medium for another 48 h. The treated cells were lysed with RIPA buffer (Thermo Fisher Scientific) containing a protease inhibitor cocktail (Roche) and centrifuged at 21,130 × g for 20 min at 4°C. The protein concentration of the supernatants was measured with the BCA protein kit (Pierce). Protein samples (30 μg) were applied to Mini-PROTEAN TGX Gels (Bio-Rad) and transferred to a polyvinylidene difluoride (PVDF) membrane. Anti-AMELX rabbit polyclonal antibody (ab153915, Abcam, 1:1,000), anti-KLK4 rabbit polyclonal antibody (PA5-109888, Thermo Fisher Scientific, 1:750), anti-MMP20 rabbit polyclonal antibody (55467-1-AP, Proteintech, 1:750), and anti-GAPDH mouse monoclonal antibody (MAB374, Millipore, 1:6,000) were used for immunoblotting. Peroxidase-conjugated anti-rabbit IgG (7074, Cell Signaling Technology, 1:100,000) and anti-mouse IgG (7076, Cell Signaling Technology, 1:100,000) were used as secondary antibodies. All immunoblotting experiments were performed three times to validate the results.

Cells were plated on 12-well cell culture plates at a density of 60,000 cells per well, or on ibiTreat 8-well μ-slides (ibidi GmbH, Munich district, Germany), at a density of 10,000 cells per well and maintained until 80% confluence. The cells were treated with mimics for a negative control, miR-16-5p, or miR-27b-3p (4.8 pmol for 12-well plates and 1.2 pmol for ibiTreat 8-well μ-slides) with a combination of overexpression vectors [100 ng (12-well plates) or 25 ng (ibiTreat 8-well μ-slides)] using Lipofectamine 3000 transfection reagent (Thermo Fisher Scientific), according to the manufacturer’s protocol, which was followed by treatment with Eda (Antibodies-online Inc., ABIN3291185), Relt (Antibodies-online Inc., ABIN4054001), or Smad3 (Antibodies-online Inc., ABIN3809504) for the miR-16-5p mimic, or Bmp2 (Antibodies-online Inc., ABIN4045152), Pax9 (Antibodies-online Inc., ABIN4216431), or Slc24a4 (Addgene, 75208) for the miR-27b-3p mimic (n = 6 per group). After 24 h of transfection, the medium was switched to differentiation medium for 2 days.

Statistical comparisons between two groups were performed with a two-tailed Student’s t-test. Multiple comparisons were conducted with one-way analysis of variance with the Tukey–Kramer post hoc test. A p-value of less than 0.05 was considered as statistically significant. For all groups, data were represented as mean ± SD.

A total of 4,846 articles were extracted from a database compilation of multiple sources through a search conducted using Rayyan (Ouzzani et al., 2016). After resolving duplicates with RefWorks, 2,306 articles were selected for further screening. A total of 2,207 articles were excluded because there was no underlying genetic mechanism dictating the gene findings or the articles did not mention any relevant study or research conducted in humans. A total of 99 articles were further reviewed and qualified through a full-text review (Figure 1A), referring to 89 studies in mice, seven in rats, two in dogs, and one in cattle. A total of 44 genes [42 genes in mice with single gene mutations and two additional genes (Bmp4 and Stim2) in compound mutant models] were identified in mice as genes associated with amelogenesis imperfecta through the systematic review (Supplementary Table S2). A search of the Mouse Genome Informatics (MGI) database identified a total of 59 mouse lines after the removal of duplicates. Upon validation of the enamel phenotype through review of the extracted articles, we identified 35 genes primarily associated with amelogenesis imperfecta (Supplementary Table S3). Among these 35 genes, 15 were uniquely found in the MGI search, and 19 were common in the systematic review and MGI search. Through a manual literature search, we identified additional 11 genes associated with amelogenesis imperfecta (Supplementary Table S4). As a result, a total of 70 genes were identified and curated [68 genes in single-gene mutant mice (Table 1) and two additional genes (after exclusion of overlapping genes in Table 1) in compound mutant mice (Table 2)] as genes associated with amelogenesis imperfecta (a.k.a. enamel hypoplasia) in mice (Figure 1B), hereafter referred as mouse amelogenesis imperfecta-related genes (mAIGenes). In addition, we found that three genes in rats, three genes in dogs, and one gene in cattle were reported in amelogenesis imperfecta (Supplementary Table S5). Among the 70 genes, mutations in 33 genes were reported in humans with amelogenesis imperfecta in isolated or syndromic cases.

These mAIGenes were further categorized into three classes of amelogenesis imperfecta based on gross anatomical observation, histological analysis, microCT, and component analyses, which all are established in human cases: hypoplastic/enamel hypoplasia/no enamel (40 genes), hypomaturation (6 genes), hypomineralized/hypocalcified (39 genes), and unknown detailed classification (4 genes) (Table 3). Some genes exhibited a combined phenotype, as seen in humans. It should be noted that different mutational strategies for deletion, overexpression, or knock-in of the same gene sometimes resulted in different tooth phenotypes. This suggests that subtle changes in the expression or deletion of non-coding genomic sequences may affect the expression and function of genes that are crucial for enamel formation.

Among the mAIGenes, 12 genes (Ambn, Amelx, Amtn, Col17a1, Csf1, Dmp1, Dspp, Enam, Lama3, Lamb3, Lamc2, and Postn) were grouped in the extracellular matrix (ECM) pathway, 11 genes (Alpl, Fam20a, Fam20c, Hras, Klk4, Map3k7, Mmp20, Plau, Rac1, Rhoa, and Tcirg1) in the enzyme pathway, and seven genes (Cftr, Cnnm4, Slc4a4, Slc10a7, Slc12a2, Slc13a5, and Slc24a4) in the ion exchanger/transporter pathway. Moreover, three genes (Stim1, Stim2, and Tmbim6) were related to a calcium ion sensor or regulator, and six genes (Cldn3, Cldn16, Ctnnb1, Gja1, Nectin1, and Nectin3) were involved in cell-cell or cell-ECM adhesions. Since ameloblasts secrete enamel proteins, mutations in genes related to ECM and enamel proteins support their causal roles in amelogenesis imperfecta. In addition, a substantial number of genes were involved in growth factor signaling cascades: four were growth factors (Bmp2, Bmp4, Gdnf, and Tgfb1), five receptors (Fgfr1, Itgb1, Itgb6, Relt, and Tgfbr2), 15 transcription factors (Ascl5, Bcl11b, Ctnnb1, Dlx3, Foxo1, Irf6, Msx2, Pax9, Pitx2, Runx1, Runx2, Sp3, Sp6, Sp7, and Tbx1), two transcriptional regulators (Hmgn2 and Med1), and one a signal mediator (Smad3). Since these factors are involved in various developmental processes, the mutations would be related to syndromic cases with various developmental defects beyond amelogenesis imperfecta (Table 4).

To further explore the functional features of mAIGenes, we performed a functional enrichment analysis and functional module cluster analysis (Figure 2A). Using a false discovery rate (FDR) < 0.01, we obtained 32 gene sets that were significantly enriched in mAIGenes, including four pathways from Kyoto Encyclopedia of Genes and Genomes (KEGG) annotations, 24 Gene Ontology (GO) Biological Process (BP) terms, three GO Cellular Component (CC) terms, and one GO Molecular Function (MF) term (Table 5). Among the top 20 most significant gene sets, genes associated with tooth development (e.g., enamel mineralization, biomineral tissue development, and odontogenesis of dentin-containing tooth) were among the most significantly enriched (Figure 2B). To investigate how these functional terms and pathways are interrelated, we used the k-means algorithm to cluster them (Supplementary Figure S1A). This analysis revealed four groups (Table 5), including the 32 gene sets mentioned above and 63 mAIGenes in the module network (Figure 2C). The groups were ordered by number of gene set, with the smaller number being named first. Group 1 had one gene set— “Protein binding” —, which included Stim1, Stim2, Slc4a4, Slc12a2, etc. (Figure 2C), whereas Group 2 was related to biomineral development and ECM. These two pathways are closely related, since most of the biomineral development process occurs in extracellular fluids (Figure 2C). Enam, Ambn, Amtn, and Amelx were commonly involved in biomineralization during tooth enamel development and located at the ECM (Figure 2C). Cell proliferation and cancer-related gene sets were clustered in Group 3, including “Cell proliferation”, “Pathways in cancer”, “Cell adhesion”, and “Positive regulation of cell migration” (Figure 2C). Group 4 highly reflected the tooth and bone development, as it contained “Odontogenesis of biomineral tissue development”, “Odontogenesis of dentin-containing tooth”, “Ossification”, “Osteoblast differentiation”, and “Skeletal system development” (Figure 2C). Bmp2, Bmp4, and Runx2 connected most of the gene sets in Group 4 (Supplementary Figure S1B), and these genes have been reported to play critical roles in bone development.

For the miRNA-mAIGene regulatory network analysis, we performed miRNA-mAIGene enrichment analysis and miRNA regulatory network analysis (Figure 3A). We identified 35 Highly Expressed MiRNAs (HEMs) in mouse incisors and 32 HEMs in molars with a frequency >1%; 26 mouse tooth HEMs were then curated by taking the intersection of the incisor and molar HEMs (Supplementary Table S6) [26]. A total of 21 of these HEMs did not have a confident -3p or -5p; therefore, we considered that these had both -3p and -5p and identified 47 HEMs, all with a certain -3p or -5p. Based on these 47 HEMs, we predicted that 32 HEMs might target the 42 mAIGenes by using our pipeline and the four miRNA-target gene databases: TargetScan, miRanda, miRTarBase, and PITA. By performing the miRNA-mAIGene regulatory relationship enrichment analysis with a cutoff adjusted p-value < 0.05, we identified 27 notable miRNAs, 41 genes, and 161 miRNA-mAIGene pairs. A total of 17 miRNAs or miRNA groups, 41 genes, and 103 miRNA-mAIGene pairs were extracted after merging the miRNAs or miRNA groups that shared the same targets (such as miR-23a/b-3p and miR-125a/b-5p) (Table 6). Three miRNAs (miR-16-5p, miR-27b-3p, and miR-23a/b-3p) were considered to be hubs in the miRNA regulatory network (Figure 3B) because they had the highest degrees (Figure 3C, Supplementary Table S7), which are defined as the number of partners that immediately interact with a node of interest in the network (Sun et al., 2012), and the lowest adjusted p-values (Table 6). The sub-network of miR-16-5p, miR-27b-3p, and miR-23a/b-3p showed that Smad3 was regulated by all the three hub miRNAs. Stim2, Csf1, Slc10a7, Bcl11b, Slc12a2, Slc12a2, Slc4a4, and Pax9 were regulated by two of these three hub miRNAs or miRNA group, whereas the other genes were regulated by one miRNA or miRNA group (Figures 3D,E). As above, miR-16-5p, miR-27b-5p, and miR-23a/b-3p were considered to be promising miRNA candidates for amelogenesis imperfecta in mice.

To evaluate the function of the miRNAs predicted by the bioinformatic analyses, we conducted ameloblast differentiation assays using mHAT9d cells, a mouse dental epithelial cell line. Although the mouse ameloblast-like cells LS8 (Chen et al., 1992) have been widely used for ameloblast studies, they are limited in their ability to differentiate. We analyzed both LS8 and mHAT9d cells under differentiation conditions and found that mHAT9d cells reacted better to the induction of differentiation (Figure 4A). For instance, the expression of the ameloblast differentiation maker genes was induced more strongly in mHAT9d cells compared to LS8 cells (Supplementary Figures S2, S3). Therefore, mHAT9d cells were used in this study. We found that expression of ameloblast differentiation marker genes (i.e., Ambn, Amelx, Enam, Klk4, and Mmp20) was induced with ameloblast differentiation medium (Figure 4B, Supplementary Figure S3). In addition, we tested whether other genes associated with amelogenesis imperfecta were induced. Among the 27 genes regulated by miR-16-5p, miR-23a-3p, miR-23b-3p, miR-27b-3p, and miR-214-3p, we found 14 genes that were upregulated under differentiation conditions (Supplementary Figure S4). miR-16-5p and miR-27b-3p were induced at relatively high expression levels in mHAT9d cells, and their expression did not change under differentiation conditions (Supplementary Figure S5A). In addition, we found that miR-16-5p and miR-27b-3p were expressed at the pre-secretion, secretion, and maturation stages of ameloblast differentiation in mouse lower incisors (Supplementary Figure S5B). Overexpression of either miR-16-5p or miR-27b-3p significantly anti-correlated with downregulation of expression of Amelx and Enam, but not Ambn, Klk4, and Mmp20, in mHAT9d cells (Figure 4B). We confirmed that the expression levels of AMELX, but not KLK4 and MMP20, were decreased by overexpression of miR-16-5p and miR-27b-3p with immunoblotting (Figure 4C). The expression of AMELX was further confirmed by immunocytochemical analysis (Figure 4D). By contrast, mimics for miR-23a-3p, miR-23b-3p, and miR-214-3p did not affect the gene expression of the ameloblast differentiation makers (Figure 4B). These results indicate that miR-16-5p and miR-27b-3p may play a critical role in ameloblast differentiation through the regulation of genes that are crucial for ameloblast differentiation.

Next, to identify the miRNA-mAIGene regulatory mechanism(s), we conducted quantitative RT-PCR (qRT-PCR) analyses for the predicted target genes for each miRNA (Bcl11b, Csf1, Eda, Fgfr1, Med1, Relt, Slc4a4, Slc10a7, Slc12a2, Smad3, Stim1, and Stim2 for miR-16-5p; Bmp2, Csf1, Foxo1, Pax9, Runx1, Slc10a7, Slc24a4, Smad3, Sp6, Sp7, Stim2, and Tmbim6 for miR-27b-3p) in mHAT9d cells. The expression of Eda, Relt, Slc4a4, and Smad3 was significantly downregulated in mHAT9d cells treated with miR-16-5p mimic (Figure 5A, Supplementary Figure S6). Similarly, the expression of Bmp2, Pax9, and Slc24a4 was significantly downregulated in mHAT9d cells treated with miR-27b-3p mimic (Figure 5B, Supplementary Figure S6). Furthermore, we confirmed that treatment of inhibitor for either miR-16-5p or miR-27b-3p had no effect on expression of Amelx and Enam, while the expression of the target genes of each miRNA was upregulated (Supplementary Figure S7). Indeed, the predicted target genes contained miRNA recognition sites for their correlated miRNAs on the 3′-UTR (Supplementary Figure S8). By contrast, there was no potential recognition site for miR-16-5p on Amelx and Enam and for miR-27b-3p on Amelx, while there was a potential recognition site for miR-27b-3p on Enam, and treatment with either mimic or inhibitor for miR-16-5p and miR-27b-3p failed to alter the expression of Amelx and Enam, suggesting that these genes are indirectly regulated by miR-16-5p and miR-27b-3p in mHAT9d cells.

Finally, to examine the functional relevance of genes that were significantly downregulated under treatment with either miR-16-5p or miR-27b-3p mimic, we conducted rescue experiments by overexpressing the target genes (Figure 6A). We found that overexpression of Eda, Relt, and Smad3 under conditions of overexpression of miR-16-5p partially restored mRNA and protein expression of Amelx and Enam (Figures 6B,C). Similarly, overexpression of Bmp2, Pax9, and Slc24a4 partially restored mRNA and protein expression of Amelx and Enam when miR-27b-3p was overexpressed (Figures 6B,C). Taken together, our results show that overexpression of miR-16-5p and miR-27b-3p inhibits ameloblast differentiation through the regulation of mAIGenes.