Autism-linked NLGN3 is a key regulator of gonadotropin-releasing hormone deficiency

Gonadotropin-releasing hormone (GnRH, encoded by GNRH1) is the master hormone regulating the hypothalamic–pituitary–gonadal (HPG) reproductive axis, and its pulsatile secretion is crucial for puberty onset, sexual maturation and fertility (Herbison, 2016). GnRH is produced by a small number of hypothalamic neuroendocrine neurons called GnRH neurons, which, during embryonic development, originate in the nasal placode and migrate along vomeronasal nerves to reach the hypothalamus (Oleari et al., 2019).

Disruption in GnRH neuron development or hypothalamic function are the leading causes of genetic reproductive disorders such as hypogonadotropic hypogonadism (HH) and Kallmann syndrome, which are characterized by absent or delayed puberty, owing to GnRH deficiency (GD) (Boehm et al., 2015). The known genes causing GD account for only 50% of cases (Boehm et al., 2015; Oleari et al., 2021), supporting the need for studying the genetic signatures of GnRH neurons as a novel strategy to expedite etiological discovery in the remaining cases. However, the study of the GnRH neuronal system is hampered by the difficulty in obtaining primary GnRH neurons, which are small in number and lack specific markers. To overcome these issues, alternative tools, such as immortalized GnRH neuron cell lines (Maggi et al., 2000) and reporter rodent and zebrafish lines (Messina et al., 2016; Kato et al., 2003; Abraham et al., 2009), have been adopted.

Recently, RNA sequencing has been applied to obtain the transcriptome of induced pluripotent stem cell (iPSC)-derived human GnRH neuron progenitors and early postmitotic GnRH neurons (Lund et al., 2020; Keen et al., 2021; Wang et al., 2022). Yet, no reports of gene expression profiles of immortalized or primary GnRH neurons during the entire developmental process are available.

Here, we obtained, for the first time, the transcriptomic profiles of rodent immortalized and primary embryonic GnRH neurons. Further, by combining filtering strategies with exome sequencing from human patients and in vitro functional experiments, we identified neuroligin 3 (NLGN3) as a new GD-disease candidate gene. NLGN3 belongs to the neuroligin family, a class of postsynaptic cell adhesion molecules that regulate synapse organization (Uchigashima et al., 2021) and dendritic outgrowth (Xu et al., 2019), and NLGN3 missense variants have previously been associated with autistic spectrum disorder (ASD) (Nguyen et al., 2020). In this work, we have described two patients carrying novel nonsense NLGN3 variants and presenting with clinical phenotypes of GD and ASD, therefore providing the first genetic correlation between these two neurodevelopmental disorders, which was previously only implied by a registry-based association study (Oleari et al., 2021). Moreover, we have demonstrated that NLGN3 is important for GnRH neuronal development, with deficiency resulting in defective neuritogenesis in a cellular model.

We first obtained the transcriptomic profiles of GN11 and GT1-7 cells, which represent well-characterized immortalized GnRH neuron cell lines, by using an Affymetrix GeneChip array (n=3; GSE174902; Fig. 1A). These neurons, which derive from mouse tumors localized in the olfactory bulb and in the hypothalamus, respectively, are characterized by different cellular properties, namely migration for GN11 cells and GnRH hormone secretion for post-migratory GT1-7 cells (Maggi et al., 2000). Sample clustering by principal component analysis (PCA) revealed a high degree of similarity within samples of each cell line (Fig. 1B). Among the two cell lines, 1612 differentially expressed probes (|log_FC|>2, adjusted P<0.05; Fig. 1C,D), mapping to 1515 unique annotated genes, were found. Specifically, 703 genes were significantly upregulated in GN11 cells, while GT1-7 cells showed a significant overexpression of 812 genes. Functional enrichment analysis, performed with reString software (Manzini et al., 2021) and based on Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotations, revealed significantly enriched pathways [false discovery rate (FDR)<0.01] that were consistent with the different biological properties of the two cell lines (Fig. 1E). In particular, genes upregulated in GN11 cells significantly populated pathways related to migration (Fig. S1A). Conversely, the most represented pathways in GT1-7 cells belonged to neuronal maturation processes (e.g. axon development and synapse assembly) (Fig. S1B). Similar results were obtained with the functional enrichment analysis of Gene Ontology (GO) biological processes (BP) (FDR<10⁻⁵; Fig. S1C). Overall, these analyses identified that a high number of maturational stage-specific genes enriched biological pathways relevant to and consistent with the different phases of GnRH neuron development, such as migration at early phases (GN11) and establishment of neuron projections and synapse development at later phases (GT-7) (Forni and Wray, 2015; Cho et al., 2019).

To refine the list of candidate genes, we leveraged an integrated approach based on functional enrichment analysis, transcriptomic profiles of primary GnRH neurons and prioritization bioinformatic tools (Fig. 2A).

We first selected genes belonging to significantly enriched pathways (FDR<0.01 for KEGG pathways and FDR<10⁻⁵ for GO BP), which related to cytoskeleton dynamics, migration and axon development for each cell line (Tables S1 and S2).

Then, to confirm the biological relevance of selected genes in vivo, we integrated these gene lists with those obtained by transcriptomic profiling of primary GnRH neurons isolated from fluorescence-activated cell sorting (FACS)-sorted Gnrh1-GFP rat embryos (Kato et al., 2003) at embryonic day (E)14, E17 and E20 (n=1 per each stage; GSE174896; Fig. S2A,B). Owing to technical constraints, RNA from primary neurons from different embryos was pooled and run as a single replicate, for each stage. By examining GFP⁺ cell transcriptomes using the E14 stage as a common baseline, we identified expression changes in 4844 detected probes, mapping to 2739 unique annotated genes (|log_FC|>1; Fig. S2C,D), including genes previously implicated in GnRH neuron development by direct experimental evidence (e.g. Sirt1, Reln, Cxcr7, Slit2, Lif) (Cho et al., 2019; Cariboni et al., 2007) or found mutated in GD patients (e.g. SEMA3A, SEMA7A, FGFR1, NRP1) (Boehm et al., 2015). Interestingly, each developmental stage was characterized by a unique transcriptomic signature (Fig. S3) that exhibited a high level of overlap with the ones observed in immortalized neurons (Fig. S4). The two datasets were then integrated to obtain two gene set clusters, containing genes from selected enriched pathways of GN11 and GT1-7 cells that exhibited a decreasing (‘early’ genes) or increasing (‘late’ genes) expression trend along with embryonic age in primary GFP⁺ cells (Fig. 2B,C).

Last, to identify, among filtered genes, those that could be potentially implicated in GD pathogenesis, ToppGene software (Chen et al., 2009) was applied to rank candidates based on functional or structural similarity with a set of known GD causative ‘input’ genes (Boehm et al., 2015; Maione et al., 2018) (Table S3). The top output genes for each list (early or late) with an overall P-value less than 0.01 were selected, thus obtaining seven early and 22 late genes (Table 1). The expression levels of two representative early and late genes were validated in GN11 and GT1-7 cells by reverse transcription quantitative PCR (RT-qPCR) (Fig. S5).

In summary, the microarray analyses combined with candidate gene-filtering strategies facilitated the identification of 29 candidate genes that, with the exception of SEMA3A, have not been previously linked with GnRH neuron biology or GD.

To validate the relevance of the top 29 genes in human GD pathogenesis, the presence of potentially deleterious variants for each gene was interrogated in exomes from our cohort of GD patients (n=47) (Saengkaew et al., 2021). Only variants that met the American College of Medical Genetics (ACMG) criteria for pathogenicity, likely pathogenicity or variants of uncertain significance were retained in the analysis. Overall, we found three variants with Combined Annotation-Dependent Depletion (CADD) >25 and minor allele frequency <1% (Table S4). Of these, only one predicted loss-of-function (LoF) variant was identified, in the X-linked gene NLGN3 (Gene ID 5441). A male patient (Case 1) with partial GD [with puberty that had initiated and then arrested (Ferreira et al., 2013)] and ASD traits was found to carry a stop-gain variant in NLGN3 (Fig. 3A). By interrogating GeneMatcher (Sobreira et al., 2015), we found an additional NLGN3 LoF variant in an unrelated patient with phenotypic features of GD and ASD with intellectual disability (Case 2; Fig. 3A).

Case 1 was found to be hemizygous for a novel stop-gain variant in NLGN3 (NM_181303.1: c.366G>A, p.W122*; VCV002443311.1) (Fig. 3A,B). He was born to healthy non-consanguineous parents, and his mother was found to be heterozygous for the same variant, as confirmed by Sanger sequencing (Fig. 3A). Clinically, the proband presented in early adulthood with a picture of pubertal arrest with biochemical evidence of GD (Table 2). The proband reported late onset of puberty and had achieved testes volume of 10 ml by the age of first review. He was commenced on treatment with testosterone esters, which were gradually increased up to a maximal dose of 250 mg every 4 weeks, but testes volume had increased only to 12 ml at the age of more than 20 years, consistent with partial GD phenotype (Boehm et al., 2015). In addition, he had associated phenotypic features including obesity [body mass index (BMI), 36; height, 178.5 cm], depression, anxiety, and social and communication difficulties consistent with ASD. Inhibin B concentrations were in the low-normal range. Mildly elevated thyroid-stimulating hormone (TSH) concentrations with low-normal free thyroxine (T4) were reported, likely secondary to obesity (Table 2). Magnetic resonance imaging scan showed a normal appearance of the pituitary gland, pituitary stalk and hypothalamus (Fig. S6). There was no family history of GD, but his mother, who carried the same variant in NLGN3, was also obese.

Case 2 (NM_181303.1: c.163C>T; p.R55*; VCV002442803.1) is a pre-pubertal boy who presented with stereotypia, learning difficulties and ASD features. In addition, he was noted to have micropenis (untreated penile length 2.5 cm in mid-childhood) and absent erections, which are known predictive signs of HH in pediatric males (Boehm et al., 2015; Swee and Quinton, 2019). At the age of 7-12 years, he was obese (BMI Z-score, +2.1; height, 148 cm) with small testes volumes (1.5 ml bilaterally). His mother also had history of learning difficulties and stereotypic behaviors, but there was no family history of GD. He had appropriately low gonadotropins and testosterone for age when first tested in mid-childhood, with normal inhibin B and low anti-Müllerian hormone (AMH) (Table 2).

In both cases, exome sequencing of the probands excluded deleterious variants in known GD genes, and no other variants of interest were identified in relevant genes, e.g. related to developmental delay. Sanger sequencing of the probands' mothers confirmed them to be the heterozygous carrier in each family, consistent with an X-linked recessive inheritance pattern. DNA from other Case 2 family members was not available.

Interestingly, both variants cause premature stop codons within the NLGN3 extracellular domain, resulting in truncated proteins that are likely to be dysfunctional. Both variants are novel, not found in the gnomAD database (v2.1.1, accessed 6 February 2022), and, notably, the NLGN3 gene is highly intolerant to protein-truncating changes (pLi, 0.98). According to the ACMG (Richards et al., 2015), these variants are both classified as pathogenic (Table 3). Bioinformatic tools predicted these variants to be damaging, deleterious and disease causing, with CADD score >35 and DANN score >0.99 (Fig. 3B and Table 3). Further, multi-species protein alignment combined with GERP++ software showed that identified variants affect highly conserved residues (Fig. 3B).

Our transcriptomic analyses revealed that the X-linked NLGN3 gene was upregulated in primary and immortalized GnRH neurons during development. Thus, to validate its developmentally regulated expression, we performed the following experiments. First, we confirmed higher NLGN3 levels in maturing GT1-7 versus immature GN11 cells (Fig. 4A,B), by RT-qPCR and immunocytochemistry, using a previously validated anti-NLGN3 antibody. Then, we analyzed NLGN3 expression in E14.5 mouse sections and confirmed colocalization of NLGN3 protein in GnRH neurons that had settled in the maturing hypothalamus (Fig. 4C).

To functionally validate the predicted pathogenicity of the identified NLGN3 variants, we overexpressed human hemagglutinin (HA)-tagged wild type (WT) (Chih et al., 2004), R55* or W122* NLGN3 in cell culture models. Initially, protein synthesis and secretion of WT versus mutated NLGN3 proteins were evaluated in COS7 cells, by immunoblotting analysis on lysates and conditioned media, respectively. As expected, a 110 kDa band corresponding to the full-length protein was observed in lysates of COS7 cells transfected with NLGN3 WT (Fig. 5A) (Bemben et al., 2019). In contrast, R55* and W122* variants led to the formation of prematurely truncated proteins of expected molecular masses (7 and 14 kDa, respectively) (Fig. 5A). In addition, a shorter form (90 kDa) of NLGN3 WT, which likely represents the ectodomain released by proteolytic shedding (Bemben et al., 2019), was detected in conditioned media. We could not detect equivalent shed or secreted forms for the NLGN3 mutant proteins (Fig. 5A). To study localization of these mutant proteins, we performed similar overexpression experiments and immunocytochemical analysis in GN11 cells, which show low levels of endogenous NLGN3 (Fig. 4A,B). Immunofluorescence staining for HA in GN11 cells revealed that mutant proteins were only detected in the endoplasmic reticulum (ER), selectively labeled with mEmerald-ER-3 vector (Fig. 5B). By comparison, the WT protein was efficiently transported to the plasma membrane. Of note, we observed only very few HA⁺ cells when transfected with R55* NLGN3, strongly supporting rapid NLGN3 protein degradation and consequent lack of detection, as previously described for other NLGN3 mutant proteins (Quartier et al., 2019). Altogether, these observations strongly support a LoF effect of the two newly identified variants.

Because NLGN3 promotes neuritogenesis in human neural progenitor cells (Gatford et al., 2022), we studied the effect of WT or mutant NLGN3 proteins on GnRH neuron morphology and neurite outgrowth by overexpression experiments in GN11 cells, as previously described (Bouilly et al., 2018). NLGN3 WT promoted the formation of several protrusions, with cells displaying an increased spreading compared to GFP-transfected controls (Fig. 5C). By contrast, the two mutants (R55* and W122*) were unable to induce such changes (Fig. 5C), with transfected cells displaying a significantly lower cell perimeter, although total area and cell polarity (i.e. aspect ratio) were not affected (Fig. 5D). Further, cells overexpressing NLGN3 mutants exhibited a significantly reduced cell shape complexity (i.e. complexity index) compared to NLGN3 WT, revealing a strong defect in cell protrusion formation (Fig. 5D). Together, these results confirmed that the identified NLGN3 variants are LoF and failed to induce cytoskeletal and membrane rearrangements typical of neuritogenesis.

In this work, we have interrogated the transcriptomic profiles of rodent immortalized and developing primary GnRH neurons and developed a new resource for understanding the key developmental regulation of GnRH neurons. Further, we have identified NLGN3 as modulator of GnRH neuron neuritogenesis and novel NLGN3 LoF variants in two unrelated patients with GD and ASD traits.

Previous works have generated transcriptomes from human iPSC- or embryonic stem cell-derived GnRH-like cells by RNA sequencing (Lund et al., 2020; Keen et al., 2021) and, very recently, at single-cell level (Wang et al., 2022). Although caution should be exercised in transcriptomic comparative analyses between different models and sequencing techniques, our transcriptomic data highlighted the presence of consistent expression trends in the majority of our 29 candidate genes across GnRH neuron clusters recently identified by a single-cell RNA-sequencing approach in human iPSC-derived GnRH neurons (Wang et al., 2022). Further, candidate genes generated from our analyses have been experimentally associated with GnRH neuron development (e.g. Bmp4, Otx2, Plxnc1) (Forni et al., 2013; Parkash et al., 2015; Diaczok et al., 2011) or with GD pathogenesis (e.g. SEMA3A) (Oleari et al., 2021), confirming the ‘bona-fide’ nature of immortalized and primary GnRH neurons as complementary models for studying physiological and pathological mechanisms of GnRH neuron development.

As an exemplar, our integrated omic analysis allowed us to identify NLGN3 as a new candidate GD gene. Specifically, the identified variant in our Case 1 proband induces a premature stop codon in NLGN3, a gene that has been linked so far to ASD (Nguyen et al., 2020; Quartier et al., 2019; Jamain et al., 2003; Redin et al., 2014) but not to GD. Notably, the second patient identified in this study (Case 2) is also hemizygous for a nonsense NLGN3 variant and, although prepubertal, displayed typical red flag signs of hypogonadism, including small testes volume and micropenis (Boehm et al., 2015; Swee and Quinton, 2019). In addition, these two probands had behavioral difficulties including ASD features or neurodevelopmental disorder (NDD), supporting the existence of a shared pathophysiological link between GD and other developmental disorders, as hypothesized by a recent phenotypic population study showing an increased risk of NDDs in patients with GD (Ohlsson Gotby et al., 2019). This implies that NLGN3 is essential for the development or function of GnRH neurons, besides its well-known role in brain neurotransmission, which is thought to be dysfunctional in ASD (Nguyen et al., 2020). Although its expression has been reported in the developing telencephalon of chick and zebrafish embryos in territories relevant to GnRH neuron development (Davey et al., 2010; Paraoanu et al., 2006), the role of NLGN3 in this system is not yet known.

In this context, because NLGN3 is developmentally upregulated both in immortalized and primary GnRH neurons and its overexpression in GN11 cells promotes neuritogenesis, NLGN3 may be required to promote neurite extension during the final phases of GnRH neuron development in vivo. Consequently, NLGN3 loss might have an impact on this process, thus causing GD. Interestingly, compound double Nlgn1/3- or Nlgn2/3-null mice display a reduced reproductive rate (Varoqueaux et al., 2006), in addition to behavioral phenotypes related to ASD (Varghese et al., 2017), strongly supporting a role for NLGN3 in the reproductive axis. Yet, future detailed phenotypic analyses of GnRH neuron development in NLGNs single and compound null mice will be necessary to define the exact function of NLGN3 in these systems.

Our study also links, for the first time, nonsense NLGN3 variants to a combined phenotype of ASD and GD. Most ASD patients with NLGN3 mutations carry missense variants, with only two previous reports of nonsense variants without functional characterization or phenotypic correlation (Nguyen et al., 2020). Further, although we cannot exclude that the identified nonsense NLGN3 variants undergo degradation in vivo due to the nonsense-mediated decay mechanism of RNA surveillance, our data still provide evidence that if this mechanism does not occur, remaining NLGN3 mutant proteins are non-functional, and their expression leads to a relevant cellular phenotype.

Although the association of truncating NLGN3 variants with the co-existence of GD with ASD or NDD is being supported so far only by the two cases reported in our study, it will be interesting to expand genotype–phenotype correlation to other cases and to understand the molecular mechanisms underlying this association. Interestingly, growing evidence supports the existence of common genetic determinants between NDDs and GD, as shown by a 2022 study linking SOX11 mutations to a novel NDD with HH (Al-Jawahiri et al., 2022). It is known that GnRH neurons also send projections to extrahypothalamic areas, including brain regions involved in intellectual functions (Casoni et al., 2016). A very recent study showed that GnRH treatment is effective in restoring cognitive performance in patients with Down syndrome (Manfredi-Lozano et al., 2022). Thus, we speculate that NLGN3 might also have a role in GnRH neuron neuritogenesis in non-hypothalamic regions, thus contributing to the neurodevelopmental phenotype directly via GnRH neurons as well as via non-GnRH networks.

Overall, our study, by exploiting the power of omic analyses, has provided a novel experimental resource that will be valuable for the identification of the elusive mechanisms underlying GnRH neuron biology and related diseases. Further, our findings provide the first proof of principle that the transcriptomic analysis of immortalized and primary GnRH neurons, combined with human genetics, in silico tools and in vitro models, may be applied to identify novel GD candidate genes. In the future, mutational screening of our top 29 genes in other patient cohorts could identify novel genetic determinants of GD; of note, PLXNC1 and CLSTN2 variants of interest in single patients from our cohort were found. Lastly, because our studies imply the existence of a shared genetic pathway between ASD and GD, it will be important for clinicians to be aware of the potential for pubertal and reproductive disorders in children with ASD/NDD to identify those who might require hormonal therapy intervention and, vice versa, to assess the psychological and social aspects of GD-affected children reporting ASD or developmental delay.

Gnrh1-GFP rat embryos (Kato et al., 2003) were used in Wistar background to isolate GnRH-GFP neurons. WT C57/Bl6 mouse embryos were used for expression studies (Italian Ministry of Health, license 5247B.N.QPE). To obtain embryos of defined gestational stages, rats and mice were mated in the evening, and the morning of vaginal plug formation was counted as E0.5.

GN11, GT1-7 cells (Maggi et al., 2000) and COS7 cells were grown as a monolayer at 37°C in a humidified CO₂ incubator in complete Dulbecco's modified Eagle medium (Euroclone) supplemented with 10% fetal bovine serum (FBS; Invitrogen). Subconfluent cells were harvested by trypsinization and cultured in 57 cm² dishes. Cells within six passages were used for all experiments.

Explants from the nasal area (E14), both nasal and basal forebrain areas (E17) and from basal forebrain area (E20) were microdissected from heads of Gnrh1-GFP embryos; nose and forebrain cells were dissociated by incubation in 0.05% trypsin with 100 μg/ml DNase I in neurobasal medium (Invitrogen) at 37°C for 15 min. Trypsinization was quenched by addition of neurobasal medium containing 10% heat-inactivated FBS (Invitrogen) at 37°C for 5 min. Cells were washed three times in neurobasal medium (without FBS) to remove serum before FACS and resuspended in neurobasal medium without Phenol Red (Invitrogen) containing L-glutamine (Invitrogen) and B-27 supplement (1:50; Invitrogen). Dissociated cells from eight to ten embryos for each age were pooled for each FACS. FACS was performed by the Wolfson Scientific Support Services (University College London, London, UK) by using a MoFlo Sorter (Dako). A non-green embryo was used as a control for fluorescence. Cells were excited by using a 488-nm Argon laser and detected by using a 530/40 (FL1) bandpass filter. A cell purity of 95-98.5% was obtained for each sort. Sorted (GFP⁺) and unsorted (GFP⁻) cells were directly collected in lysis buffer (Qiagen) and used for RNA extraction.

For GN11 and GT1-7 cells, total RNA was isolated by acid guanidinium thiocyanate-phenol-chloroform extraction (Trizol, Invitrogen) followed by a Qiagen RNeasy kit clean-up procedure (Qiagen). Total RNA from GFP⁺ and GFP⁻ FACS-purified cells was extracted immediately after collection by using a Qiagen RNeasy Plus kit (Qiagen). RNA integrity was verified with an RNA 6000 Nano and Pico kit, and only high-quality RNAs, with RNA integrity number greater than 7 were used. A WT Expression Kit (Invitrogen) and Ovation Pico WTA System V2 (Nugen) were respectively used to prepare sense-strand cDNA. Biotin labeling was performed with a WT Terminal Labeling kit Affymetrix (Invitrogen) and Encore Biotin Module kit (Nugen). Labeled cDNA was hybridized (45°C for 17 h) to the Affymetrix Mouse Gene 1.0 ST and Rat Gene 1.0 ST Gene chips according to the manufacturer's protocol.

Microarray data analysis, including quality controls, normalization and gene filtering was conducted using AMDA software (Pelizzola et al., 2006). The identification of differentially expressed genes was addressed using linear modeling approach and empirical Bayes methods (Smyth, 2004) together with FDR correction of the P-value (Benjamini–Hochberg); the adjusted P-value that has been selected in the analysis is <0.05 and log_FC>2. GO analyses were performed as previously described with reString software (Manzini et al., 2021). Briefly, hierarchical clustering was performed with Euclidean metric on log-transformed ratios of probes fluorescence intensity to average probe intensity (Virtanen et al., 2019 preprint). PCA was performed with Scikit-learn (Pedregosa et al., 2011). Data visualization was performed with SciPy (Virtanen et al., 2019 preprint), matplotlib (Hunter, 2007) and seaborn (Waskom et al., 2022) libraries for the Python programming language.

For case 1, exome-sequencing data from 47 patients from a UK GD cohort (CPMS ID 30730) were analyzed. Exome sequencing was performed on DNA extracted from peripheral blood leukocytes, using an Agilent V5 platform and Illumina HiSeq 2000 sequencing. The exome sequences were aligned to the UCSC hg19 reference genome using Burrows-Wheeler Aligner software (BWA-MEM; bwa-0.7.12). Picard tools (picard-tools-1.119) was used to sort alignments and mark PCR duplicates. The genome analysis toolkit (GATK-3.4-46) was used to realign around indels and recalibrate quality scores using dbSNP, Mills and 1000 genomes as reference resources. Variant calling and joint genotyping using pedigree information was performed using HaplotypeCaller in genomic variant call format (GVCF) mode from the genome analysis toolkit. The resulting variants were filtered using the variant quality score recalibration function from GATK. An analysis of the called variants was performed using Ingenuity Variant Analysis (Qiagen). Filtering for potential causal variants was carried out using filters for quality control, allele frequency and predicted functional annotation. Potentially pathogenic variants in candidate genes were verified by Sanger sequencing. Libraries of genomic DNA samples were prepared using a Sureselect Human All Exon v5 kit (Agilent Technologies) and were sequenced on a HiSeq instrument (Illumina) according to the manufacturer's recommendations for paired-end 76-bp reads. The bioinformatics pipeline, alignment processes and quality procedures were as previously described (Yang et al., 2014). Version 3.4-46 of the Genome Analysis Toolkit was used for this study.

For case 2, exome sequencing was realized as previously published (Nambot et al., 2018). Briefly, libraries of genomic DNA samples were prepared using a Human Core Exome kit (Twist Biosciences) and were sequenced on a NovaSeq 6000 instrument (Illumina) according to the manufacturer's recommendations for paired-end 151-bp reads. A mean depth of 86.96× was reached, and 97.2% of the refseq exons were covered at least by ten reads. Variants were identified using a computational platform of the Fédération Hospitalo-Universitaire (FHU) TRANSLAD, hosted by the University of Burgundy Computing Cluster. Raw data quality was evaluated by FastQC software (v0.11.4). Reads were aligned to the GRCh37/hg19 human genome reference sequence using the Burrows-Wheeler Aligner (v0.7.15) (Li and Durbin, 2009). Aligned read data underwent the following steps: (1) duplicate paired-end reads were removed by Picard software (v2.4.1), and (2) base quality score recalibration was done by the Genome Analysis Toolkit (GATK v3.8) Base recalibrator. Using GATK Haplotype Caller, single-nucleotide variants with a quality score >30 and an alignment quality score >20 were annotated with SNPEff (v4.3) (Cingolani et al., 2012). Rare variants were identified by focusing on nonsynonymous changes present at a frequency less than 1% in the GnomAD database. Copy number variants were detected using xHMM (v1.0), annotated using in-house python scripts and then filtered regarding their frequency in public databases (Database of Genomic Variants, International Standards for Cytogenomic Arrays, Deciphering Developmental Disorders).

To introduce the c.163C>T and c.366G>A variants into human NLGN3 gene, the WT human HA-tagged NLGN3 expression vector (Addgene, 59318) was mutagenized using a QuickChange Lightning Site-Directed Mutagenesis Kit (Agilent Technologies) and specific oligonucleotides for NLGN3 R55* (fw, 5′-GGCAGTGGTACTCAGGCACCCCTTAGC-3′; rev, 5′-GCTAAGGGGTGCCTGAGTACCACTGCC-3′) and W122* (fw, 5′-GTCATGCTGCCGGTCTGATTCACTGCCAACTTGGATATCG-3′; rev, 5′-TCAGACCGGCAGCATGACTTCGGGCACAGCTGTGTGGATG-3′). To visualize the ER, the mEmerald-ER-3 vector (a gift from Prof. Diego De Stefani; Addgene, 54082) was co-transfected. Expression vectors were transiently transfected (for 24 h) into GN11 cells using Lipofectamine 3000 (Invitrogen).

Cells were lysed in 150 mM NaCl, 50 mM Tris-HCl (pH 7.4) and 1% Triton X-100, supplemented with protease and phosphatase inhibitors (Roche). Then, 20 μg of proteins were transferred to nitrocellulose membranes (Bio-Rad), immunoblotted with rabbit anti-HA (1:1000; Cell Signaling Technology, 3724) and rabbit anti-GAPDH (1:5000; Cell Signaling Technology, 2118), followed by horseradish peroxidase-conjugated anti-rabbit antibody (1:10,000; Santa Cruz Biotechnology, sc-2301).

Paraformaldehyde (PFA)-fixed GN11 cells were incubated with PBS containing 10% normal goat serum and 0.1% Triton X-100. Rabbit anti-HA (1:800; Cell Signaling Technology, 3724) and chicken anti-GFP (1:1000; Abcam, ab13970) were used as primary antibodies. Secondary antibodies used were Alexa Fluor 488-conjugated donkey anti-rabbit and anti-chicken and Cy3-conjugated donkey anti-rabbit Fab fragments (1:200; Jackson ImmunoResearch). To detect F-actin, cells were stained with TRITC-conjugated phalloidin (1:400; Sigma-Aldrich, P1951) for 30 min at 37°C (Cannarella et al., 2021). Nuclei were counterstained with DAPI (1:10,000; Sigma-Aldrich, D9542). For immunoperoxidase labeling, cells were incubated with hydrogen peroxide to quench endogenous peroxidase activity before incubation with rabbit anti-NLGN3 (1:200; Budreck and Scheiffele, 2007), followed by biotinylated goat anti-rabbit antibody (1:400; Vector Laboratories, BA-1000), and then developed with an ABC kit (Vector Laboratories, PK6100) and DAB (Sigma-Aldrich, D4293).

Total RNA from GN11, GT1-7 and FACS-sorted cells was retrotranscribed into cDNA as previously described (Busnelli et al., 2020). Briefly, 1 μg of total RNA was reverse transcribed with random hexamers and MultiScribe reverse transcriptase (Applied Biosystems) following the manufacturer's Ginstructions. The expression level of murine Nlgn3 (fw, 5′-GAAGATGGATCCGGCGCTAA-3′; rev, 5′-ACGATGACGTTGCCGTAACT-3′), Sema3a (fw, 5′-CGTCTTCCGGGAACCAACAA-3′; rev, 5′-TGCACAGGCTTTGCCATAGA-3′), Snai2 (fw, 5′-GAACTGGACACACACACAGTTATT-3′; rev, 5′-TGCCGACGATGTCCATACAG-3′), L1cam (fw, 5′-GCTCCTCATCCTGCTCATCC-3′; rev, 5′-TCTCCAGGGACCTGTACTCG-3′) and Sema3c (fw, 5′-GAACCCATGTTTGTGGACGC-3′; rev, 5′-CCACCAGTGTCATTAGGGCA-3′) genes was quantified by quantitative PCR (qPCR) on a Bio-Rad CFX Connect thermal cycler with Luna Universal qPCR Master Mix (NEB) in 10 μl reactions, with final concentration of 0.25 μM for each primer. The cycling conditions were 95°C for 1 min, followed by 40 cycles of 15 s at 95°C, 30 s at 60°C and 30 s at 72°C. A final melting curve analysis assured the authenticity of the target product. Triplicate samples were run in all reactions; first-strand DNA synthesis reactions without reverse transcriptase were used as controls. The ΔCq value and the ΔΔCq were calculated relative to control samples using quantification cycle (Cq) threshold values that were normalized to the reference gene, Gapdh (fw, 5′-CATCCCAGAGCTGAACG-3′; rev, 5′-CTGGTCCTCAGTGTAGCC-3′).

WT C57/Bl6 embryos at E14.5 were fixed for 3 h in 4% PFA and then cryoprotected overnight in 30% sucrose, embedded in OCT and cryosectioned for immunostaining. PFA-fixed tissue sections (25 μm) or cells were incubated with serum-free protein block (Dako, X0909). For immunofluorescence staining, goat anti-NLGN3 (1:25; Santa-Cruz Biotechnology, sc-14091), rabbit anti-GnRH (1:400; Immunostar, 20075), rabbit anti-NLGN3 (1:100; Budreck and Scheiffele, 2007) and goat anti-PLXND1 (1:200; R&D Systems, AF4160) (Cariboni et al., 2015) were used as primary antibodies. Secondary antibodies used were Alexa Fluor 488-conjugated donkey anti-rabbit and Cy3-conjugated donkey anti-goat Fab fragments (1:200; Jackson ImmunoResearch). Nuclei were counterstained with DAPI (1:10,000; Sigma-Aldrich, D9542).

Cells were examined with a Zeiss LSM 900 confocal laser scanning microscope equipped with a Zeiss Axiocam 305 color and using a Plan Apochromat 40×1.3 Pil DIC (VIS-IR M27) oil immersion objective (Zeiss). DAPI, Alexa Fluor 488 and Cy3 were excited at 353, 493 and 548 nm and observed at 400-500, 500-540 and 540- 700 nm, respectively. We captured 1024×1024 pixel images in a stepwise fashion over a defined z-focus range corresponding to all visible fluorescence within the cell. Maximum projections of the z-stack with 0.25 μm optical section were obtained post-acquisition by using Zeiss ZEN System. Adobe Photoshop CS6 software was used to prepare the presented images.

Morphological analyses were performed on confocal images using ImageJ (v1.52a) as previously described (Bouilly et al., 2018) and summarized in Fig. 4D.

Statistical tests employed are outlined in figure legends and were conducted when the experiment had been performed a minimum of three times on a minimum of three individual samples. Data are presented as mean±s.d., and results were considered significant with a P-value less than 0.05 (GraphPad Prism 7.0).

Ethical approval for Case 1 was granted by the London-Chelsea National Research Ethics Service committee (13/LO/0257) and the UK NHS Health Research Authority (IRAS 95781). Ethical approval for Case 2 was granted by Ethics Committees of Centre Hospitalier Universitaire de Caen and Centre Hospitalier Universitaire Dijon Bourgogne. All participants provided written informed consent prior to study participation. The study was conducted in accordance with the guidelines of the Declaration of Helsinki. All individual-level data, including clinical data, were de-identified. Participants or their legal representatives gave consent to the publication of the results of this research work in the present study. The animal work was approved by the University of Milan Animal Welfare Body and by the Italian Ministry of Health to A.C. and was conducted in accordance with the EU Directive 2010/63/EU.

We thank Prof. Maria Foti and Genopolis gene-sequencing facility (University of Milano-Bicocca) for their help with Affymetrix microarrays and analyses; Dr Thomas Adejumo for help with FACS at the University College London Wolfson Institute for Biomedical Research; and Prof. John Parnavelas for laboratory resources. We also thank the imaging facility NOLIMITS at the University of Milan for help with confocal imaging.