Tissue-Targeted Transcriptomics Reveals SEMA3D Control of Hypoglossal Nerve Projection to Mouse Tongue Primordia

Abstract

The mouse hypoglossal nerve originates in the occipital motor nuclei at embryonic day (E)10.5 and projects a long distance, reaching the vicinity of the tongue primordia, the lateral lingual swellings, at E11.5. However, the details of how the hypoglossal nerve correctly projects to the primordia are poorly understood. To investigate the molecular basis of hypoglossal nerve elongation, we used a novel transcriptomic approach using the ROKU method. The ROKU algorithm identified 3825 genes specific for lateral lingual swellings at E11.5, of which 34 genes were predicted to be involved in axon guidance. Ingenuity Pathway Analysis-assisted enrichment revealed activation of the semaphorin signaling pathway during tongue development, and quantitative PCR showed that the expressions of Sema3d and Nrp1 in this pathway peaked at E11.5. Immunohistochemistry detected NRP1 in the hypoglossal nerve and SEMA3D as tiny granules in the extracellular space beneath the epithelium of the tongue primordia and in lateral and anterior regions of the mandibular arch. Fewer SEMA3D granules were localized around hypoglossal nerve axons and in the space where they elongated. In developing tongue primordia, tissue-specific regulation of SEMA3D might control the route of hypoglossal nerve projection via its repulsive effect on NRP1.

I.  Introduction

Tongue movement has crucial roles during swallowing, speech, and mastication and is controlled exclusively by the hypoglossal nerve [23, 29, 34, 40, 47]. The hypoglossal nerve innervates two main tongue muscle groups, defined anatomically as intrinsic and extrinsic. The intrinsic muscles comprise the bulk of the tongue and control its movement and morphology, whereas the extrinsic muscles are further subdivided into muscles that control protrusion or retraction of the tongue [7, 11].

In mice, the hypoglossal nerve originates in the hypoglossal nerve nuclei also known as the occipital motor nuclei. The nuclei differentiate from rhombomere 8 of the motor nuclei and are located in the dorsomedial medulla oblongata [13, 20, 30, 32]. By embryonic day (E) 10.5, the hypoglossal nerve ventrally leaves the medulla oblongata and begins to emerge as multiple branches heading towards the tongue primordia. The branches fasciculate to form a nerve bundle that reaches its target tongue location at E14.5 under proper axon guidance [5, 16, 17, 27, 37]. Axon guidance factors are molecules that provide spatiotemporal cues to growth cones of the developing nerves by their specific localization during embryonic development [33]. Secreted factors that function in the long range include netrin, slit, class 3 semaphorin, and various growth factors. Contact-mediated factors such as class 4–7 semaphorins and ephrin function in the short range [15]. Each of these factors has attractive or repulsive effects [9]. However, these molecules can trigger diverse and sometimes opposite axon responses, depending on the receptor composition and intracellular properties of the responding growth cone [48].

The details of how the hypoglossal nerve finds its route to the tongue primordia are poorly understood. Hepatocyte growth factor, produced from branchial arches, is considered a major chemoattractant for the guidance of the hypoglossal nerve [6, 36, 41]. However, because the hypoglossal nerve elongates in a restricted region in the developing tissue where many biological events occur concurrently, the involvement of other guidance factors has also been suggested [6, 36, 41]. While costly spatial transcriptomics have been developed, bulk-sample transcriptomics still allows transcriptional dynamics to be examined and is widely applicable in particular to survey tissue-specific events in small and scarce samples. In this study, we investigated the molecular events specific to hypoglossal nerve innervation by detecting tissue-specific genes in spatiotemporally distinct tissue samples.

II.  Materials and Methods

 Mice

Pregnant ICR mice were obtained from Charles River (Yokohama, Japan). We designated the morning of the day when a vaginal plug was observed as E0.5. All mice were kept under a 12 hr light/12 hr dark cycle with free access to standard laboratory feed and water. Animal care and use protocols were performed in accordance with approved institutional use and Nippon Dental University guidelines. To collect embryos, pregnant mice were humanely euthanized by cervical dislocation after inhalation anesthesia. Embryos were collected into fresh Hanks’ balanced salt solution (Invitrogen, Carlsbad, CA) and maxillofacial tissues were dissected from the head under a stereomicroscope. After collection, specimens were immediately fixed in 4% paraformaldehyde and then embedded in paraffin wax. These samples were used for histological examination. For gene expression analyses, each target sample, namely the LLS, the LMA, the MMA, and the T, was dissected en bloc and immediately immersed in RNAlater^TM Stabilization Solution (Invitrogen) to protect RNAs from RNase. Note that E11.5-LLS samples contained axon tips of the hypoglossal nerves and the hypoglossal nerve target region (Figure 2).

 Hematoxylin & eosin (H&E) staining and immunohistochemistry

Paraffin-embedded 4-μm thick sections were prepared and subjected to H&E staining or immunohistochemistry. H&E staining was performed using a standard protocol. For immunohistochemistry, sections were deparaffinized, made hydrophilic in ethanol, washed with distilled water and 0.01 M PBS, and then subjected to heat-induced epitope retrieval in Tris-EDTA (pH 9.0). Sections were then blocked with Protein Block Serum-Free Ready-to-use (X0909; Agilent, Santa Clara, CA) and incubated with primary antibodies at 4°C overnight. The primary antibodies used were: 1:5000 monoclonal mouse anti-GAP43 (G9264; Sigma-Aldrich, Burlington, MA), 1:200 polyclonal goat anti-NRP1 (AF566; R&D Systems, Minneapolis, MN), 1:100 polyclonal rabbit anti-SEMA3D (NBP1-85517; Novus Biologicals, Centennial, CO), and 1:500 Endomucin (V.7C7) mouse monoclonal antibody (sc-65495; Santa Cruz Biotechnology, Santa Cruz, CA). Immunoreactivity was detected using a VECTASTAIN^® Elite^® ABC-HRP Kit (Vector Laboratories, Newark, CA) and a DAB Substrate Kit (Vector Laboratories) or a fluorescent secondary antibody in combination with the TrueVIEW^® Autofluorescence Quenching Kit (Vector Laboratories) to suppress endogenous fluorescence. Nuclei were stained using SlowFade^TM Gold Antifade Mountant with DAPI (Thermo Fisher Scientific). The fluorescent secondary antibodies used were as follows: 1:100 Donkey anti-Rabbit IgG H&L Alexa Fluor^® 555 (ab150066; Abcam), 1:100 Donkey anti-Goat IgG (H + L) Cross-Adsorbed Secondary Antibody Alexa Fluor^TM 488 (A11055; Invitrogen), and 1:100 Donkey Anti-Mouse IgG H&L Alexa Fluor^® 647 (Abcam; ab150107). Fluorescence images were captured, then displayed and merged in pseudo-color using ImageJ software (NIH, Bethesda, MD) manipulated with ImageJ macros. We have tested the positive reactivity of each primary antibody by using tissue sections that contain molecules known to be reacted. For each detection, we also set negative controls by applying PBS to tissue sections instead of the primary antibody.

 Histology-based 3D reconstruction

The digitization of serial width-ranged images was performed as described previously. Briefly, the high-resolution digitization of consecutive immunolabeled serial sections was achieved using a NanoZoomer HT virtual microscope (Hamamatsu Photonics, Hamamatsu, Japan). Image registration and segmentation were performed using ImageJ software, and a stack of serial images was obtained for 3D reconstruction. The architecture of hypoglossal nerve axons and other cranial nerves was highlighted by painting them with different pseudo-colors in the 3D images. The images obtained were superimposed on tissue outlines.

 Tissue collection for tissue-specific gene expression analysis

Tissues were dissected and pooled together in a single tube for microarray analysis. The following were collected into single tubes: E9.5 median mandibular arches (E9.5-MMA; 12 embryos/tube), E11.5 lateral lingual swellings (E11.5-LLS; 4 embryos/tube), E11.5 lateral mandibular arches (E11.5-LMA; 4 embryos/tube), and E14.5 median tongue (E14.5-T; 4 embryos/tube).

 Microarray analysis

Total RNA was extracted from each sample using TRIzol^® Reagent (Thermo Fisher Scientific) and a miRNeasy Mini Kit (Qiagen, Hilden, Germany) according to the manufacturers’ instructions. RNA samples were quantified using an ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, Germany) and their quality was confirmed using an Agilent 4200 TapeStation (Agilent Technologies, USA). For microarray analysis, cRNA was amplified, labeled, and then hybridized to a 60K Agilent 60-mer oligomicroarray according to the manufacturer’s instructions. Microarray slides were scanned using an Agilent scanner and relative hybridization intensities and background hybridization values were calculated using Agilent Feature Extraction software (Ver.9.5.1.1). Our microarray data can be obtained in the Gene Expression Omnibus database (accession #GSE249490).

 Tissue-specific gene analysis

Quantile normalization was performed on the raw data using Python (version 3.7.9) and lincRNAs present in the data were removed. To detect tissue-specific genes in spatiotemporal comparisons, the ROKU method [19] was used in R (version 4.0.3). Duplicate gene symbols were combined into single symbols. QuickGO (EMBL-EBI) [4] was used to obtain the gene list and was assigned the annotation “GO: 0007411: axon guidance”. An intersection of the tissue-specific genes and their gene lists was then performed. Ingenuity Pathway Analysis (Qiagen) was performed for downstream analysis and Python scripts were used to visualize heatmaps by calculating Z-score normalized signal values of candidate genes.

 Quantitative RT-PCR (qPCR)

The median region of the tongue primordia was collected between E10.5 and E18.5. Total RNA was isolated and purified from each sample using an miRNeasy Micro Kit (Qiagen) according to the manufacturer’s protocol. The quantity and purity of isolated RNA were assessed using a spectrophotometer and cDNA was synthesized from <1 μg RNA using ReverTra Ace^® qPCR RT Master Mix with gDNA Remover (Toyobo, Osaka, Japan). Custom primers were designed across exon-exon junctions using the National Center for Biotechnology Information (NCBI) Primer-BLAST program. The primer sequences are described in Table 1. The relative expression levels of target genes were quantified using the StepOne Real-time PCR system (Thermo Fisher Scientific) and THUNDERBIRD^® Probe qPCR Mix (Toyobo). Target genes were normalized against an internal control gene (Gapdh) and the expression level was calculated using the delta-delta C_T method. The experiment was repeated at least three times.

Table 1. 

PCR primers

Gene/Direction	Sequence	Length (bp)
mSema3d Fw	TTGCTAGCAGGAAGGGTGAAC	21
mSema3d Rv	ACTTCTGGATCTTGGGTTCTCATC	24
mNrp1 Fw	AGGAGGCACCACTGTCCTG	19
mNrp1 Rv	AACCGTATGTCGGGAACTCTGA	22
mPlxnd1 Fw	CAGTACTGCCCTCGGAGATTG	21
mPlxnd1 Rv	TCGAACGCCATTTCCATAGTC	21
mPlxna2 Fw	AGCCCAGAGATCCCAGTGAA	20
mPlxna2 Rv	CTTGGCGCCACTCCAGATCC	20

 Statistical analysis

Statistical analysis was performed for time-series qPCR data. The Tukey-Kramer test was performed using Python library statsmodels (version 0.14.0) [39] and P-values < 0.01 were considered statically significant and denoted with an asterisk in bar plots.

III.  Results

 The target region of the hypoglossal nerve is defined by the spatial configuration of cranial nerves at E11.5

The spatial localization of cranial nervous system axons was reconstructed from serial E11.5 embryo sections immunostained for the neuronal marker GAP43 (Figure 1A). GAP43 is a vertebrate neuron-specific protein involved in axon growth [22]. As shown in Figure 1B, the hypoglossal nerve developed from several nerve nuclei, bent in a ventral direction, and then projected in a cranial direction to reach the central area close to the tongue primordia. Only the GAP43-positive hypoglossal nerve (shown in red) extended its axons to that area. Other GAP43-positive nerves (trigeminal and facial nerves; shown in yellow) remained in lateral areas and did not approach the tongue primordia. On the basis of these observations, we defined the central area of the first branchial arch at E11.5, which does not include the mandibular arch itself, as the “hypoglossal nerve target region”.

Fig. 1.

Three-dimensional view of mouse cranial nerves at E11.5. (A) Data collection from histological sections. Digitized serial images of GAP43-immunostained sections were processed to produce a steric structure. (B) Histology-based 3D reconstruction of the hypoglossal nerve (HgN; shown in red) and other cranial nerves (shown in yellow). The embryo surface is shown in purple. The hypoglossal nerve target region defined in this study is indicated by an asterisk. pLLS, posterior LLS; HgN, hypoglossal nerve; CT, chorda tympani; MN, mandibular nerve; MA, mandibular arch.

 E11.5-LLS expresses 34 tissue-specific genes related to axon guidance

We hypothesized that the hypoglossal nerve target region presets guidance factors in the tissue space to exclusively attract the hypoglossal nerve. To identify such factors, we performed DNA microarray analysis and used the ROKU method, an algorithm that detects tissue-specific genes by comparing expression patterns across different periods and tissue regions. As illustrated in Figure 2A, we regarded each of these regions as an independent tissue unit for ROKU analysis. We set the lateral lingual swellings (LLS) including the hypoglossal nerve target region at E11.5, “E11.5-LLS”, as the main target. We also set the lateral mandibular arches (LMA) at E11.5, “E11.5-LMA”, on the spatial axis for comparison with tissues from different sites over the same period to cancel out the influence of guidance factors for other cranial nerves. In addition, to compare tissues from the same site at different time points, the medial mandibular arch (MMA) at E9.5 “E9.5-MMA” and the tongue (T) at E14.5 “E14.5-T” were set on the temporal axis of E11.5-LLS.

From the signal data of 52,146 microarray probes, ROKU detected 3825 tissue-specific genes in E11.5-LLS. In conjunction with 191 entries in the “axon guidance” Gene Ontology (GO) terms, 34 genes were designated as E11.5-LLS-specific and relevant to axon guidance (Figure 2B and Table 2).

Fig. 2.

Experimental settings for tissue-specific gene expression analysis. (A) Preparation of tissue samples. In the schematic representation of bulk samples, the mandibular arch tissue is shown in yellow, and the tongue primordia and developing tongue are colored pink. Each tissue was collected from several mouse embryos and pooled before processing. E11.5-LLS includes developing hypoglossal nerves as illustrated and is the main target of this study. As represented in the scheme, the temporal axis (E9.5-MMA and E14.5-T) and spatial axis (E11.5-LMA) were set for microarray analysis. (B) The number of tissue-specific genes detected by the ROKU method is listed. A Venn diagram shows genes that are common for the E11.5-LLS-specific and Gene Ontology term “axon guidance” parameters.

Table 2. 

List of E11.5-LLS-specific genes in relation to axon guidance

ProbeID	Gene Symbol	E9.5-MMA	E11.5-LLS	E11.5-LMA	E14.5-T	modH	ranking
A_52_P30152	Robo1	0	1	0	0	0.094	434
A_55_P2731571	Nrcam	0	1	0	0	0.216	1316
A_66_P122015	Atoh1	0	1	0	0	0.435	2723
A_52_P295822	Robo2	0	1	0	0	0.521	3284
A_66_P100845	Etv1	0	1	0	0	0.831	6056
A_55_P2736540	Sema3d	0	1	0	0	0.841	6157
A_52_P371237	Nrp1	0	1	0	0	0.871	6521
A_52_P128134	Foxd1	0	1	0	0	1.105	10953
A_55_P2032649	Kif5b	0	1	0	0	1.19	14036
A_52_P515857	Reln	0	1	0	0	1.201	14614
A_55_P2722020	Tenm2	0	1	0	0	1.458	25185
A_55_P2745553	Epha3	0	1	0	0	1.535	27768
A_55_P2715809	Lgi1	0	1	0	0	1.551	28341
A_55_P1969058	Epha7	0	1	0	0	1.903	45013
A_55_P2144556	Flrt3	0	1	0	0	1.919	46029
A_51_P189943	Dscam	0	−1	0	0	0.644	4184
A_51_P365369	Reln	0	−1	0	0	0.892	6807
A_51_P305320	B4gat1	0	−1	0	0	1.119	11368
A_55_P2931015	Sema6d	0	−1	0	0	1.323	20536
A_52_P661	Plxna4	0	−1	0	0	1.369	22357
A_55_P2100555	L1cam	0	−1	0	0	1.391	23067
A_65_P06809	Klf7	0	−1	0	0	1.402	23414
A_55_P2805510	Gas1	0	−1	0	0	1.448	24890
A_51_P237040	Nog	0	−1	0	0	1.554	28470
A_51_P494992	Runx3	0	−1	0	0	1.604	30324
A_51_P489935	Pax6	0	−1	0	0	1.607	30423
A_52_P282000	Enah	0	−1	0	0	1.67	32882
A_52_P258116	Wnt3	0	−1	0	0	1.782	37893
A_55_P2076489	Rac1	0	−1	0	0	1.794	38521
A_55_P2908134	Kif5b	0	−1	0	0	1.797	38660
A_55_P2168023	Sema3g	0	−1	0	0	1.805	39045
A_52_P597791	Robo2	0	−1	0	0	1.833	40522
A_52_P219314	Vasp	0	−1	0	0	1.833	40569
A_51_P342877	Scn1b	0	−1	0	0	1.844	41245

Tissue-specific genes identified by the ROKU algorithm were combined with the Gene Ontology term “axon guidance” and ranked for specificity in samples by Shannon entropy. The resulting 34 genes are shown with their rankings. In the E11.5-LLS column, 1 and −1 indicate upregulation and downregulation in the corresponding tissue, respectively.

A hierarchical clustering heatmap demonstrated these 34 genes to be E11.5-LLS-specific (Figure 3A). Ingenuity Pathway Analysis (IPA)-assisted enrichment revealed that semaphorin signaling and ephrin signaling were especially activated (Figure 3B). The network of these enriched terms indicated that upregulated pathways formed a semaphorin-related cluster that is connected to the “axon guidance” hub node, whereas downregulated pathways exhibited a confounding network (Figure 3C). Notably, genes for the ligand, semaphorin 3D (Sema3d), and its receptor, neuropilin 1 (Nrp1), were upregulated in the same period (Figure 3D) and ranked within the top 20% of E11.5-LLS-specific genes (Table 2). These data indicated that Sema3d and Nrp1 genes are specifically expressed at E11.5-LLS and are therefore candidates for axon guidance of the hypoglossal nerve.

Fig. 3.

Detailed profiles of E11.5-LLS-specific and axon guidance-related genes. (A) Heatmap showing 34 genes (normalized by Z-score) that were specifically upregulated or downregulated in E11.5-LLS. (B) Enrichment analysis of the 34 genes by Ingenuity Pathway Analysis (IPA). The “Axonal guidance signaling” pathway, as well as sub-pathways involved in neural guidance, are enriched. N, number of E11.5-LLS-specific genes attributed to the corresponding pathway. (C) Networks representing the association between the enriched terms (nodes) and the overlapping genes (edges). The number in each node corresponds to the number in B. (D) Schematic representation of the semaphorin signaling pathway enriched in IPA. The ligand (Sema3d) and receptor (Nrp1) are upregulated.

 Significant upregulation of Sema3d and Nrp1 genes in the median region of tongue primordia at E11.5

The expression profiles of Sema3d and Nrp1 genes in the developing tongue (median region of the tongue primordia at E10.5–E18.5) were validated by real-time qPCR. Sema3d and Nrp1 exhibited significant peak expressions at E11.5. We also quantified the expression levels of Plxnd1 and Plxna2, which encode the SEMA3D co-receptors PLXND1 and PLXNA2, respectively, and found their peak expressions were also at E11.5 (Figure 4).

Fig. 4.

Chronological expression profiles of semaphorin signaling molecules in the developing mouse tongue. E11.5-LLS-specific Sema3d and Nrp1 expression was validated by real-time qPCR. Plxnd1 and Plxna2, which belong to the semaphorin signaling pathway, were also analyzed. Their encoded proteins form a heteromeric complex with NRP1. The expression level relative to that at E10.5 is indicated for each time point. *p < 0.01 (Tukey-Kramer); bars, mean ± 1 SD.

 Localization of SEMA3D around the developing route of the NRP1-positive hypoglossal nerve

The distribution of NRP1 and SEMA3D proteins was verified by immunohistochemistry using frontal and sagittal sections of E11.5 embryos (Figures 5A and 6A). NRP1-positive signals were observed in endomucin-expressing endothelial cells (Figure 5B), stromal cells around the ridges (Figure 5C), and GAP43-positive cranial nerves (Figures 5D, 5F, 6B, and 6C). SEMA3D-positive signals were scattered in the extracellular space as dot-like granules, some of which were distributed along the NRP-positive cell surfaces (Figure 5C).

Fig. 5.

Localization of GAP43 (neuronal marker), NRP1, and SEMA3D in tongue primordia at E11.5 (frontal view). (A) Schematic illustration of the E11.5 mouse embryo. The mandibular arch is shown in yellow and lateral lingual swellings, i.e., tongue primordia, are colored pink. Antero-posterior positions of cross-sections are shown. Each section corresponds to panels D to F below. (B) Double immunostaining shows the colocalization of NRP1 with endomucin (endothelial cell marker). Bar = 20 μm. (C) Double immunostaining image of NRP1 and SEMA3D. SEMA3D is observed as dot-like granules in the intercellular space or NRP1-positive cell surface. Bar = 20 μm. (D–F) Frontal sections of (D) the posterior part of the LLS, (E) middle part of the LLS, and (F) anterior part of the mandibular arch. NRP1 is detected in GAP43-positive nerves and stroma. SEMA3D is observed as dot-like granules and large cell-like spots (arrowheads). Bars in D–F = 150 μm; XII, hypoglossal nerve; V, trigeminal nerve.

In the frontal view, the GAP43/NRP1-double-positive hypoglossal nerve was present in the posterior part of the LLS. A distinct NRP1-positive cell population was observed beneath the lateral epithelium of the LLS (Figure 5D and 5E) but was absent in the anterior part of the mandibular arch (Figure 5F). A substantial number of SEMA3D-positive granules were distributed in the LLS, just underneath the epithelium, but fewer such granules were observed in the vicinity of and below the hypoglossal nerve (Figure 5D and 5E). SEMA3D was broadly localized to the anterior part of the mandibular arch and its upper epithelium (Figure 5F).

In the sagittal view, the NRP1-positive trigeminal nerves were lateral to the tongue primordia (Figure 6B). The hypoglossal nerve approaching the LLS from below was observed on the inner side of the LLS (Figure 6C), but NRP1-positive signals were scarce on the primordia midline (Figure 6D). SEMA3D signals were spread around the slightly dorsolateral area of the mandibular arch (Figure 6B). Consistent with the frontal view, there were fewer SEMA3D signals around the hypoglossal nerve in the LLS but they were abundant in the sub-epithelial region of the LLS and mandibular arch (Figure 6C). In the near-midline region, i.e., the hypoglossal nerve target region, very few SEMA3D signals were observed, except for in the sub-epithelium (Figure 6D, asterisk). Of note, there were some spots of concentrated SEMA3D signal (Figure 5D, 5F, 6B–6D, arrowheads).

Fig. 6.

Localization of GAP43, NRP1, and SEMA3D in tongue primordia at E11.5 (sagittal view). (A) Schematic illustration of the E11.5 mouse embryo. The mandibular arch is shown in yellow and lateral lingual swellings, i.e., tongue primordia, are colored pink. Mesio-distal position of sagittal sections. Each section corresponds to panels B to D below. Lateral part of the LLS (B), the hypoglossal nerve target region at the site of hypoglossal nerve entry (C), and a site close to the midline (D). NRP1 is detected in GAP43-positive nerves, including the hypoglossal nerve (XII), and stroma. SEMA3D is observed as tiny granules in the extracellular space or in large cell-like spots (arrowheads). V, trigeminal nerve; VII, facial nerve. *Central portion of the hypoglossal nerve target region. Bar = 150 μm.

IV.  Discussion

We performed a histology-based three-dimensional (3D) reconstruction of the mouse embryonic head to acquire the trajectory of hypoglossal nerve axons. The anatomical location of the hypoglossal nerves at E11.5 was consistent with previous immunohistological observation using the GAP43 marker [42]. The axon tips (growth cones) of the hypoglossal nerves reached the lower central region of the tongue primordia. None of the other cranial nerves penetrated this region; therefore, we focused on it as the hypoglossal nerve target region and investigated the molecular mechanism that underlies the projection of the hypoglossal nerve to this tissue space.

In this study, we used a novel approach to analyze gene expression in a target tissue by preparing distinct tissues that sandwich the target tissue spatially and temporally. Generally, to screen for candidate genes in complex molecular embryology [31], differentially expressed genes (DEGs) are detected by comparing two groups via transcriptome analysis [3]. However, this method cannot be applied to detect tissue-specific genes directly [49]. Multiple pairwise comparisons should be avoided because the results are likely to be significantly different [12]. Furthermore, the makeshift adjustment of thresholds for DEG detection should be avoided. Therefore, we performed tissue-specific gene analysis using the ROKU method. The ROKU method calculates entropy for individual genes and detects outliers for each tissue using Akaike’s information criterion without adjusting thresholds. Attention should be paid to dispersion among samples when replicates cannot be prepared because of limited sample availability [19, 49].

We obtained genes that were E11.5-LLS-specific and axon guidance-related. Sema3d and Nrp1 in the semaphorin pathway were significantly upregulated. Semaphorin signaling pathways have diverse roles in inflammation, immune response, bone metabolism, tumor progression and invasion, angiogenesis and neurogenesis [21, 24, 28]. SEMA3D is a secreted protein that belongs to class 3 semaphorins [1, 18]. After binding to its receptor NRP1, the SEMA3D signal alters the cytoskeleton and cell adhesion to exert a repulsive effect and redirect cell migration [2, 26, 43, 44]. It is particularly interesting to find a repulsive signal that predominates in the hypoglossal nerve target region rather than attractive factors, which are expected in axon guidance. We speculate that attractive factors are commonly active in oral tissues, where many nerves develop, and therefore these did not fall into the tissue-specific framework of our ROKU analysis.

Immunohistochemical detection of the protein products revealed that GAP43-positive cranial nerves, including the trigeminal, facial, and hypoglossal nerves, were all positive for NRP1. NRP1 was also detected in endothelial cells and stromal cells. However, the SEMA3D protein was distributed in the embryonic tissues as dot-like granules, which resided in intercellular space. Another class 3 semaphorin cognate, SEMA3A, has been detected as similar dot-like granules trapped in the extracellular matrix or neurons in the brain [10, 50]. Therefore, secreted SEMA3D was probably also trapped in the extracellular matrix.

At E11.5, SEMA3D-positive granules were broadly distributed beneath the epithelium of the tongue primordia, as well as lateral and anterior to the mandibular arch. However, few SEMA3D-positive granules were localized around the hypoglossal nerve axons. Importantly, in the hypoglossal nerve target region, SEMA3D-positive granules were scattered in sub-epithelial connective tissue, but many fewer granules were observed in areas deeper than the target region. SEMA3D acts as a repulsive factor; therefore, it is reasonable that the hypoglossal nerve projects along the tissue space where SEMA3D granules are poorly localized.

Although SEMA3D-positive granules had a definite spatial configuration in the primordia, some questions remain, such as which cells produce SEMA3D and how the axon innervation route is constructed. A plausible candidate for the SEMA3D-producing cell lineage is cranial neural crest cells, which are involved in the formation of tongue connective tissue and vasculature [35]. A previous in situ hybridization study reported that neural crest cells in the chick branchial arch expressed Sema3d [8]. Our immunohistochemistry identified some large spots of SEMA3D immunoreactivity, indicating another possible source of SEMA3D if these spots represent living cells.

Our qPCR results revealed a substantial level of Sema3d expression from E10.5 to E11.5, indicating that the spatial arrangement of SEMA3D might be established within this short period. The axons from ventral motor nuclei, including those of the hypoglossal nerve, extend ventrolaterally with reduced sensitivity to repulsive factors due to the transient expression of Cxcr4 until E10.5 [25]. Therefore, E11.5 is the likely timepoint where the NRP1-expressing hypoglossal nerve is driven in the direction of the hypoglossal nerve target region by the repulsive effect of SEMA3D. Concurrently, this environment might sequester the hypoglossal nerve from other NRP1-expressing nerves because E11.5 is also the timepoint when other cranial nerves enter the tongue primordia from the lateral side [42]. The selective removal of NRP1 from motor neurons leads to defects in the fasciculation of the hypoglossal nerve and a decrease in neural crest-derived Schwann cells [16]. These findings indicate that the mutual interaction between an SEMA3D-equipped environment and NRP1-expressing hypoglossal nerve axons is indispensable for tongue functions.

The NRP1 receptor does not function alone but with a plexin family molecule, PLXND1, which interacts with SEMA3D to transduce its signal downstream. PLXNA2 also acts as a SEMA3D receptor where it functions as a co-receptor with NRP1 [14, 38, 45, 46]. Although Plxnd1 and Plxna2 were not prominent in E11.5-LLS in our ROKU-assisted analysis, qPCR revealed the significant upregulation of both genes in the median tongue primordia at E11.5.

Not all of the above findings were expected at the beginning of this study. The use of the ROKU method in our assay resulted in too many tissue-specific genes being obtained with the outlier detection of Akaike’s information criterion. However, in combination with a ranking calculation from the Shannon entropy and Gene Ontology analysis, tissue specificity can be determined as demonstrated. As long as a balance is maintained between sensitivity and specificity, this approach is a powerful tool that can identify potentially important genes and shed light on fundamental properties of the tissue environment.

V.  Conflicts of Interest

The authors declare that there are no conflicts of interest.

VI.  Acknowledgments

This work was supported by JSPS KAKENHI (grant numbers 18K09530 and 21K09822). We thank Jeremy Allen, PhD, and J. Ludovic Croxford, PhD, from Edanz (https://jp.edanz.com/ac) for editing a draft of this manuscript.

VII. References

• 1   Aghajanian,  H.,  Choi,  C.,  Ho,  V. C.,  Gupta,  M.,  Singh,  M. K. and  Epstein,  J. A. (2014) Semaphorin 3d and semaphorin 3e direct endothelial motility through distinct molecular signaling pathways. J. Biol. Chem. 289; 17971–17979.
• 2   Alto,  L. T. and  Terman,  J. R. (2017) Semaphorins and their signaling mechanisms. Methods Mol. Biol. 1493; 1–25.
• 3   Anjum,  A.,  Jaggi,  S.,  Varghese,  E.,  Lall,  S.,  Bhowmik,  A. and  Rai,  A. (2016) Identification of Differentially Expressed Genes in RNA-seq Data of Arabidopsis thaliana: A Compound Distribution Approach. J. Comput. Biol. 23; 239–247.
• 4   Binns,  D.,  Dimmer,  E.,  Huntley,  R.,  Barrell,  D.,  O’Donovan,  C. and  Apweiler,  R. (2009) QuickGO: a web-based tool for Gene Ontology searching. Bioinformatics 25; 3045–3046.
• 5  Catala, M. and Kubis, N. (2013) Gross anatomy and development of the peripheral nervous system. in “Handbook of Clinical Neurology,” 1st ed., pp. 29–41, Elsevier B.V., 115; 29–41.
• 6   Caton,  A.,  Hacker,  A.,  Naeem,  A.,  Livet,  J.,  Maina,  F.,  Bladt,  F., et al. (2000) The branchial arches and HGF are growth-promoting and chemoattractant for cranial motor axons. Development 127; 1751–1766.
• 7   Chen,  X.,  Wang,  J. W.,  Salin-Cantegrel,  A.,  Dali,  R. and  Stifani,  S. (2016) Transcriptional regulation of mouse hypoglossal motor neuron somatotopic map formation. Brain Struct. Funct. 221; 4187–4202.
• 8   Chilton,  J. K. and  Guthrie,  S. (2003) Cranial Expression of Class 3 Secreted Semaphorins and Their Neuropilin Receptors. Dev. Dyn. 228; 726–733.
• 9   Comer,  J. D.,  Alvarez,  S.,  Butler,  S. J. and  Kaltschmidt,  J. A. (2019) Commissural axon guidance in the developing spinal cord: From Cajal to the present day. Neural Dev. 14; 1–16.
• 10   de Wit,  J.,  Toonen,  R. F.,  Verhaagen,  J. and  Verhage,  M. (2006) Vesicular trafficking of semaphorin 3A is activity-dependent and differs between axons and dendrites. Traffic 7; 1060–1077.
• 11   Fregosi,  R. F. (2011) Respiratory related control of hypoglossal motoneurons—Knowing what we do not know. Respir. Physiol. Neurobiol. 179; 43–47.
• 12   Grandhi,  A.,  Guo,  W. and  Peddada,  S. D. (2016) A multiple testing procedure for multi-dimensional pairwise comparisons with application to gene expression studies. BMC Bioinformatics. 17; 1–12.
• 13   Guthrie,  S. (2007) Patterning and axon guidance of cranial motor neurons. Nature reviews. Neuroscience 8; 859–871.
• 14   Hamm,  M. J.,  Kirchmaier,  B. C. and  Herzog,  W. (2016) Sema3d controls collective endothelial cell migration by distinct mechanisms via nrp1 and plxnD1. J. Cell Biol. 215; 415–430.
• 15   Hinck,  L. (2004) The versatile roles of “axon guidance” cues in tissue morphogenesis. Dev. Cell. 7; 783–793.
• 16   Huettl,  R. E. and  Huber,  A. B. (2011) Cranial nerve fasciculation and Schwann cell migration are impaired after loss of Npn-1. Dev. Biol. 359; 230–241.
• 17   Huettl,  R. E.,  Eckstein,  S.,  Stahl,  T.,  Petricca,  S.,  Ninkovic,  J.,  Götz,  M., et al. (2016) Functional dissection of the Pax6 paired domain: Roles in neural tube patterning and peripheral nervous system development. Dev. Biol. 413; 86–103.
• 18   Jurcak,  N. R.,  Rucki,  A. A.,  Muth,  S.,  Thompson,  E.,  Sharma,  R.,  Ding,  D., et al. (2019) Axon Guidance Molecules Promote Perineural Invasion and Metastasis of Orthotopic Pancreatic Tumors in Mice. Gastroenterology 157; 838–850.
• 19   Kadota,  K.,  Ye,  J.,  Nakai,  Y.,  Terada,  T. and  Shimizu,  K. (2006) ROKU: A novel method for indentification of tissue-specific genes. BMC Bioinformatics 7; 1–9.
• 20   Kanjhan,  R.,  Fogarty,  M. J.,  Noakes,  P. G. and  Bellingham,  M. C. (2016) Developmental changes in the morphology of mouse hypoglossal motor neurons. Brain Struct. Funct. 221; 3755–3786.
• 21   Kanth,  S. M.,  Gairhe,  S. and  Torabi-Parizi,  P. (2021) The Role of Semaphorins and Their Receptors in Innate Immune Responses and Clinical Diseases of Acute Inflammation. Front. Immunol. 12; 1–9.
• 22   Kawasaki,  A.,  Okada,  M.,  Tamada,  A.,  Okuda,  S.,  Nozumi,  M.,  Ito,  Y., et al. (2018) Growth Cone Phosphoproteomics Reveals that GAP-43 Phosphorylated by JNK Is a Marker of Axon Growth and Regeneration. iScience 4; 190–203.
• 23   Koshino,  H.,  Hirai,  T.,  Ishijima,  T. and  Ikeda,  Y. (1997) Tongue motor skills and masticatory performance in adult dentates, elderly dentates, and complete denture wearers. J. Prosthet. Dent. 77; 147–152.
• 24   Li,  Z.,  Hao,  J.,  Duan,  X.,  Wu,  N.,  Zhou,  Z.,  Yang,  F., et al. (2017) The role of semaphorin 3A in bone remodeling. Front. Cell. Neurosci. 11; 1–8.
• 25   Lieberam,  I.,  Agalliu,  D.,  Nagasawa,  T.,  Ericson,  J. and  Jessell,  T. M. (2005) A Cxcl12-Cxcr4 chemokine signaling pathway defines the initial trajectory of mammalian motor axons. Neuron 47; 667–679.
• 26   Liu,  Y.,  Berndt,  J.,  Su,  F.,  Tawarayama,  H.,  Shoji,  W.,  Kuwada,  J. Y., et al. (2004) Semaphorin3D Guides Retinal Axons along the Dorsoventral Axis of the Tectum. J. Neurosci. 24; 310–318.
• 27   Loh,  C.,  Maya,  M. M. and  Go,  J. L. (2002) Cranial nerve XII: The hypoglossal nerve. Semin. Ultrasound CT MRI. 23; 256–265.
• 28   Mastrantonio,  R.,  You,  H. and  Tamagnone,  L. (2021) Semaphorins as emerging clinical biomarkers and therapeutic targets in cancer. Theranostics 11; 3262–3277.
• 29   Melong,  J.,  Bezuhly,  M. and  Hong,  P. (2021) The Effect of Tongue-Tie Release on Speech Articulation and Intelligibility. Ear, Nose Throat J. 0; 1–5.
• 30   Müller,  F. and  O’Rahilly,  R. (2011) The initial appearance of the cranial nerves and related neuronal migration in staged human embryos. Cells Tissues Organs. 193; 215–238.
• 31   Oligny,  L. L. (2001) Human molecular embryogenesis: An overview. Pediatr. Dev. Pathol. 4; 324–343.
• 32   O’Rahilly,  R. and  Müller,  F. (1984) The early development of the hypoglossal nerve and occipital somites in staged human embryos. Am. J. Anat. 169; 237–257.
• 33   Onesto,  M. M.,  Short,  C. A.,  Rempel,  S. K.,  Catlett,  T. S. and  Gomez,  T. M. (2021) Growth Factors as Axon Guidance Molecules: Lessons From in vitro Studies. Front. Neurosci. 15; 1–17.
• 34   Ono,  T.,  Hori,  K.,  Tamine,  K. and  Maeda,  Y. (2009) Evaluation of tongue motor biomechanics during swallowing-From oral feeding models to quantitative sensing methods. Jpn. Dent. Sci. Rev. 45; 65–74.
• 35   Parada,  C.,  Han,  D. and  Chai,  Y. (2012) Molecular and cellular regulatory mechanisms of tongue myogenesis. J. Dent. Res. 91; 528–535.
• 36   Raper,  J. and  Mason,  C. (2010) Cellular strategies of axonal pathfinding. Cold Spring Harb. Perspect. Biol. 2; 1–21.
• 37   Reissig,  L. F.,  Geyer,  S. H.,  Winkler,  V.,  Preineder,  E.,  Prin,  F.,  Wilson,  R., et al. (2022) Detailed characterizations of cranial nerve anatomy in E14.5 mouse embryos/fetuses and their use as reference for diagnosing subtle, but potentially lethal malformations in mutants. Front. Cell Dev. Biol. 10; 1–17.
• 38   Schwarz,  Q. and  Ruhrberg,  C. (2010) Neuropilin, you gotta let me know: Should I stay or should I go? Cell Adhes. Migr. 4; 61–66.
• 39  Seabold, S. and Perktold, J. (2010) Statsmodels: Econometric and Statistical Modeling with Python. Proc. 9th Python Sci. Conf.; pp. 92–96.
• 40   Abd-El-Malek,  S. (1955) The part played by the tongue in mastication and deglutition. J. Anat. 89; 250–254.
• 41   Short,  C. A.,  Onesto,  M. M.,  Rempel,  S. K.,  Catlett,  T. S. and  Gomez,  T. M. (2021) Familiar growth factors have diverse roles in neural network assembly. Curr. Opin. Neurobiol. 66; 233–239.
• 42   Sugimoto,  T.,  Taya,  Y.,  Shimazu,  Y.,  Soeno,  Y.,  Sato,  K. and  Aoba,  T. (2015) Three-Dimensional Visualization of Developing Neurovascular Architecture in the Craniofacial Region of Embryonic Mice. Anat. Rec. 298; 1824–1835.
• 43   Takahashi,  K.,  Ishida,  M. and  Takahashi,  H. (2009) Expression of Sema3D in subsets of neurons in the developing dorsal root ganglia of the rat. Neurosci. Lett. 455; 17–21.
• 44   Takahashi,  K.,  Tomizawa,  K.,  Ishida,  M.,  Hirokawa,  K. and  Takahashi,  H. (2009) Identification and tissue-specific expression of a novel isoform of Semaphorin 3D. Biochim. Biophys. Acta - Gen. Subj. 1790; 395–400.
• 45   Toledano,  S.,  Nir-Zvi,  I.,  Engelman,  R.,  Kessler,  O. and  Neufeld,  G. (2019) Class-3 Semaphorins and Their Receptors: Potent Multifunctional Modulators of Tumor Progression. Int. J. Mol. Sci. 20; 556.
• 46   Visser,  J. J.,  Cheng,  Y.,  Perry,  S. C.,  Chastain,  A. B.,  Parsa,  B.,  Masri,  S. S., et al. (2015) An extracellular biochemical screen reveals that FLRTs and Unc5s mediate neuronal subtype recognition in the retina. Elife 4; 1–27.
• 47   Yoshida,  M.,  Kikutani,  T.,  Tsuga,  K.,  Utanohara,  Y.,  Hayashi,  R. and  Akagawa,  Y. (2006) Decreased tongue pressure reflects symptom of dysphagia. Dysphagia 21; 61–65.
• 48  Zang, Y., Chaudhari, K. and Bashaw, G. J. (2021) New insights into the molecular mechanisms of axon guidance receptor regulation and signaling. in “Current Topics in Developmental Biology”, pp. 147–196, Academic Press Inc., 142; 147–196.
• 49   Zhao,  G.,  Li,  Q.,  Wang,  I.-M. M.,  Liu,  X.,  Fang,  X. and  Zhang,  X. D. (2015) An effective analytic method for detecting tissue-specific genes in RNA-seq experiments. Pharmacogenomics 16; 1769–1779.
• 50   Zimmer,  G.,  Schanuel,  S. M.,  Bürger,  S.,  Weth,  F.,  Steinecke,  A.,  Bolz,  J., et al. (2010) Chondroitin sulfate acts in concert with semaphorin 3A to guide tangential migration of cortical interneurons in the ventral telencephalon. Cereb. Cortex. 20; 2411–2422.