Human umbilical cord blood was obtained from the donors with normal term delivery after signing informed consent. Cord blood sample collection was approved by the Ethics Committee of Guangzhou First People’s Hospital (Approval No K-2021-008-02) and complies with the Declaration of Helsinki.

hESCs, H9 line cells (WA09), were obtained from the WiCell Research Institute (Madison, WI, USA) with the Material Transfer Agreements (No.: 19-W0512, 24-W0162 and 24-W0163), and cultured and passaged on the feeder cells of mouse embryonic fibroblasts (MEF). The culture medium was DMEM-F12 (Gibco) basal medium supplemented with 20% knockout ™ Gibco), 4 ng/ml fibroblast growth factor bFGF (Peprotech), 100 mM L-glutamine (Gibco), 1% MEM non-essential amino acids solution (Thermo), 0.0007% β Mercaptoethanol. The medium was changed daily.

H9 cells were cultured on gelatin-coated plates with a feeder layer of MEF for 5 days, then the cells were subjected for endodermal differentiation. At the first day, differentiation medium was X-VIVO 20 medium (Lonza) supplemented with 100 ng/ml Activin A (Peprotech), 50 ng/ml WNT3A (Peprotech) or 3 μM CHIR99021(Selleck), 10 μM ROCK inhibitor Y-27632 (Selleck). The following day, the medium was changed to X-VIVO 20 (Lonza) medium containing 100 ng/ml Activin A, 0.5% B27 (Gibco). At day 4 after the differentiation, endodermal cells were dissociated with TrypLE (Gibco), and reseeded on Matrigel (Corning)-coated plates with DMEM-F12 medium (Gibco)or EGM2 medium (Lonza), supplemented with 0.5% B27, 1 μM RA (Sigma), 10 μM SB431542 (MCE), 20 ng/ml BMP4 (Peprotech), 2.5μM IWP2 (Selleck) or 2.5μM IWR1 (Selleck). The medium was changed daily, and the differentiation lasted until day 9 for the induction of 3PPE.

For the differentiation of TECs, the differentiation from 3PPE cells were carried out from day 9 after the differentiation of 3PP cells employing DMEM-F12 medium supplemented with WNT3A (50 ng/mL), BMP4 (10ng/mL), FGF8 (50 ng/mL), FGF10 (5ng/mL) (all from Peprotech), Heparin sodium (10 ng/mL, GlpBio), SB431542 (10μM), SANT-1 (2.5 μmol/L) (selleck), 2-Phosphate-L-ascorbic acid trisodium salt (50 μg/mL, Sigma), B27 (0.5%, Gibco), RA (1μM, Sigma) for a total of 4 days, and the medium were changed once a day.

Total RNAs were extracted by the Trizol (Takara) from hESCs that were differentiated at a specified time using a specified culture method. RNA concentration was measured using NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific). According to the manufacturer’s instructions, 1000ng of RNAs were reversely transcribed into cDNAs using HiScript III RTSuperMix for qPCR (+gDNA wiper) kit (Vazyme). qRT-PCR was performed by using ChamQ Universal SYBR qPCR Master Mix (Vazyme) on QuantStudio (Applied Biosystems Life Technologies, Thermo Fisher Scientific). With GAPDH as internal parameter, the relative mRNA expression of samples was calculated by the 2^-ΔΔCt method. The primers used were listed in Supplementary Table 3 .

The cells were inoculated and cultured in 12 well plates with cell slides which were placed before the inoculation for 1 days, then cells were fixed with 4% paraformaldehyde (PFA, biosharp) at room temperature for at least 15 minutes, and permeated with 0.5% Triton X-100 (Solarbio) for 20 minutes. Next the cells were blocked in goat serum (Boster) for 1 hour at room temperature, then incubate with the primary antibody at 4°C overnight. Next day, after rewarming at 37°C for 1 hour and washing, the cells were incubated with the secondary antibody at room temperature in dark for 2 hours. The cells were added with DAPI (Beyotime) to stain the nucleus and sealed with a fluorescent antiaging quenching agent (Beyotime). The images were collected by Zeiss 880 multiphoton laser scanning microscope or single photon confocal microscope (Nikon). All antibodies used were listed in Supplementary Table 4 .

Cells were dissociated into single cells using TrypLE ™ Express Enzyme (1X) (Gibco), and the single cell suspension was collected and washed with PBS once for immunostaining. For the staining of the cell surface marker, the cells were directly stained and incubated at 4°C in dark for 30 minutes, and the dead cells were eliminated by DAPI or 7AAD (BD) staining. For the detection of nuclear transcription factors, FOXP3/Transcription factor staining buffer set (invitrogen) was used to fix and perforate cells. After cells were kept away from light at room temperature for 60 minutes, 1 x perm buffer was directly added to each tube to stop the reaction. The cells were then incubated at room temperature for 45 minutes with the primary antibody, afterwards PBS (Solarbio) was filled, the stained cells were centrifuged with 300g for 5 minutes, and with 8500g for 2 minutes after mixing. After cells were resuspended and incubated at room temperature for 45 minutes with the second-fluorescence antibody, then PBS (Solarbio) was filled, the stained cells were centrifuged with 300g for 5 minutes, and with 8500g for 2 minutes after mixing. Flow cytometry was performed with BD Celesta and BD Fortessa. Folwjo (V10) software was used for data analysis.

Cells were lysed by ultrasound in RIPA cell lysate buffer containing protease inhibitor and PMSF. The protein was harvested at 4°C by the centrifugation with 14000g for 15 minutes. According to the manufacturer’s protocol, BCA protein assay kit (Boster) was used for the determination of protein concentration. After protein quantification, protein denaturation was carried out and 20 µg proteins were added to 10% SDS-PAGE separation gel, and then transferred to PVDF membrane (Merck). The protein strips were sealed with 5% defatted milk powder (Solarbio) at room temperature for 1 hour, then incubated with the primary antibody, and incubated at 4°C overnight. Subsequently, the protein strips were washed by TBST (Solarbio) and incubated the second antibody at room temperature. One hour later, the expression of the target proteins was detected by the chemiluminescent ECL kit (Affinity).

Small interfering RNAs (siRNAs) against HOXA3 and EPHB2 were synthesized by Guangzhou Ige BIOTECHNOLOGY Co., Ltd, and the siRNA sequences for the negative control, HOXA3 and EPHB2 were listed in Supplementary Table 3 . At day 7 after the differentiation, when the cell density reached 70% – 80%, the differentiated cells were transduced with HOXA3 siRNA or EPHB2 siRNA and negative control siRNA using Lipofectamine 3000 (Thermo Fisher Scientific), respectively. Six to eight hours after the transduction, the differentiation medium was replaced with opti-MEM medium (Thermo Fisher Scientific), and cells were harvested at days 9 after the differentiation. siRNAs were introduced into cells for transient transduction following the manufacturer’s instructions.

The effects of HOXA3 and EPHB2 on cell proliferation were determined by a cell counting kit-8 assay (CCK-8; Dojindo). There were three groups including blank control group, NC group, siRNA-HOXA3 group or siRNA-EPHB2 group. Briefly, 5 × 10³ cells per well (in quadruplicate per group) were seeded into 96-well plates and cultured in differentiation medium. At 24, 48 and 72 hours after the culture, 10 μL of CCK-8 solution and 100 μL of fresh culture medium were added to each well and incubated at 37°C for additional 4 hours. The absorbance was observed at 450 nm by using a microplate reader (Amershan Iamger 600). The ratio of the survival/proliferation of cells was represented through the absorbance of the test wells minus the optical density of the blank wells.

For the scratch wound assay, 2 × 10⁵ cells per well (in triplicate per group) were plated into 12-well plates and incubated to reach confluence. The monolayer was scratched using a tip and washed with serum-free medium to remove detached cells. Cells were photographed at 0 and 24 hours post wounding. The closure area of wound was calculated by which migration area (%) = (A₀ - An)/A₀ × 100, where A₀ represented the area of initial wound area, An represented the area of wound at the metering point.

The migration assays were performed with the transwells, briefly 5 × 10⁴ cells per well were suspended in serum free medium and seeded into the upper chambers of 24-well transwell plates (Corning) with 8 μm pore filters, then the lower chambers were added with differentiation medium. At 48 and 72 hours after the culture, the cells attached on the upper surface of the filter membranes were removed and migrated cells of the lower surface were stained with 0.5% crystal violet for several minutes. The number of migrated cells was observed under an optical microscope (Nikon).

Cell cycle distribution was determined by a cell cycle and apoptosis analysis kit (Meilunbio) according to manufacturer’s protocols. In brief, 3PPE cells were transduced with HOXA3 siRNA or EPHB2 siRNA and negative control siRNA for 48 hours, then these cells were fixed with 70% cold ethanol at 4°C overnight. The next day, the cells were washed with cold PBS and stained with the mixture of propidium iodide and RNase A. Afterwards, the cell cycle was detected by flow cytometry (BD Biosciences).

RNAs were extracted from cells of NC group or siRNA-HOXA3 group with triplicate for each group. 100ng of total RNAs was used to prepare sequencing libraries using the Ligation Mediated Sequencing (LM-Seq) protocol (21) and quantified with the Qubit 2.0 fluorometer (ThermoFisher). Final cDNA libraries were quantitated with the Agilent 2100 bioanalyzer, multiplexed and loaded at a final concentration of 2.5 nM, and sequenced as single reads on the Illumina NovaSeq 6000 (Illumina).

DESeq2 software was used to identify differentially expressed genes (DEGs). Empirical Bayesian methods are to estimate prior values of log-fold changes and deviations and to calculate posterior values of these statistics. DEGs were evaluated by meeting the two criteria, one included more than onefold change in expression and the other included pvalue <.05. For functional enrichment analysis, all DEGs were mapped to terms in the GO databases, and then significantly enriched GO terms were searched for among the DEGs using p <.05 as the threshold. GO term analysis was classified into three subgroups, namely biological process (BP), cellular component (CC) and molecular function (MF).

GSEA (http://software.broadinstitute.org/gsea/index.jsp) was performed as previously described[26]. Analyzed gene sets (GOBP_NEGATIVE_REGULATION_OF_EPITHELIAL_CELL_MIGRATION, GOBP_POSITIVE_REGULATION_OF_EPITHELIAL_CELL_PROLIFERATION and WP_WNT_SIGNALING_PATHWAY_NETPATH) were downloaded from GSEA. To assess the enrichment, the absolute value of the GSEA metric (Pearson correlation) was considered for deregulated genes.

The ChIP was achieved using ChIP Assay Kit (Chromatin Immunoprecipitation (ChIP) Assay Kit) (Beyotime). In brief, DNA and proteins were cross-linked with 37% formaldehyde at 37°C for 10 minutes. The 3PPE cells were added into glycine solution for 5 minutes at room temperature, then the cells were washed and collected using cold PBS solution supplemented with 1mM PMSF (Solarbio). The cell precipitates obtained after the centrifugation were resuspended by SDS lysis buffer containing 1mM PMSF, and fully lysed on ice for 10minutes. Next, cell lysates were dealt with pulses ultrasonication to break nuclear membrane and shear DNA into fragments with 200–1000 bp fragments. The final lysates were incubated with 1 μg HOXA3 antibody (Santa Cruz) or 1 μg mouse immunoglobulin G (IgG) antibody (Beyotime) at 4°C overnight. All told, 10 μg chromatin was used as input control and mouse anti-IgG was used as negative control. The precipitated DNAs were subjected to PCR analysis using a primer pair specific for EPHB2. The amplified fragments were analyzed by 1.5% (w/v) agarose gel analysis and verified by qRT-PCR.

Raw data for single-cell sequencing provided by Margaret E. Magaletta et al (22). The scRNA-seq short reads and count matrices data were used in this study from the Gene Expression Omnibus (GEO) database under accession codes “GSE182135”. The atlas had a total of 57717 cells and ten samples (E9.5:2; E10.5:3; E11.5:3; E12.5:2). A typical Seurat workflow was followed, using gene selection (with FindVariableGenes), scaling (with ScaleData), and principal components analysis (PCA with RunPCA). For this initial analysis, 50 principal components (PCs) were retained for downstream steps, and the Louvain graph clustering (via FindClusters) was used with a resolution of 2 graph-based clustering with the Louvain algorithm with FindClusters, and UMAP visualization with RunUMAP. Out of the resulting 47 clusters, we removed cluster 12, 23, 27, 29, 31, 32, 36, 38, 39, 41, 42 and 44 due to low Pax9 and Epcam (associated with neuronal signature). The 49,087 cells passing quality control were analyzed using Seurat to produce overview figures. 60 principal components (PCs) were retained for downstream steps, and the Louvain graph clustering (via FindClusters) was used with a resolution of 0.5. The code and documentation are available at the GitHub repositories [https://github.com/fyj3600/HOXA3/tree/fyj3600-patch-1].

All statistical analyses were performed using the SPSS 13.0 package (SPSS International, Chicago, IL, USA). Continuous variables were presented as mean ± SD and analyzed using the Student’s t-test or ANOVA (analysis of variance).

To investigate the molecular processes underlying the development of pharyngeal endoderm, we analyzed a single-cell transcriptomic catalog (22) of mouse pharyngeal endoderm between embryonic days E9.5 and E12.5 which covers the transition from pharyngeal endoderm to pharyngeal organ primordia (23). To specifically isolate the pharyngeal endoderm, Epcam/Pax9^VENUS cells were sorted as previously described (22). After the initial data processing and quality control of the single-cell RNA (scRNA) profiles (see “Methods”), the resulting dataset yielded 49,087 single-cell transcriptomes (9,626 from E9.5, 11,230 from E10.5, 17,155 from E11.5 and 11,076 from E12.5). First, the single cell transcriptome data of mouse pharyngeal endoderm between embryonic days E9.5 and E12.5 were dimensional reduction and unsupervised clustering ( Figure 1A ; Supplementary Figure 1A ). We performed UMAP display of the pharyngeal endoderm specific genes Pax9 and Epcam, and the results showed that almost all cells expressed Pax9 and Epcam ( Figure 1B ). Importantly, known markers were used to identify 3PP (24), including Bmp4, Six1, Pax1, Hoxa3 and Gcm2 genes. We identified 27 distinct clusters ( Supplementary Figure 1A ), then we analyzed the expression levels of related markers of 3PPE in each cluster using violin plots ( Supplementary Figure 1B ). 3PPE cluster was identified by manual inspection of differentially expressed genes for Hoxa3, Bmp4, Foxn1 and Gcm2. Finally, 8, 13 and 14 clusters were defined as 3PPE cluster ( Figures 1C, D ; Supplementary Figure 1B ). The results indicated that Bmp4 and Six1 were not specific to 3PP, Pax1 was also not restricted to 3PP, and Gcm2 was expressed in the parathyroid domain of the 3PP or parathyroid, although all these markers were expressed in 3PP ( Figure 1D ). In parallelly, we also characterized key genes of the thymic epithelium, and the results showed that Cldn3, Ly75, Krt8 were more widely distributed in the pharyngeal pouches, but Foxn1, Psmb11 and Krt5 were specifically highly expressed at day E12.5 in mouse embryos ( Figure 1E ). We subsequently performed further analysis of these genes for Hoxa3, Eya1, Six1, Epcam as well as Pax9. It was worth noting that these genes presented a dynamic change pattern ( Figure 1F ). In particular, the expression of Hoxa3 in 3PPE cell cluster was increased rapidly to the highest level at day E10.5 in mouse embryos, followed by a rapid decrease ( Figure 1F ). At the same time, the dynamic changes in the expression levels of Pax9 and Six1 as well as Eya1 were also similar to those of Hoxa3 ( Figure 1F ). These results indicated that the key transcription factors of the third pharyngeal pouch showed a dynamic change trend of initial increase and later decrease in a very short period of time during mouse embryonic development.

To investigate the role of HOXA3 during the development of 3PPE, we established an approach to differentiate 3PPE from hESCs. We developed four different methods to differentiate hESCs and thus explored the most efficient one for the differentiation for 3PPE ( Figure 2A ). During the differentiation of definitive endoderm (DE), we used the classical DE differentiation approach with different combinations of Activin A (19), WNT3A and Activin A (11), or CHIR99021 and Activin A (25). At day 4 after the differentiation, the cells were examined by RT-qPCR, and the results showed that mesodermal marker Brachyury (BRA) was significantly downregulated ( Figure 2B ), and DE markers CXCR4, SOX17, FOXA2, and GATA4, were significantly upregulated in cells treated with the combination with WNT3A and Activin A when compared to those with combination with CHIR99021 and Activin A ( Figure 2C ). Thus, these results indicated that the combination with WNT3A and Activin A was more oriented toward definitive endodermal commitment.

Immediately afterwards, we continued the induction of endodermal cells toward 3PPE using EGM2 medium or DMEM-F12 medium supplemented with 0.5% B27, 1 μM RA, 10 μM TGF-β Inhibitor SB431542, 20ng/ml BMP4, 2.5 μM IWP2 or 2.5 μM IWR1. During the induction, the morphology of cells induced by the combination of DMEM-F12 and WNT3A was more homogeneous ( Figure 2D ), while cells with multiple morphologies were produced using EGM2 medium. Whereupon the cells were assayed at day 9 after the differentiation by qRT-PCR, and results showed the highest expression levels of key transcription factors TBX1, SIX1, PAX8 and PBX1 were obtained in 3PPE under the condition with DMEM-F12 and WNT3A ( Figure 2E ). Meanwhile, flow cytometry analysis also exhibited that the culture condition with the combination of DMEM-F12 and WNT3A produced 95% of SIX1⁺EpCAM⁺ cells at day 9 after the differentiation, while only 17% of the same cells were produced with the combination of EGM2 and WNT3A ( Figure 2F ; Supplementary Figure 2A ). Additionally, the similar proportion of PAX9⁺EpCAM⁺ cells produced by the combination of DMEM-F12 and WNT3A was also much higher than that of the combination of EGM2 and WNT3A ( Figure 2G ; Supplementary Figure 2B ). Therefore, an efficient approach for directed differentiation into 3PPE from hESCs was eventually determined.

The successful differentiation of the endoderm lays a tamped-down foundation for the formation of 3PPE. Under the microscope, we observed the uniform morphology of endodermal cells, which had clear cell borders, gray brown cytoplasm, a uniform shape in size and were arranged in a single layer ( Figure 3A ). To further characterize the endodermal cells, we collected cells at days 2 and 4 for analysis, and found that positive cell proportion for the pluripotent OCT4 gene was dramatically decreased ( Supplementary Figure 3A ), and the expression levels of endodermal genes CXCR4, GATA4, SOX17, CER1, FOXA2 all were gradually increased over time after the differentiation ( Figure 3B ), whereas the expression level of mesodermal gene BRA was first up-regulated then down-regulated ( Figure 3C ). In addition, the proportions of endodermal SOX17⁺ cells and CXCR4⁺ cells were increased from 16.1% to 18.0% ( Supplementary Figure 3B ), and 62.2% to 92.0% ( Supplementary Figure 3C ) respectively, and the proportion of CXCR4⁺SOX17⁺ double positive cells were increased from 14.7% to 17.8% ( Figure 3D ; Supplementary Figure 3D ). Moreover, co-immunofluorescent staining analysis further verified that the endodermal cells expressed the endodermal marker FOXA2, but did not express the pluripotent gene SSEA4 at day 4 after the differentiation ( Figure 3E ).

Next, the endodermal cells were further induced toward 3PPE cells, and qRT-PCR was performed to determine the gene expressions of cells at days 0, 4, 7, and 9 during the induction. As expected, the expression levels of 3PPE marker genes EYA1, PAX9, HOXA3, and SIX1 were all highly expressed over time during the entire differentiation from day 0 to day 9 ( Figure 3F ). Moreover, the positive cell proportions of SIX1, a key marker of 3PPE, and EpCAM, a surface marker of epithelial cells, were determined by flow cytometry, and the results showed that SIX1⁺ cell proportions were increased from 31.4% to 96.6% at day 9 ( Supplementary Figure 4A ), and the proportions of SIX1⁺EpCAM⁺ double positive cells were increased from 30.4% to 95.1% at day 9 ( Figure 3G ; Supplementary Figure 4C ), indicating that pure population of 3PPE were generated. Meanwhile, these results were also verified by co-immunofluorescent staining for SIX1 and EpCAM ( Figure 3H ). Next, we also found that the proportions of cells positive for PAX9, another key marker of 3PPE, was increased from 40.0% to 95.3% ( Supplementary Figure 4B ), while the proportions of double positive cells for PAX9 and EpCAM was increased from 39.6 to 92.5% at day 9 ( Figure 3I ; Supplementary Figure 4D ), and these results were also verified by co-immunofluorescent staining for PAX9 and EpCAM ( Figure 3J ), highly consistent with those of SIX1 and EpCAM. Overall, we established a highly efficient approach to differentiate hESCs into 3PPE.

We developed six differentiation culture conditions to induce TECs ( Supplementary Table 1 ), and screened out the optimal culture one for the efficient differentiation towards TECs. Briefly, the differentiation medium for TECs was employed directly to 3PPE cells for a 4-day period at day 9 after the differentiation. At day 13 after the differentiation to TECs, the cells were collected for the analysis by RT-qPCR, and we found that the relative expression of FOXN1, a key transcription factor of TECs, was extremely higher under TEC-2 culture condition with over 10% of FOXN1⁺ cells in the differentiated cells when compared to those in other five culture conditions ( Figure 4A ; Supplementary Figure 5A ). MHC-II, a specific marker of mature TECs, was also up-regulated with the highest expression levels in differentiated cells under TEC-2 culture condition ( Figure 4B ; Supplementary Figure 5B ). Importantly, the proportion of MHC-II⁺EpCAM⁺ mature TECs could be reached to 16.8% ( Figure 4B ). TECs are usually divided into medullary TECs and cortical TECs. The former plays a role in the positive selection of T cells, while the latter makes negative selection of T cells and retains T cells that have no immunogenicity to host cells. We further found that the specific marker K5 of medullary TECs and the autoimmune regulator AIRE had the highest expression levels respectively in differentiated cells under TEC-2 culture condition ( Figure 4C ). Furthermore, we determined over 28% cells positive cells for K5 and 11% for AIRE by flow cytometry ( Figure 4D ; Supplementary Figure 5C ), and a large number of cells co-expressed MHC-II and AIRE using immunofluorescent staining ( Figure 4E ; Supplementary Figures 5F, G ), demonstrating that 3PPE cells had the ability to differentiate into medullary TECs. In addition, three major cortical TEC markers were also determined. Using flow cytometry, we found that the differentiated cells had 25% proportion of K8⁺ cells, and 18% of CD205⁺ and EpCAM⁺ cells, respectively ( Figure 4F ; Supplementary Figures 5A, E ), and the expression of K8 was further verified by immunofluorescent staining ( Figure 4G ). To further demonstrate the role in promoting T cell development, we co-cultured hESC-derived TECs with human hematopoietic stem/progenitor CD34⁺ cells (HSPC), and the cells were analyzed by flow cytometry one week after the co-culture. Surprisingly, over 20% of double positive CD45⁺CD3⁺ T cells were generated ( Figure 4H ; Supplementary Figure 5H ), almost 60% for CD4 single positive (SP) T cells, and over 30% for CD8 single positive T cells were included in CD45⁺CD3⁺ T cells ( Figure 4I , bottom panel). On the contrary, without the microenvironment support provided by TECs, only 8% of CD45⁺CD3⁺cells were differentiated from CD34⁺cells ( Figure 4H ), and all of them were induced to CD4⁺SP cells, and none of CD8⁺SP cells was generated ( Figure 4I , upper panel). Moreover, we further found that the proportions for TCRαβ⁺CD4⁺ double positive (DP) T cells was 56.8% ( Figure 4J , bottom panel), and for TCRαβ⁺CD8⁺DP T cells was 31.0% in CD45⁺CD3⁺ T cells in co-culture condition ( Figure 4K , bottom panel). However, without the microenvironment support provided by TECs, either CD4⁺SP cells or CD8⁺SP cells did not express TCRαβ ( Figures 4J, K , upper panel). Surprisingly, we found that mature CD3⁺TCRαβ⁺ DP T cells accounted for 27.6% in the total CD45+ cells ( Figure 4L ; Supplementary Figure 5H ), and CD3 positive cells were further verified using immunofluorescent staining ( Figure 4M ). Consistently, without the microenvironment support provided by TECs, only 2% CD3⁺TCRαβ⁺ were produced ( Figure 4L ). Collectively, our results not only demonstrated that 3PPE could be sufficient for further differentiation into functional TECs, also laid the fundamental basis for our further mechanism studies.

Given that HOXA3 expression was associated with the differentiation and formation of 3PPE, we wondered whether HOXA3 inhibition would affect the development of 3PPE. To this end, we used siRNA against HOXA3 during the differentiation of hESCs into 3PPE at day 7 ( Figure 5A ), and we found that HOXA3 was significantly reduced both at mRNA expression level ( Figure 5B ) and protein expression level ( Figures 5C, D ). This result was also verified by immunofluorescence analysis in which HOXA3 protein expression was relatively weak ( Supplementary Figures 6A, B ). Next, total RNAs of the cells was subjected to RNA sequencing, and the Heatmap showed that the expression levels of key transcription factors of 3PPE were all downregulated when HOXA3 was repressed ( Figure 5E ). In addition, we also further validated the results by qRT-PCR, and found that the repression of HOXA3 led to the downregulation of BMP4, SIX1, TBX1, PAX9, PBX1 ( Figure 5F ), and that protein levels also declined in synchrony ( Figures 5G, H ). Immunofluorescent staining was performed to further verify these findings, and we observed diminished immunofluorescent intensity of transcription factors SIX1 ( Figures 5I, J ) and PAX9 ( Supplementary Figures 6C, D ) upon HOXA3 inhibition. Thus, these lines of evidence indicated that HOXA3 acted upstream of these genes and a critical regulator of these key genes in the development of endoderm of the third pharyngeal pouch.

Next, we performed GO functional enrichment on the differential genes using transcriptome sequencing data which were obtained before and after HOXA3 inhibition in 3PPE, and found enrichments of biological processes such as cell population proliferation, animal organ development, and cell migration ( Figure 5K ), and also found that HOXA3 was involved in regulating migration ( Figure 5L ) and proliferation ( Figure 5M ) of epithelial cells by GSEA analysis. To this end, healing assays and transwell migration assays were performed to determine the effect of HOXA3 inhibition on the migration of 3PPE. Notably, the migratory ability of 3PPE cells was enhanced after HOXA3 inhibition with siRNA for 24 hours in healing assays ( Figure 5N ), and the number of cells that migrated through transwells 72 hours after the culture was much higher than those in the NC group in transwell migration assays ( Supplementary Figure 6E ). The Heat map results also showed significant upregulation of cell migration related genes after HOXA3 inhibition with siRNA ( Figure 5O ). CTHRC1 (26, 27) and INSL3 (28, 29) promote cell migration, the results of RT-qPCR assay exhibited the increased migratory capacity of cells by increasing the expressions of CTHRC1 and INSL3 following HOXA3 inhibition ( Figure 5P ). Moreover, we further performed cell cycle assay on 3PPE 48 hours after the transduction with siRNAs, and found that the inhibition of HOXA3 led to G1 arrest in the cell cycle ( Supplementary Figure 6F ), and the expression levels of proliferation related genes PCNA and Cdc25C were also significantly downregulated ( Supplementary Figure 6G ). Additionally, we also evaluated the role of HOXA3 in cell proliferation of 3PPE. Interestingly, using cell counting kit-8 assay, we found that the reduction of HOXA3 inhibited the proliferation of 3PPE cells ( Figure 5Q ), qRT-PCR results further verified that the expression levels of MCM2 and CDC45L as well as c-Myc, which are related to cell proliferation, were also significantly downregulated ( Figure 5R ). Taken together, our results revealed that the up-regulation of HOXA3 promoted the differentiation and proliferation and inhibited the migration, on the contrary, the down-regulation of HOXA3 enhanced the migration and suppressed the differentiation and proliferation, thus, these findings indicated that HOXA3 functioned as the on-off switch to regulate the differentiation and proliferation as well as the migration during the development of hESC-derived 3PPE.

The Wnt signaling pathway is inextricably linked to cell migration (30, 31), proliferation (32, 33), as well as differentiation (34, 35), thus, we next explored whether HOXA3 had a direct relationship with Wnt signaling pathway. As expected, we found that Wnt signaling pathway was significantly downregulated upon HOXA3 inhibition in GSEA analysis ( Figure 6A ), and these significant repressions included CTNNB1, LRP5, DVL3, APC, LRP6, TCF4, AXIN1, CCND1, WNT9A, WNT8B, WNT7A, WNT7B, WNT4, FRAT1, as well as GSK3β in heat map ( Figure 6B ). In order to clarify whether the Wnt pathway was influenced by HOXA3, we examined Wnt signaling related genes at both the protein and mRNA expression levels. qRT-PCR analysis showed that the knockdown of HOXA3 induced significant repressions of a large number of Wnt signaling related genes, including WNT9A, WNT1, FZD3, FZD8, FZD5, FRAT1, GSK3B, TCF12, CNND1 ( Figure 6C ). In addition, Western blotting was performed to further confirm that protein expression levels of FZD8, p-GSK3β, β-catenin, Cyclin D1 and c-Myc were obviously decreased in 3PPE after the treatment with HOXA3 siRNA ( Figures 6D, E ). EPHB2 is expressed in intestinal epithelial cells, and EPHB receptor is the target gene of Wnt signal (36), and EPHB receptors are thought to be positive regulators of proliferation in colon crypts (32). Thus, we further explored whether HOXA3 regulated EPHB2. First, we detected the relative expression of EPHB2 by qRT-PCR after inhibiting HOXA3, and found that EPHB2 was down-regulated ( Figure 6F ), and this down-regulation was verified at protein expression level by Western blot ( Figures 6G, H ). Importantly, we found that HOXA3 had a binding site upstream of the promoter of EPHB2 through JASPAR transcription factor motif database, and relative score is greater than 0.8 ( Figure 6I ; Supplementary Figure 7 ; Supplementary Table 2 ). Next, we determined the mediator which linked HOXA3 and Wnt signaling pathway by identifying the direct target of HOXA3, then chromatin immunoprecipitation (ChIP) assays were performed to use different primers that encompassed the EPHB2 promoter region for detecting the interaction between HOXA3 and EPHB2 promoter in 3PPE cells, and the results of ChIP PCR analysis ( Figure 6J ) and Chip-qPCR ( Figure 6K ) using primary antibody against anti-HOXA3 showed that DNA fragments treated with anti-HOXA3 antibody were significantly higher than those in IgG group, and that HOXA3 physically interacted with EPHB2, indicating EPHB2 as the target of HOXA3.

In order to clarify whether HOXA3 regulated Wnt signal through EPHB2, we designed rescue experiments to verify that EPHB2 served as vital downstream target gene of HOXA3. First, we used EPHB2 siRNA to determine the effect of EPHB2 on Wnt signaling pathway. To this end, we used siRNA against EPHB2 during the differentiation of hESCs into 3PPE at day 7, and we found that the expression of EPHB2 was significantly reduced ( Figure 7A ). Next, we further evaluated the consequence by qRT-PCR, and found that the repression of EPHB2 led to downregulation of GSK3B, FZD3, FZD8, TCF4, TCF12 ( Figure 7B ). Moreover, we further performed cell cycle assay on 3PPE cells 48 hours after siRNA transduction, and the results showed that the inhibition of EPHB2 led to G1 arrest in the cell cycle ( Supplementary Figures 8A, B ), and the expression levels of proliferation related genes PCNA and Cdc25C were also significantly downregulated ( Supplementary Figures 8C, D ). Additionally, we also evaluated the role of EPHB2 in the proliferation of 3PPE cells. Interestingly, using cell counting kit-8 assay, we found that the reduction of EPHB2 inhibited the proliferation of 3PPE cells ( Figure 7C ). Thus, healing assays were performed to determine the effect of EPHB2 inhibition on the migration of 3PPE cells. Notably, the migratory ability of 3PPE cells was enhanced after EPHB2 inhibition with siRNA for 24 hours in healing assays ( Figure 7D ), and the results of RT-qPCR assay exhibited the increased migratory capacity of cells by increasing the expressions of INSL3 following EPHB2 inhibition ( Figure 7E ). In order to determine targeting specificity of HOXA3 to the mediator of Wnt signaling pathway, the rescue experiments were performed. During the differentiation process of endodermal cells of 3PP, we first employed different concentrations of activator EFNB2 of EPHB2 at 100ng/mL, 200ng/mL and 500ng/mL to evaluate the enhancement of EPHB2 expression. The results showed that the addition of 200 ng/mL EFNB2 significantly increased the expression of EPHB2 and CNND1 in 3PPE cells ( Supplementary Figures 8E, F ). Next, we conducted qRT-PCR analysis for examining the rescue effect of activator EFNB2 of EPHB2 on the expressions of EPHB2 ( Figure 7F ) and markers of Wnt signaling pathway (FZD8, FZD5, TCF12, GSK3B, and CNND1) after the inhibition of HOXA3 with siRNA ( Figure 7G ). Consistently, HOXA3 siRNA profoundly down-regulated the expressions of these genes, however, the expression of these genes was increased or rescued by the addition of activator EFNB2 of EPHB2 ( Figures 7F, G ). Thus, these results demonstrated that HOXA3 regulated WNT pathway specifically through EPHB2. In order to further verify this mechanism, wounding healing assays were performed, and we observed a faster wound healing in the HOXA3 siRNA-treated group versus the NC group, indicating that the HOXA3 inhibition enhanced the migration of the cells ( Figure 7H ). Of note, the addition of EPHB2 activator EFNB2 protein was capable to markedly reduce cell migration which was promoted by HOXA3 inhibition ( Figure 7H ). Meanwhile, we also found that the migration gene INSL3 was significantly upregulated in response to HOXA3 inhibition, but was suppressed by the addition of EFNB2 ( Figure 7I ), indicating that EFNB2 prevented the cell migration under HOXA3 or EPHB2 inhibition. In addition to this phenomenon, cell cycle assays also were further performed using flow cytometry, and the results illustrated that EFNB2, an activator of EPHB2, was able to increase S phase in the cell cycle after the treatment with HOXA3 siRNA ( Figures 7J, K ), indicating that the extra expression of EFNB2 promoted the cell division (growth) after the prevention of the migration. Moreover, we conducted qRT-PCR analysis for examining the effects of EPHB2 on the expression of mitotic cell cycle markers PCNA and Cdc25c and cell proliferation markers CDC45L and MCM2. However, the addition of EFNB2 increased the expression cell cycle markers PCNA and Cdc25c ( Figure 7L ) and cell proliferation markers CDC45L and MCM2 ( Figure 7M ) in the cells of 3PP after the treatment with HOXA3 siRNA, further indicating that the extra expression of EFNB2 promoted the cell proliferation coupled with the inhibition of the migration under HOXA3 or EPHB2 inhibition. Taken together, these results demonstrated that the upregulation of HOXA3 by growth factors and small molecules induced the differentiation of 3PPE from hESCs and the expression of EPHB2 which activated Wnt signaling pathway to promote the cell proliferation of 3PPE, then the downregulation of HOXA3 switched its function to promote the cell migration, and suppressing the cell proliferation also through the inhibition of EPHB2-mediated Wnt pathway ( Figure 7N ).