All animal care and use procedures were approved by the Animal Care and Use Committee of Xinjiang University (IACUC-20210709) and performed according to the guidelines of the U.S. National Institutes of Health. Female CD-1 mice (6–8 weeks old) and male CD-1 mice (8–9 weeks old) were purchased from the Animal Resource Centre of Xinjiang Medical University and housed in a specific pathogen-free facility of Xinjiang University under a standard 12-h light/12-h dark cycle with free access to food and water.

Reverse transcription reagents and qPCR reagents were purchased from Applied Biological Materials Inc (abm, Richmond, Canada). Unless otherwise specified, other reagents were purchased from Sigma-Aldrich (MO, United States).

Female mice were treated with pregnant mare serum gonadotropin (PMSG) and human chorionic gonadotropin (high) to induce superovulation and were mated with male mice as previously described (Wu et al., 2013). With the observation of a copulatory plug the next morning, the pregnancy stage was considered 0.5 days post-coitus (dpc). The uterus with a short Fallopian tube was recovered at 3.5 dpc. Subsequently, a blunt 30-gauge needle was inserted into the Fallopian tube, and expanding blastocysts were flushed from the uterus using 200 μL of M2 medium. Embryos were collected and cultured in KSOM medium under mineral oil. After culture for 6–8 h, blastocyst hatching was classified based on the hatching site (A-site, B-site, or C-site) (An et al., 2021). After 16 h of culture, blastocysts were divided into hatched and hatching failure groups (hatching embryos, H; non-hatching embryos, N; Figure 1A).

Total RNA extracted from embryos using TRIzol method. Expanding blastocysts (E), blastocysts with A-site, B-site, or C-site hatching sites (A, B, C), hatched blastocysts (H), and non-hatching embryos (N) were analyzed by RNA sequencing (Figure 1A). In each group, a panel of 30 embryos was collected and stored in 100 μL of TRIzol (Thermo Fisher Scientific, MA, United States) with three replicates (total 90 embryos). Transcriptome sequencing with Smart-Seq (Goetz and Trimarchi, 2012) was conducted by the Guangzhou GENE DENOVO Company. Differentially expressed genes (DEGs) were analyzed by EdgeR and normalized in fragments per kilobase of transcript per million mapped reads (FPKM, see Supplementary Table 1). DEGs and gene sets were analyzed using Gene Ontology (GO) term enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis using the online analysis system (Guangzhou GENE DENOVO Company, Guangzhou, China). In transcription factor-target gene regulatory network analysis, transcription factors (TFs) were identified in the selected gene set, the JASPAR database was used to obtain transcription factor binding motifs, and MEME FIMO software was used to predict the transcription target genes in the selected set. A transcription factor targeting network showed the association between TFs and target genes.

A pool of 10 embryos in each group was collected, and cDNA was synthesized by All-In-One 5 × RT MasterMix according to the manufacturer’s instructions (abm). For single embryo measurements, we lysed one embryo for cDNA synthesis using a scaled-down reaction based on the manufacturer’s instructions (abm). Simply, cDNA synthesis was modified by treating an embryo in a micro-drop (5 μL) scaled-down reaction solution on a plastic dish covered with mineral oil to prevent evaporation of fluids. The cDNA was stored at −20°C for qPCR. Through transcriptome analysis we identified ten key genes (Cd36, Ccl9, C3, Cyp17al, Il-1a, Ccl5, Susd4, Cfb, Ptgs1, and Lyz2) for analysis and designed qPCR primers for them (primer sequence see Supplementary Table 2). Those 10 key genes selected from GO and KEGG enrichment in immunity and with relatively high expression (FPKM>1). The housekeeping gene Gapdh was used for RNA normalization. The 20-μL qPCR reaction comprised 10 μL of Blasaq 2 × qPCR Master Mix (abm), 0.5 μL of each upstream and downstream primer, and 2 μL of cDNA. The products were electrophoresed on a 2% agarose gel to confirm amplification. For RT-qPCR analysis, data were obtained from three technical replicates in each experiment and three biological replicates in each blastocyst group.

Expression of the protein products of C3, Cdx2, Il-1β, Lyz2, and Plac1 was determined using immunofluorescence (IF) as previously described (An et al., 2021). Briefly, the embryos were washed three times in 0.1% polyvinyl alcohol (PVA) for 3–5 min, fixed in 4% paraformaldehyde buffer for 20 min, followed by three washes with 0.1% PVA. Embryos were permeabilized with 0.5% Triton X-100 at room temperature for 15 min, washed three times with 0.1% PVA, and blocked in 4% fetal bovine serum at room temperature for 1 h. Embryos were incubated with the primary antibody C3, Cdx2, IL-1β, Lyz2 and Plac1 respectively (diluted 1:200 in 2% phosphate-buffered saline, Abcam, Cambridge, United Kingdom), at 37°C for 1 h or at 4°C overnight. Embryos were washed with 0.1% PVA and incubated in the corresponding secondary antibody (goat anti-rabbit IgG AF488, Abcam) in the dark at 37°C for 1 h, washed three times with 0.1% PVA, and stained with 4′,6-diamidino-2-phenylindole (DAPI) for 5 min. Images were acquired using a laser-scanning confocal microscope (Nikon, Japan).

Using the R language survival and glmnet package, we prepared the input data and then built the LASSO regression model. We selected the appropriate λ value and remodeled the data using this λ value. We performed cross-validation on the LASSO analysis results to determine the gene data to be included in the model, ensuring a rational and effective model.

Blastocysts were classified into five groups according to their hatching status, including expanding blastocysts (E), hatching blastocysts with A-site (1-2 o’clock), B-site (3 o’clock), C-site (4-5 o’clock) hatching sites (A, B, C), hatched blastocysts (H) and non-hatching embryos (N). Blastocysts hatching at O-site (12 o’clock) and D-site (6 o’clock) were excluded due to the rarity of hatching events and the difficulty of collecting enough embryos for analysis. First, gene expression profiles were analyzed to describe the general state of development during blastocyst hatching. Then, to explain how blastocyst hatching determines its implantation, a comparative analysis was performed in blastocysts hatching from the B-site (good pregnancy outcome) vs. the C-site (poor pregnancy outcome). To reveal how hatching occurs and is achieved, blastocysts were analyzed at the stages from expansion to hatching and hatched. To reveal the intrinsic determinants of different hatching outcomes, hatched and non-hatched blastocysts were analyzed comparatively. Based on the above data analysis, a predictive model of implantation success was developed. A modified single blastocyst gene expression detection approach was established to confirm the expression profiles of key genes in the prediction model.

The data obtained for RT-qPCR and IF are presented as mean ± standard error of the mean (SEM). One-way analysis of variance with a post hoc two-sided Sidak t-test was used to determine differences between treatment groups. All statistical analyses were performed using GraphPad Prism 9.01 for Windows, version 24.0 software (GraphPad Software, La Jolla, CA, United States). A p < 0.05 difference between groups was considered significant.

We previously found a strong association between the birth rate and the hatching states in mice (An et al., 2021). To further investigate this, we performed Smart-Seq on pre-hatching or expanding blastocysts at 3.5 dpc (E); blastocysts hatching from A-site (A), B-site (B), and C-site (C) at 3.5 dpc plus 8 h; and hatched blastocysts (H) as well as non-hatching blastocysts (N) at 3.5 dpc plus 16 h (Figure 1A),consistent with research groups as previous study (An et al., 2021). Global gene expression, analyzed by principal component analysis (PCA), showed distinct groups based on hatching states (Figure 1B). A and B, which resulted in higher embryo implantation (55.6%, 65.6%) (An et al., 2021), clustered closely at the left side of the plot, whereas C and N, with very low birth rates of 21% and 5.2%, respectively, clustered closely. In addition, fully hatched blastocysts formed a distinct cluster. Hierarchical cluster analysis of gene expression segregated the developing blastocysts into a cluster that comprised expanding blastocysts and hatching A, B and C blastocysts. The other cluster comprised H and N blastocysts (Figure 1C). Thus, blastocysts with different pregnancy outcomes had marked differences in gene expression profiles. Moreover, as a blastocyst developed from expanding to hatching to hatched, its gene expression profile changed dramatically (Figure 1D), reflecting blastocyst development during this process.

Comparing the blastocysts hatching from the B-site (65.6% birth rate) vs. the C-site (21.3% birth rate) (An et al., 2021), we found that 75 genes were upregulated, and 103 genes were downregulated (Figure 2A). Among DEGs, 1) biology process of inflammatory response and immune response, and 2) pathway of infectious disease were highly enriched by GO and KEGG analysis (Figures 2B, C). As immune genes Cd36, Cfb, Ccl9, Il-1β, Ccl5, Il-1α, Ptgs1, Lyz2, C3, and Susd4 showed significant differences between B-site and C-site blastocysts (Figure 2D), immunological changes in the hatching blastocyst likely allow it to recognize the endometrium at embryo implantation.

For B-site and C-site hatching blastocysts, DEGs fell into groups with > 5-fold change or < 5-fold change (Figure 2A), suggesting that genes with a > 5-fold change may be transcriptionally coregulated during the 8 h of blastocyst hatching. Indeed, we found by transcription factor-target gene regulatory network analysis that 2 TFs, TCF24 and DLX3, controlled over 82 out of 178 DEGs (46%) (Figure 2E). We selected four DEGs, Lyz2, Ptgs1, Il-1β, and Susd4, and two previously investigated genes significantly correlating with implantation, Lifr and Upa (An et al., 2021), and analyzed the correlation of their expression level with birth rate. Lyz2, Ptgs1, and Il-1β were positively correlated with birth rate, whereas Lifr, Upa, and Susd4 were negatively correlated with birth rate (Figures 2F, G). Differentiation genes coordinate and determine developmental fate during blastocyst hatching. Therefore, these gene products may be involved in maternal-fetal interactions during implantation, which helps to determine pregnancy outcomes.

We used RT-qPCR to measure transcription of immune-related genes Ptgs1, Lyz2, Il-α, Cfb, Ccl9, Cd36, Ccl5, C3, Cyp17 al, and Susd4 in five blastocysts each from expanding, hatching, hatched, and non-hatching groups. We found upregulation of Lyz2, Il-α, Cfb, Ccl9 and Susd4, and downregulation of Cd36 expression during blastocyst hatching (Figure 3A). Ccl5, C3, Cyp17al, and Ptgs1 exhibited low expression levels or were undetectable. There was a different level of gene expression in the non-hatching blastocyst. The RT-qPCR gene expression patterns were consistent with the transcriptome analysis.

We used immunofluorescence (IF) microscopy to measure genes Cdx2 and Plac1 (regulating trophoblast differentiation), C3, IL-1β, and Lyz2 (immune-relative genes) in each hatching state (Figure 3). IF signals for Plac1, C3, IL-1β, and Lyz2 were mainly in the cytoplasm of trophoblast cells, whereas Cdx2 was expressed in the nuclei of trophoblast cells (Figure 3). Notably, C3 and IL-1β were expressed on the extra-luminal surface of the TE of the hatched blastocysts, implying that these proteins play a role in maternal-fetal interactions when the embryo comes into contact with the endometrium. However, in non-hatched blastocysts, this protein translocation was not detected. Based on mean fluorescence intensity, we found a decrease in the expression of Plac1, Cdx2, and C3 and an increase in the expression of Lyz2 and IL-1β from expanding to hatching blastocysts (Figure 3B) consistent with the transcriptomic and RT-qPCR results. Non-hatching blastocysts had lower levels of Lyz2 and IL-1β compared to hatched blastocysts, showing a dysregulation of these proteins.

We analyzed the DEG expression patterns during the hatching process for E, B, and C blastocysts by Venn diagrams (Figures 4A, B) and found that the expression of 307 genes correlated with blastocyst hatching, thereby determining blastocyst development from expanding to hatching (A-site, B-site, and C-site). By trend analysis, the 307 DEGs showed four significant gene expression patterns, with 141 genes in profiles 6 and 7 showing an upregulated trend and 133 genes in profiles 0 and 1 showing a downregulated trend (Figure 4C). Thus, there were two dominant sets of DEGs that were associated with blastocyst development during hatching. By transcription factor-target gene regulatory network analysis, we found that 1) TF ATOH8 regulated 41.8% of DEGs with an upregulated trend (59/141, Figure 4D), 2) TFs SPIC, GM9044, GM9046 and GM9040 regulated 97% of DEGs with a downregulated trend (129/133, Figure 4E), and 3) SPIC was a major regulator of coregulated Gm9044, Gm9046 and Gm9040.

By GO analysis of DEGs, the top 10 enriched terms included the regulation of cell differentiation, cell proliferation, and developmental and apoptotic processes (Figure 4F). In the KEGG analysis, the top 10 enriched signaling pathways included cytokine-cytokine receptor interaction, viral infection, and a Toll-like receptor signaling pathway (Figure 4G). Thus, remodeling of cell processes and the shaping of immunological properties were the major developmental events during blastocyst hatching.

PCA analysis of DEGs for E, H, and N blastocysts revealed clear clusters for E, H, and N blastocysts, showing a low correlation between hatched and non-hatched blastocysts (Figure 5A). In non-hatched blastocysts, 1,433 genes were upregulated, and 743 genes were downregulated compared with hatched blastocysts (Figure 5B). Venn diagram analysis of DEGs for E, N, and H blastocysts showed that 3,428 genes correlated with differential hatching outcomes (Figure 5C), thereby determining whether blastocysts hatched. By trend analysis, we determined that 1) 868 genes showed a decreasing trend, 2) 1,217 genes showed a V-shape trend, and 3) 821 genes showed an inverted V-shaped trend (Figure 5D). For these three trends, transcription factor-target gene regulatory network analysis identified 17 TFs (852 target genes), 8 TFs (472 target genes) and 10 TFs (668 target genes) (Figure 5E), indicating complex patterns of transcriptional gene regulation that determine whether the blastocyst will hatch. Changes in the expression of these genes, which were enriched in processes related to cell surface receptor signaling, binding, and the plasma membrane (Figure 5F), affect the hatching outcome. KEGG analysis showed that the top enriched signaling pathways were related to cancer, phenylalanine metabolism, and the cAMP signaling pathway (Figure 5G). These GO and KEGG results demonstrate that during hatching, these molecular events and their transcriptional regulation control the processes that lead to blastocyst development.

We used RT-qPCR to measure the transcription of genes associated with these trends, including Podxl, Lama1, Ptk2b, Pramel17, Pramel31, and Amigo2. Lama1 and Ptk2b were upregulated, whereas expression of Pramel17, Amigo2, Pramel31, and Podxl was low or undetectable in blastocyst hatching (Figure 5I), consistent with their transcriptomic expression patterns (Figure 5H).

We ranked the DEGs associated with blastocyst development from the transcriptome analysis and identified 13 candidate genes to develop a model to predict blastocyst implantation success. The nine DEGs that were downregulated in three out of four measurements and the four DEGs that were upregulated in two out of three measurements (Figures 6A, B) are shown in a heatmap (Figure 6C). Predictive aggregate effects were identified using LASSO regression (Figure 6D). To determine the gene variables to include in the model and to ensure a rational, effective model, we used LASSO regression analysis using the R language survival and glmnet package, followed by cross-validation of the LASSO analysis results (Figures 6D–F). We screened Lyz2, Cd36, Cfb, and Cyp17a1 as variables using LASSO regression and determined by LASSO Cox regression analysis that these four genes, the minimum value for multivariate Cox modeling, were statistically significant. To determine the predictive performance of the model, we evaluated these genes using receiver operating characteristic (ROC) curve analysis, which showed that they made a significant contribution to the prediction of the model, with area under the ROC curve values generally greater than 0.8 (Figure 6G).

Using these validated genes as the variables, we plotted a nomogram graph based on the LASSO regression and cross-validation data (Figure 6H) to produce a predictive model of implantation success for hatching blastocysts. In this predictive model, the level of expression corresponded to the score for each gene. For example, if the relative expression level of Lyz2 was 1, this gene was given a score of 32. The scores for each gene were added together to produce a total score, which corresponded to the predicted value, which indicated the potential implantation outcome for a blastocyst.

We measured the expression of the genes in the prediction model in single blastocysts by RT-qPCR. For a single blastocyst hatching at the B- and C-sites, we found that the expression of Lyz2 and Cfb was consistent with data obtained from pooled blastocyst RT-qPCR, confirming the validity of the single-blastocyst technique (Figures 7A, B). Differences in the expression of Lyz2 were observed across 27 populations of expanding blastocysts, suggesting that the expression of Lyz2 is regulated by the intrinsic programming fate before hatching (Figure 7C). PCA analysis of Lyz2 expression profiles in blastocyst populations from the expanding, hatching, hatched, and non-hatched stages showed strong clustering for blastocysts in the expanding state, weaker clustering for blastocysts at the hatching state, and poor clustering for cells at the hatched stage (Figure 7D). However, based on Lyz2 expression, blastocysts that failed to hatch were dispersed (Figure 7D).