A pair of 4 week old male and female gerbils (Meriones unguiculatus; RjTub:MON) was purchased from Janvier Labs, France. Upon arrival animals were euthanized and brains were dissected and collected in TRIzol^® reagent (Invitrogen, Carlsbad, CA, United States) for RNA extraction. Liver biopsies were collected in tissue lysis buffer for DNA extraction, and subsequent DNA isolation was performed using the phenol-chloroform extraction method (Maxeiner et al., 2019). All animal handling followed standards set by and protocols approved by local authorities as well as Saarland University, Homburg, Germany.

To determine the transcriptional start site of neuroligin-4 in gerbils, RNA from the female gerbil was analyzed by 5′ RACE (5′ rapid amplification of cDNA ends; SMARTer^® RACE 5′/3′ Kit, Takara Bio Europe SAS, Saint-Germain-en-Laye, France) following the supplier’s instructions for the amplification of 5′ cDNA ends. RNA was transcribed into cDNA and an artificial DNA sequence was added. The cDNA then served as template for gene specific primers and a universal primer provided by the kit. Amplification of gene products was performed using SeqAmp DNA Polymerase from the kit or Q5 High-Fidelity 2x Master Mix (New England Biolabs GmbH, Frankfurt, Germany). Amplification products of the 5′ RACE experiment and other PCR amplicons used for genomic sequencing were separated on agarose gels and relevant bands were extracted using the Zymo Gel extraction kit (Zymo Research Europe GmbH, Freiburg, Germany). Amplicons were cloned using the In-Fusion cloning kit (Takara Bio Europe SAS). Amplification of plasmid DNA and bacterial culture was performed according to standard procedures, and plasmid DNA for subsequent applications was prepared using resin-based DNA purification systems (Zymo Research). Nucleic acid quantification was performed using the NanoDrop One spectrophotometer (VWR International). Gene specific oligonucleotides were purchased from IDT and are listed in Supplementary Data Sheet 1. Sequencing reactions were performed at the AGCT-Lab sequencing core facility of the Max Planck Institute for Multidisciplinary Sciences, Göttingen, Germany.

Oligonucleotides were designed based on sequence homology between gerbil NLGN4X and NLGN4Y (MX18723/MX18724; cf. Supplementary Data Sheet 1) to amplify DNA in PCR reactions from the gene pair in a single run. PCR reactions were set up using 5 ng of genomic template DNA and Q5 High-Fidelity 2x Master Mix (New England Biolabs) with the following cycle conditions: 98°C for 30 s, 35x cycle (98°C for 10 s, 65°C for 20 s, 72°C for 30 s), and final extension at 72°C for 1 min. All PCR reactions were performed on a T100 Thermal Cycler (Bio-Rad Laboratories, Feldkirchen, Germany). Amplification products were separated on an agarose gel as published previously (Maxeiner et al., 2019).

All bioinformatic analyses were performed essentially as described previously (Maxeiner et al., 2020). Searches for sequence similarities were done using the BlastN suite (NCBI),¹ with default settings for “blast” and “megablast” within specified nucleotide collections and/or WGS for certain species or families. Sequences were aligned using Multalin (Corpet, 1988).² Alternatively, Clustal Omega³ was used in conjunction with the use of MEGA X to infer phylogenetic relationships (Kumar et al., 2018). The following web-based interfaces were used to assess additional features of proteins or nucleic acid sequences: determining the cleavage site of signal peptides (Almagro Armenteros et al., 2019),⁴ the GC-content of any given DNA sequence,⁵ the GC3-content and relative codon usage,⁶ a visual representation of the GC-profile,⁷ and the translation of a DNA sequence into a protein sequence.⁸

English and Latin names of different species were taken from the classification system used along the deposited sequences in GenBank (NCBI), which follows the systematics and classification of Mammal Species of the World (Wilson and Reeder, 2005). Statements on the relative abundance of mammal species derive from the list of different species presented by Burgin et al. (2018). The number of rodent families and their evolutionary relationship was taken from D’Elía et al. (2019), who combined the formerly separate suborders myomorpha, castorimorpha, and anomaluromorpha under the umbrella term supermyomorpha.

All sequences were deposited with GenBank (NCBI) under the following accession numbers: MN045396 (NLGN4Y, closes gap in NHTI01010755.1), MN045397 (NLGN4X, connects NHTI0102956.1, and NHTI01041465.1), MN045398 (NLGN4X, connects NHTI01041465.1 with NHTI01023488.1), and MN045399 (NLGN4X, 5′RACE amplicon covering exon1a, 1b and most of exon 1c).

We previously reported changes in neuroligin-4 genes in a subset of rodents, the clade eumuroida, that were caused by an erosion of the pseudoautosomal region (Maxeiner et al., 2020). In parallel to this work, a draft of the genome of the Mongolian gerbil (Meriones unguiculatus) was released, which belongs with rats and mice to the Muridae family (Wilson and Reeder, 2005). This genome release contained several fragments of incompletely annotated “neuroligin-4-like” genes (Supplementary Data Sheet 2), reminiscent of the diversification of neuroligin-4 into X- and Y-chromosome-specific genes seen in humans and horses (Raudsepp and Chowdhary, 2015). The number of sequences and the fact that the genome of a male gerbil had served as the source of genomic DNA led to the hypothesis that both genes might reflect a yet unprecedented case of neuroligin-4 diversification into NLGN4X and NLGN4Y in a rodent species.

To assign potential X- and Y-chromosome-specific copies of neuroligin-4 genes, we isolated genomic DNA from a male and female Mongolian gerbil. We identified a length polymorphism on “exon 6” (second to last coding exon in all neuroligin genes, Bolliger et al., 2008) that was identified applying a comparable PCR approach adopted from previous sex-typing strategies employing neuroligin-4 genes in horses and humans (Supplementary Data Sheet 3; Maxeiner et al., 2019; Zaffalon et al., 2019; Figure 1B). This attempt resulted in a single band for the neuroligin-4 gene information covered by multiple sequence fragments, thus suggesting that this must represent the gerbil NLGN4X gene. Consequently, we used genomic DNA from the female gerbil to fill in the NLGN4X sequence gaps, and the genomic DNA from the male gerbil to complete the NLGN4Y sequence. Sequencing gaps were closed by PCR cloning and subsequent sequencing of the PCR fragments (Figure 1A).

To determine the transcriptional start site, we performed a 5′RACE experiment (5′ rapid amplification of cDNA ends) with RNA from the brain of the female gerbil, resulting in the identification of two additional exons upstream of the exon with the start codon (Figure 1A). Based on this sequence information, we annotated the corresponding NLGN4Y gene by sequence comparison. Both genes, NLGN4X and NLGN4Y, consist of eight exons, respectively (two exons coding for 5′UTR region, six exons coding for the neuroligin-4 ORF). In line with the previously applied nomenclature for neuroligin genes (Bolliger et al., 2008; Maxeiner et al., 2020), we define the ATG-bearing (start codon-bearing) exon as exon 1 and the final exon encoding the transmembrane and C-terminal region as exon 7. Compared to other neuroligin isoforms, a homolog of exon 2 is absent from neuroligin-2 and -4 genes (Bolliger et al., 2008), which needs to be distinguished from a secondary loss of exon 3 absent in neuroligin-4 genes in certain hamster species (Maxeiner et al., 2020). To acknowledge the presence of any non-coding exons upstream of the ATG-carrying exon 1, we propose to further add alphabetical lettering (i.e., exon 1a, 1b, 1c). The exon base counts for gerbil NLGN4X/Y are as follows: exon 1a (5′UTR coding), 219/227; exon 1b (5′UTR coding), 114; exon 1c (5′UTR/protein coding), 799/792; exon 3 (protein coding), 60; exon 4 (protein coding), 153; exon 5 (protein coding), 186; exon 6 (protein coding), 1,063/1,030; exon 7 (protein coding/3′UTR), 925/937. Whereas the overall GC-content of both neuroligin-4 genes is slightly above 50% (NLGN4X, 56.12%; NLGN4Y, 54.11%), the GC-content is substantially increased in all protein-coding exons judged by the GC-profile along all exons and introns (Figure 1A).

The organization of the gerbil neuroligin-4 genes is characterized by a series of features that are also seen in neuroligin-4 genes of the clade eumuroida, a subset of the suborder myomorpha (Maxeiner et al., 2020). These include the accumulation of repetitive sequences in introns and exons, and a general increase of the total GC and GC3 content, i.e., the increase of GC in the third position of synonymously coding codons. Given the presence of bona fide splice donor and acceptor sites as well as the lack of preliminary stop codons or any frameshifts resulting in non-sense mutations, we conclude that both gerbil neuroligin-4 genes are likely transcribed and translated into proteins. Alignment of NLGN4X and NLGN4Y amino acid sequences (Supplementary Figure 1) showed none or rare amino acid differences in exons 1, 4, and 5, and more substantial differences in exons 3, 6, and 7 comparable to our previous findings (Maxeiner et al., 2020). Functionally important regions, such as the neurexin-binding site and the PDZ-domain binding site are identical. Alignment with the respective human and horse NLGN4X and NLGN4Y protein sequences indicates that the gametolog pair in M. unguiculatus likely arose independently (Figure 2).

The observation of neuroligin-4 divergence into gametologs in gerbils prompted us to revisit our previous assessment of the neuroligin-4 phylogeny (Maxeiner et al., 2020). Including new genome releases of rodent, primate, and lagomorph species as well as those of treeshrews and the colugo, we extended our analysis to the entire mammal superorder euarchontoglires (Murphy et al., 2001) and assessed the phylogeny of neuroligin-4 using a total of 66 sequences (Figure 3). Many of these were not fully annotated in the databases and needed to be inferred by sequence comparison (Supplementary Data Sheets 2, 4).

Further, we aimed at including at least one representative of each of the 35 rodent families. Rodents separate into three major branches, the “porcupine-like” hystricomorpha (17 families), the “squirrel-like” sciuromorpha (3 families), and the “mouse-like” supermyomorpha (15 families) (D’Elía et al., 2019), with the latter one being the largest as regards the number of different species. Supermyomorpha represent approximately 75% of all rodents (Burgin et al., 2018; Maxeiner et al., 2021). Primates, treeshrews (Tupaia spec.), the colugo (G. variegatus), all families within the suborder hystricomorpha and sciuromorpha, as well as Nannospalax galili and Castor canadensis from the suborder supermyomorpha display a very close evolutionary relationship as judged by exon sizes and phylogenetic analysis (Figure 3 and Supplementary Data Sheets 2, 4).

The sequences of Arvicola amphibius, Arvicanthis niloticus and gerbils (M. unguiculatus and Rhombomys opimus) confirmed the neuroligin-4 gene changes within the clade eumuroida described above, i.e., an increase in GC and GC3 content, an incorporation of repetitive sequences, and a reduction of the total gene size. Interestingly, though, we confirmed similar changes now also outside of the eumuroida clade (Supplementary Data Sheets 2, 5). This was, for instance, the case for the kangaroo rat, Dipodomys spectabilis, a member of the Heteromyidae family, which is closely related to beavers (Castoridae). The male kangaroo rat genome data also indicate the presence of two “neuroligin-4-like” genes (Supplementary Data Sheet 2). Both are remarkably similar, with only one amino acid difference (Supplementary Figure 2). Thus, the two isoforms are depicted as “a” and “b” variants as the sequences do not provide conclusive information about neighboring genes that would indicate an X- or Y-chromosomal location. In fact, both genes are identical in their 5 prime regions up to the intron following the exon carrying the presumptive start codon, indicating that D. spectabilis neuroligin-4 genes might be in the process of transitioning from pseudoautosomal to gametologous genes. Remarkably, exon 3 is also absent from both neuroligin-4 copies in D. spectabilis, reminiscent of some hamster species (Maxeiner et al., 2020).

Lagomorphs represent an order that is most closely related to rodents. Analysis of their neuroligin-4 sequences revealed changes in the neuroligin-4 coding sequence that appear to have happened independently to those observed within rodents, or more specifically within supermyomorpha. Our analysis of four representatives from three different genera (Oryctolagus, Sylvilagus and Lepus) of lagomorphs reveal that they cluster tightly, but clearly apart from the group of highly similar neuroligin-4 orthologs present in primates and hystricomorpha (Figure 3, for branching statistics see Supplementary Figure 3). Gene and protein sequence analyses indicate that apart from their phylogenetic distance to rodents likewise changes as observed in supermyomorpha have emerged. These include an increase of the GC- and GC3-content, a reduction of the overall gene size, and an accumulation of repetitive sequences. These changes within the coding region are immediately evident by plotting the GC-and GC3-content for each species (Figure 4 and Supplementary Data Sheet 5). All investigated sequences from primates and hystricomorpha cluster in a linear fashion, with hystricomorpha species displaying a generally elevated GC- and GC3-content as compared to primates. The squirrel-like sciuromorpha show further an increase in their respective GC content, partially exiting the linear GC/GC3 relationship. Most of the supermyomorpha, however, excluding C. canadensis and N. galili, form a loose, non-linear cloud of data points with exceptionally high GC- and GC3-contents. This “supermyomorpha island” also harbors a tight cluster representing the four investigated lagomorph species.

Many genome releases derive from female source DNA, which yields limited insights into the actual size of the pseudoautosomal region and its boundary to X-specific genes since a direct comparison to the sequence of the Y-chromosome is missing. However, this information is necessary to determine whether there might be a potential diversification of neuroligin-4 into gametologs. Interestingly, very few genome releases are from males, and a respective NLGN4Y gametolog can been identified. This is the case for some primates from the clade catarrhini, including humans, rhesus macaque (Macaca mulatta), and chimpanzee (Pan troglodytes), and the clade platyrrhini with the common marmoset (Callithrix jacchus). Beyond this, NLGN4X/Y gametologs were identified in horses (Equus caballus), gerbils (M. unguiculatus, this paper), and D. spectabilis (this paper) (for a list of the coding sequences, see Supplementary Data Sheet 6). Given the limited number of species’ genomes from male specimen, it appears feasible to assume that in an early stage of mammalian evolution NLGN4X and NLGN4Y developed as gametologs, like AMELX/Y, and that NLGN4Y still needs to be identified in the corresponding male samples. This, in turn, would mean that all sequenced NLGN4Y orthologs should have a single common ancestor. To scrutinize this notion and to increase the resolution of our analysis, we used the underlying protein coding nucleotide sequences and not the amino acid sequences to infer an evolutionary tree for given NLGN4X and NLGN4Y pairs of six different species and the gene variants in D. spectabilis (Figure 5 and Supplementary Data Sheet 6; for branching statistics cf. Supplementary Figure 4). Strikingly, the corresponding tree indicates five independent incidences of neuroligin-4 diversification into gametologs. Within the primate order, the NLGN4X/Y diversification occurred within the clades catarrhini and platyrrhini independently. A comparable scenario occurred in rodents (gerbils and D. spectabilis), and an additional independent separation happened during the evolution of horses in the order perissodactyla.

Given the close evolutionary relationship between gerbils and mice, we hypothesized that some of the PAR genes conserved in mice (i.e., Sts, Asmt or Akap17a, Maxeiner et al., 2020) might still be present in the close vicinity of the neurolgin-4 gene orthologs in gerbils. Indeed, the contigs (VFHZ01026354.1[X-chr.], VFHZ01013942.1[Y-chr.]) with the sequence information of NLGN4X and NLGN4Y contained ASMT, AKAP17A, ASMT, DHRSX, CD99, XG, ARSE, and STS within a stretch of approx. 80 kilobases, respectively (Figure 6). However, these genes appeared to be conserved as potentially functional proteins only on the contig that also contained NLGN4X. On the corresponding contig with NLGN4Y, solely AKAP17A appears to be translatable into a functional protein (“AKAP17AY”). Incomplete sequence information from a sister species of the mangolian gerbil, the fat sand rat Psammomys obesus, indicates that a contig (NESX01053770.1) with the neuroligin-4 gene ends downstream of exon 7 with the repetition of a telomere motif. Thus, gerbils might not retain any PAR-like homology among the remaining ancestral autosomal genes, possibly explaining the observation that synaptonemal complexes do not occur during male meiosis (De la Fuente et al., 2007).