Our aim was to amplify and sequence the region of genomic DNA encompassing the exons 1–9 of the Zp4 gene in several species of the subfamily Murinae (Table 1). Sequences were aligned to the corresponding genomic portions of DNA of Mus (chromosome 13) and Rattus (chromosome 17). Sequences of the mRNA obtained from Mus mattheyi, Mus pahari, and Mastomys coucha ovaries were added to determine the limits of the exons. The complete sequences (covering exons 1–9) were obtained only for Rattus rattus and Mus minutoides. For the other species, the coverage ranged between 70% (Apodemus flavicollis and Mus minutoides) and 13% (Lemniscomys striatus and Malacomys longipes).

These new sequences were aligned with Zp4 sequences of other muroid rodents found in Genbank or ENSEMBL (Table 1). The full length alignment (exons 1–9) comprises 36 sequences and 5,140 bp (genomic DNA) and 1,301 bp (coding portion). Two portions of the sequences were scrutinized: the beginning of the gene (from exons 1–3) and the end of our alignment (exons 8 and 9), in both of which stop codons were found in Mus musculus (Goudet et al., 2008). The first fragment in all our samples was successfully amplified (except Lemniscomys striatus) and the second one in thirteen samples.

The results showed that stop codons are present in the first three exons of Zp4 in eight species of the subgenus Mus: M. caroli, M. cypriacus, M. cookii, M. famulus, M. macedonicus, M. musculus, M. spicilegus, and M. spretus (Figure 1); however, they were also present in exons 8 and 9 in M. musculus and M. spretus.

The phylogenetic tree (Figure 2) confirms the monophyly of this subgenus and indicates that the pseudogenization took place after the divergence of the subgenus Mus and before species diversification. Previous studies reported that the four subgenera of Mus diverged between 6 and 7 MYA (Lecompte et al., 2008; Pagès et al., 2012; Meheretu et al., 2015). Within the subgenus Mus the earliest offshoot is estimated to have appeared at around 5 MYA (Pagès et al., 2012), indicating that the pseudogenization took place between 5 and 7 MYA.

To determine whether Zp4 pseudogenization affects only the subgenus Mus, it was necessary to confirm the expression of the four ZP genes and proteins in other subgenera belonging to the genus Mus, in our case Nannomys and Coelomys (Musser and Carleton, 2005; Pagès et al., 2012). Individuals of two species—Mus mattheyi (subgenus Nannomys) and Mus pahari (subgenus Coelomys)—were studied. Furthermore, a species from another genus of the Murinae subfamily, Mastomys coucha, was also analyzed. The species were selected according to their availability and phylogenetic interest for this study.

Using RT-PCR analysis, full-length cDNAs of Mus mattheyi and Mus pahari Zp1 (Supplementary Figure 1) and Zp4 (Supplementary Figure 2) were obtained from total RNA prepared from ovaries. The sequence analysis indicated that they have a complete coding region. The open reading frames (ORFs) encode polypeptides with a theoretical molecular weight of 68.61 and 59.51 kDa (Mus mattheyi Zp1 and Zp4) and 68.37 and 59.54 kDa (Mus pahari Zp1 and Zp4).

These genetic sequences would translate a predictive protein in both species with a high degree of similarity to ZP1 and ZP4 proteins of other mammals (Supplementary Figures 3, 4). In the N-terminal region, a signal peptide is present, whose peptidase cleavage site was predicted by means of the Bendtsen et al. (2004) algorithm. The C-terminal region corresponds to the transmembrane domain (TMD), and is followed by a cytoplasmic tail (Krogh et al., 2001). Moreover, a basic amino acid domain upstream of the TMD may serve as a consensus furin cleavage site (CFCS) (Arg-Arg-Arg-Arg/RRRR). The molecules have a conserved ZP domain, which is present in most sequences of envelope glycoproteins in many species. Upstream of the ZP domain, there is a trefoil domain, characteristic of ZP1 and ZP4 proteins, with six Cys residues, as reported for ZP proteins (Bork, 1993), and between the signal peptide and the trefoil domain a single ZP-N domain at their N-termini, as reported previously (Callebaut et al., 2007; Nishimura et al., 2019) (Supplementary Figures 3, 4). The presence of all these domains indicates that the Zp1 and Zp4 of these rodents could share an apparent similar molecular structure with other ZP proteins.

Thus, taking into consideration that ZP2 and ZP3 are present in all vertebrates described to date and that they have never suffered pseudogenization, the presence of a complete ORF of Zp1 and Zp4 mRNA in these murine species suggests a ZP consisting of four glycoproteins. Nevertheless, partial amplification of Zp2 and Zp3 were made in both species to demonstrate the presence of the four transcripts (Figure 3). Furthermore, Zp1, Zp2, Zp3, and Zp4 mRNAs from Mastomys coucha were also partially amplified (Figure 3).

The mRNA of Zp1, Zp2, Zp3, and Zp4 are effectively translated in proteins, as confirmed by the detection of several peptides belonging to Zp1 and Zp4 in Mus mattheyi and Mus pahari and the four proteins (Zp1, Zp2, Zp3, and Zp4) in Mastomys coucha (Table 2) by MS/MS analyses.

A total of 12, 6, and 13 different peptides corresponding to Zp4 were identified in the different analyses, yielding a sequence coverage of 30.99, 16.74, and 49.08% for Mus mattheyi, Mus pahari, and Mastomys coucha, respectively (Figure 4).

Taken together, these data indicate that four ZP proteins are expressed in Mus mattheyi, Mus pahari, and Mastomys coucha ovaries.

The next question was whether the ZP composition of the egg could interfere with fertilization, for this reason we performed in vitro fertilization experiments with species differing in their ZP composition. Four rodents with different ZP composition were used: Mus musculus with three ZP proteins (Zp1, Zp2, and Zp3) and Mus mattheyi, Mus pahari, and Mastomys coucha with four ZP proteins (Zp1, Zp2, Zp3, and Zp4). The in vitro fertilization rates in a non-competitive context were analyzed (Table 3).

Oocytes from the four species were co-incubated with spermatozoa in conspecific or heterospecific reciprocal crosses. Fertilization success in the conspecific crosses was high only in Mus musculus (79.16%), whereas in the other species the rate was zero or very low (0% in Mus mattheyi, 0.81% in Mastomys coucha, and 3.5% in Mus pahari). When Mus musculus spermatozoa participated in the fertilization the rates were 67.5% in co-incubation with ova from Mus pahari, 11.7% with ova from Mus mattheyi, and 0% with ova from Mastomys coucha. The fertilization rate was very low when Mus mattheyi or Mus pahari spermatozoa were used, except for Mastomys coucha spermatozoa with Mus musculus oocytes (67.14%). Our observations showed that Mus mattheyi and Mus pahari sperm were still able to adhere to the heterologous ZP, but we observed that the number of spermatozoa that adhere to the ZP was much lower than in the control group. Besides, we observed that when the sperm adhered to the ZP, they remained immobile (data not shown). Taking into consideration the results obtained in the in vitro fertilization crosses between Mus musculus vs. Mus pahari (67.5%) and Mastomys coucha vs. Mus musculus (67.14%), it can be concluded that the presence of the 3 or 4 ZP proteins is not a limiting factor for in vitro fertilization.

Genes with a role in fertilization show a common pattern of rapid evolution, which can be attributed to positive selection. An analysis of such selection in Zp4 gene was made in order to know the level of interspecific divergence and the existence of positive selection sites.

The selection analysis pointed to significant positive selection for both muroid datasets. For the dataset of the 28 muroids, both tests (M1a vs. M2a and M7 vs. M8) were significant (p < 0.05 in the former case and p < 0.01 in the latter), and one site (141 R) showed a posterior probability > 95% (Table 4). For the dataset of the 13 muroids, both tests (M1a vs. M2a, M7 vs. M8) were significant (p < 0.01), and one site (547 L) showed a posterior probability > 95% (Table 4).

The DNA and mRNA analyses of Mus mattheyi and Mus pahari Zp4 sequences indicated the presence of a coding sequence for a full-length protein. The alignment showed a high degree of conservation: 87.32 and 88.05% with Mastomys coucha, 76.47 and 77.35% with hamster, 82.94 and 83.49% with rat, and 63.82 and 64.38% with human, for Mus mattheyi and Mus pahari, respectively. At the ZP domain, the 10 cysteines found were conserved in all the species. The furin cleavage site, described in human (Kiefer and Saling, 2002) and rat (Hoodbhoy et al., 2005) coincided with the potential sites for the rodents analyzed (Supplementary Figure 4).

Six (Asn50, Asn74, Asn122, Asn209, Asn226, and Asn299) and five (Asn50, Asn74, Asn122, Asn209, and Asn299) potential N-glycosylation sites were identified in the mature protein in Mus mattheyi and Mus pahari, respectively. Of which, Asn122 and Asn209 were identified by proteomics in Mus mattheyi, indicating that these sites are not glycosylated or not always occupied (Figure 4). The potential N-glycosylation sites Asn50, Asn74, and Asn226 have been identified in rat Zp4 (Hoodbhoy et al., 2005), so they seem to be conserved in Murinae; however, Asn299 has not been detected in the rest of the species analyzed, implying that there are differences in the glycosylation pattern between these species and the rat. Further studies are necessary to identify the glycosylation sites and the type of oligosaccharide chain present.

In mature Zp4 protein, a total of 76 and 73 potential O-glycosylation sites were found in Mus mattheyi and Mus pahari, respectively. Some of the peptides identified contained some of these O-glycosylation sites, 25 in Mus mattheyi and 11 in Mus pahari, so that they would be free or partially occupied in the mentioned proteins (Figure 4). O-glycosylation data are only available for sow and rat ZP4 (Kudo et al., 1998; Hoodbhoy et al., 2005). An O-glycosylated region (Thr296, Ser298, Ser301, Ser304, and Thr312) has been described in rat (Hoodbhoy et al., 2005), and is conserved in Mus mattheyi and Mus pahari (Ser298, Ser301, and Ser304), in addition to Thr296, which is only present in Mus mattheyi (Supplementary Figure 4).

Taking into account that HPLC/MS analysis can be considered a semi-quantitative technique, the fact that the coverage for Zp1 in Mus mattheyi and Mus pahari was 7.22 and 7.98%, respectively, compared to the coverage for Zp4 of 30.99% in Mus mattheyi and 16.74% in Mus pahari, this suggests that Zp1 is less abundant in the ZP than Zp4. These results coincide with those published by Boja et al. (2003) for house mouse (ZP coverage of 56% for Zp1, 96% for Zp2 and 100% for Zp3), and for rat with a protein coverage of 52% for Zp1 and 70% for Zp4 (Boja et al., 2005; Hoodbhoy et al., 2005). They would also agree with the results published for other species, such as rabbit (Stetson et al., 2012) and cat (Stetson et al., 2015), where the coverage of ZP1 is lower than that of ZP4; however, they do not coincide with the results published for hamster, where the coverage of ZP1 was 12.6% and that of ZP4 11.2% (Izquierdo-Rico et al., 2009). Future studies using quantitative proteomics are needed to clarify the ZP protein ratios in different species.

The DNA analysis in different taxa showed the presence of stop codons in the eight species belonging to subgenus Mus and no stop codons in the other Mus species. The gene expression analysis in Mus pahari (subgenus Coelomys) and Mus mattheyi (subgenus Nannomys) clearly demonstrated the presence of four transcripts (Zp1, Zp2, Zp3, and Zp4) using RT-PCR and four ZP proteins using proteomic analysis. These results agree with previous studies that reported the existence of four proteins in the ZP in other placental species like human (Lefièvre et al., 2004), rat (Hoodbhoy et al., 2005), hamster (Izquierdo-Rico et al., 2009), rabbit (Stetson et al., 2012), cat (Stetson et al., 2015), or ferret (Moros-Nicolás et al., 2018c). These four proteins were also present in some marsupials, even though in this group the evolution of the ZP proteins is more complex as there are several copies of the ZP3 gene (Moros-Nicolás et al., 2018a). The four-protein composition could be considered as the ancestral state in eutherian mammals, and ZP1 or ZP4 being lost in some lineages during their evolution.

Until now the ZP composition in the Murinae subfamily was only known in the rat (with four glycoproteins) and house mouse (with three glycoproteins), suggesting that Zp4 was lost after their divergence around 11.2 MYA (Aghová et al., 2018). Using comprehensive taxonomic sampling within the subfamily, especially within the genus Mus (with representative taxa of the four subgenera), we were able to narrowed down the approximate date of the loss of Zp4. First, all other murine genera included in our study seem to have a functional Zp4, suggesting a more recent loss meaning that it occurred within the genus Mus, which diverged around 7.2 MYA (Pagès et al., 2012). Second, we found that all Mus species from the Coelomys, Nannomys and Pyromys subgenera have four glycoproteins, while all species from the Mus subgenus have only three, suggesting that Zp4 pseudogenized early in the Mus subgenus lineage. Previous studies have reported that the four subgenera of Mus diverged between 6 and 7 MYA (Chevret et al., 2005; Lecompte et al., 2008; Pagès et al., 2012), initially the subgenus Coelomys, then Nannomys, and finally Mus and Pyromys (Veyrunes et al., 2006). Within the subgenus Mus, the earliest offshoot is estimated to have appeared at around 5 MYA (Pagès et al., 2012), indicating that Zp4 pseudogenization took place 5–7 MYA.

Pseudogenization events are not rare in the ZP family. ZP1 has been identified as a pseudogene in several species (Goudet et al., 2008; Tian et al., 2009; Stetson et al., 2012; Moros-Nicolás et al., 2018c). The ZPAX gene was lost in mammals before the divergence between marsupials and placentals (Tian et al., 2009). These multiple and independent loss events may be explained in part as the consequence of a duplication event, since ZP1 and ZP4 are paralogue genes (Bausek et al., 2000). After duplication, three evolution events may occurred to the duplicate copies: (a) the ancestral function is partitioned and shared by both copies (subfunctionalization); (b) one gene acquires a new function and the other retains the original one (neofunctionalization), or (c) one gene conserves the original function and the other degenerates to a pseudogene (Cañestro et al., 2013).

The duplication took place in vertebrates before the mammalian divergence (Goudet et al., 2008), and after the duplication event, ZP4 or ZP1 may have pseudogenized in some mammalian lineages, while the remaining ZP protein continued to perform the function of the ancestral gene. In the house mouse, the pseudogenization of Zp4 indicates that Zp1 retained the function of the ancestral gene.

The evolution of male and female reproductive proteins was probably promoted by positive Darwinian selection. Moreover, comparative sequencing studies among taxonomic groups have led to the discovery that reproductive proteins evolve more rapidly than other genes expressed in other tissues (Swanson et al., 2001; Torgerson et al., 2002). This positive selection has been described in seminal plasma proteins (Kingan et al., 2003; Dorus et al., 2004), oviductal proteins (Moros-Nicolás et al., 2018b), and also in other proteins related with fertilization, such as ZP3, CatSper1 or CD9 (Swanson et al., 2001, 2003) These proteins would have been under a selective pressure that may be related to male-female interaction, in this case, sperm-egg interaction. Our analysis identified two sites in Zp4 that are under positive selection (141R and 547L). These results strongly suggest that Zp4 gene has been subjected to positive selection during evolution. These amino acids appear to be important for the protein function. Future studies using direct mutagenesis will be useful to unravel the specific function of these amino acids in ZP4 protein.

Several studies have analyzed the function of different ZP proteins. Recent studies in human demonstrate that mutations in ZP1 gene are related to infertility (Huang et al., 2014; Sun et al., 2019; Yuan et al., 2019), these mutations could affect the shuttling of glycoproteins to the secretory pathway, which would prevent the formation of the ZP around the ova, but also the formation and development of eggs (Huang et al., 2014; Sun et al., 2019; Yuan et al., 2019). The house mouse has provided interesting information on the functions of the different ZP proteins thanks to the use of animals genetically modified as KO and transgenic (Liu et al., 1996; Rankin et al., 1996, 1999, 2001; Dean, 2004). It was demonstrated that Zp1 offers stability and structural integrity to the matrix. KO mice for Zp1 have an abnormal ZP, being more porous; however, these mice are fertile, although their litter sizes are low (Rankin et al., 1999). On the other hand, KO mice for Zp2 or Zp3 present oocytes that are not surrounded by a ZP (Liu et al., 1996; Rankin et al., 1996, 2001). In the case of rodents with four ZP proteins, the roles played by Zp1 and Zp4 remain unresolved.

The study of ZP4 was not possible until now in animal models because, while the KO technology was well-developed in the house mouse (Mus musculus), which has a pseudogenized Zp4, the technique was very difficult to perform in rat, hamster and rabbit. However, development of CRISPR-cas9 technology made it possible to develop the KO technique in species that possess ZP4 (Fan et al., 2014; Bae et al., 2020). In fact, we have recently reported the phenotype of the female rabbit without the ZP4 gene (Lamas-Toranzo et al., 2019). The female rabbit is subfertile and it was observed that this protein is crucial for the embryo development but not for fertilization (Lamas-Toranzo et al., 2019). Moreover, the ZP was significantly thinner, more permeable, and exhibited a more disorganized and fenestrated structure suggesting a structural role (Lamas-Toranzo et al., 2019). The development of KO animals for ZP4 in other species with four ZP proteins, like the rat or the mice presented in this work could also be a useful tool to study the function of this gene. Furthermore, transgenic mice showing a humanized ZP4 have provided valuable information (Yauger et al., 2011). Indeed, these transgenic mice are fertile; however, their ZP is not recognized by human sperm, which means that ZP4 is not sufficient to support human sperm binding to the ZP (Yauger et al., 2011).

A previous study has shown that heterologous fertilization between different species of rodents is possible, although the success is directly related to the phylogenetic proximity of the species (Roldan et al., 1985): the heterologous fertilization rate in vivo and in vitro is considerably lower than the homologous fertilization rate (Roldan et al., 1985; Roldan and Yanagimachi, 1989; Dean and Nachman, 2009; Martín-Coello et al., 2009). Furthermore, in those cases in which embryo culture was carried out, cleavage arrest or embryo degeneration was observed (Roldan et al., 1985). In our study, heterologous IVF was possible when the spermatozoa from Mus musculus had to fertilize the oocytes from Mus mattheyi and Mus pahari, demonstrating that the presence of Zp4 is not involved in the species-specific binding of sperm. This also means that a ZP formed of four glycoproteins is neither a physical nor biological barrier for the spermatozoa of species with a ZP formed of three glycoproteins, at least in in vitro conditions, indicating that Zp4 does not produce any steric hindrance that impedes the specific gamete interaction.

All animal procedures were performed following the Spanish Animal Protection Regulation, RD 1201/2005 which conforms to European Union Regulation 2003/65. All animal experiments were approved by the institutional review board of the University of Murcia according to the guide for Care and Use of Laboratory Animals as adopted by the Society for the Study of Reproduction.

Adult females and males of four species of murine rodents were used: the house mouse (Mus musculus), Matthey's mouse [Mus mattheyi (subgenus Nannomys)], Gairdner's shrewmouse [Mus pahari (subgenus Coelomys)], and the southern multimammate mouse Gairdner's shrewmouse (Mastomys coucha). Mus musculus specimens were of the hybrid strain C57CBAF1, purchased from Harlan Ibérica, (Barcelona, Spain), while Mus mattheyi and Mus pahari were obtained from the “Institut des Sciences de l'Evolution de Montpellier” (Montpellier, France) and Mastomys coucha specimens were obtained from Hobbyzoo (Madrid, Spain), whose species were verified by PCR amplification of the cytochrome b. Animals were kept under standard laboratory mouse conditions in an environmentally controlled room with a 14 h light:10 h darkness photoperiod under constant temperature and relative humidity conditions. Animals were provided with food (Harlan Ibérica, Barcelona, Spain) and water, both available ad libitum. Animals to be used for the experiments were weaned when they were 4 weeks old. Males were kept in individual cages and were used when they were >12 weeks old; after weaning, females were housed together. Mus musculus females were used when they were 6–8 weeks old, and Mastomys coucha, Mus mattheyi, and Mus pahari females were used when they were 6–12 weeks old.

Ovaries were obtained from three different species of mouse: Mastomys coucha, Mus mattheyi, and Mus pahari. The animals were sacrificed by CO₂ overdose, and the ovaries were obtained and frozen in liquid nitrogen and kept at −80°C (for molecular and proteomic analysis) or washed in PBS and used directly (to obtain the ZPs).

Before oocytes were obtained from Mastomys coucha, Mus mattheyi, Mus musculus, and Mus pahari, the animals were subjected to a hormonal treatment to induce superovulation. Females were injected intraperitoneally with 7.5 IU of equine Chorionic Gonadotrophin (eCG) (Sigma-Aldrich, St. Louis, USA), followed 48 h later by 5 IU of human Chorionic Gonadotrophin (hCG) (Lepori Pharma, Spain). The animals were sacrificed 14 h after hCG injection by cervical dislocation and their oviducts were removed. Cumulus-oocyte-complexes (COCs) were obtained from the ampulla of the uterine tube and placed in PBS (for proteomic analysis) or HTF medium (for IVF analysis), COCs were removed or not by gently pipetting into 0.5% hyaluronidase (Sigma-Aldrich, St. Louis, USA).

To obtain the isolated ZPs, animals were subjected to ovarian stimulation and the oocytes were obtained as explained above. Cumulus cells (CCs) were removed by using 0.5% hyaluronidase (Sigma-Aldrich, St. Louis, USA) in PBS, and the ZPs were obtained after vigorous pipetting of each oocyte using a narrow-bore micropipette in PBS, followed by four washes in PBS to eliminate the oocyte debris. ZPs were solubilized for 30 min at 65°C (Accu Block^TM, Labnet, USA), and kept at −20°C until use.

As mentioned above females of the four species were subjected to a hormonal treatment to induce superovulation. In vitro fertilization was performed as previously described by Hourcade et al. (2010). Males were killed by cervical dislocation. The epididymides and vasa deferentia were removed from males of the four mouse species and placed in 1,000 μl of M2 medium, and adipose tissue and blood vessels were removed. The clean structures were placed in a 500 μl drop of Human Tubular Fluid (HTF) medium (BSA supplemented) covered with mineral oil (Sigma-Aldrich, St. Louis, USA), from which spermatozoa were collected. Concentrations were determined with a Bürker hemocytometer. Spermatozoa were incubated in HTF for 30 min at 37°C with 5% CO₂ in air for capacitation.

Fourteen hours after hCG injection, females were sacrificed by cervical dislocation and their oviducts were removed. COCs were obtained from the ampulla of the uterine tube and using a wide-bore pipette tip, placed in 500 μl of HTF medium. Each sample was inseminated with a final concentration of 1 × 10⁶ spermatozoa/ml, and 30 min after, each well was observed under an inverted microscope to assess sperm-oocyte binding. Gametes were co-incubated for 5 h at 37°C under 5% CO₂ in air, after which, remaining CCs and attached sperm were removed by washing in HTF medium with a fine pipette; oocytes were then washed three times in potassium simplex optimized medium (KSOMaa) and placed in culture drop for 24 h at 37°C under 5% CO₂ in air. Nine hours after insemination, the extent of successful fertilization was assessed in a group of presumptive zygotes by pronucleus visualization under a microscope. Another group of presumptive zygotes was confirmed 24 h after, by analyzing 2-cell stage embryos.

Nine hours after insemination, oocytes were incubated in a 100 μl drop of 10 pg Hoechst-33342 dye (bisbenzimide trihydrochloride, Sigma, Madrid, Spain), before being placed on a coverslip for viewing by a fluorescent microscope (Nikon Optiphot-2). The DNA + H-33342 complex was excited with UV at 355 nm light and epifluorescence emission at 465, and photographed. For this, a G 365 excitation filter, an FT 395 dichromatic beam splitter, and an LP 420 barrier filter were used. Both epifluorescent and brightfield photographs were taken using a Coolpix MDC Lens, Nikon, Japan.

Sequences were analyzed to determine the degree of homology with other known sequences using the “BLAST program” (Basic Local Alignment Search Tool) (http://www.ncbi.nlm.nih.gov/blast/). Multiple sequence alignment was carried out using “Clustal Omega” (http://www.ebi.ac.uk/Tools/msa/clustalo/).

The amino acid sequences were analyzed using the software packages “signalP” (www.cbs.dtu.dk/services/SignalP/), “smart genome” (http://smart.embl-heidelberg.de/) to predict the signal peptide and different domains and “NetOGlyc” (www.cbs.dtu.dk/services/NetOGlyc) and “NetNglyc” (www.cbs.dtu.dk/services/NetNGlyc) to predict potential O-linked and N-linked glycosylation sites, respectively. The theoretical protein molecular weight and mature protein molecular weight were calculated with “PeptideMass” from “ExPASy” (http://web.expasy.org/peptide_mass/).

To complete the dataset, ZP4 sequences from different muroid rodents (Table 1) were retrieved from GenBank and Ensembl gene predictions.

All these predictions were checked manually to detect annotation errors especially close to splicing sites. Similarity searches were performed using BLAST and BLAT against assembled genomes in Ensembl followed by a manual compilation of data to predict further genes or exons missing from the Ensembl predictions. It was also checked that the new sequences corresponded to a syntenic region of the corresponding chromosome. Only the exonic portions were kept for the phylogenetic analysis. Translated sequences were aligned using Muscle in Seaview (Gouy et al., 2010). The pseudogene sequences were added afterwards to the nucleotide alignment and manually aligned. The best-fit model of evolution (SYM+G) was determined using the Akaike information criterion (AIC; Akaike 1973), as implemented in jModelTest v2.1.7 (Darriba et al., 2012). Phylogenetic trees were reconstructed using two probabilistic approaches: maximum likelihood (ML) and Bayesian inferences (BI). The ML phylogeny was reconstructed with PhyML (Guindon et al., 2010). The robustness of each node was assessed with 1,000 bootstrap replicates. BI was performed using MrBayes v3.2 (Ronquist et al., 2012). Four independent runs of 10,000,000 generations sampled every 500th generation were performed. A burn-in period was determined graphically using Tracer1.7 (Rambaut et al., 2018). It was also checked that the effective sample sizes (ESSs) were above 200 and that the average SD of split frequencies remained <0.05 after the burn-in threshold. We discarded 10% of the trees and visualized the resulting tree with FigTree v1.4 (Rambaut, 2016). The robustness of the nodes was estimated with Posterior Probabilities (PP).

Selection analyses were made with the muroid datasets, modified to remove the pseudogene sequences, leading to the first alignment of 28 taxa with 711 bp (237 codons) and a second alignment with only the complete Zp4 sequences including 13 taxa (1,644 bp). The analyses were performed with CODEML from PAML4 (Yang, 2007). Data were analyzed under different models: M1a (neutral model), M2a (selection), M7 (beta distribution), and M8 (beta distribution and selection). The likelihood ratio test (LRT) of positive selection was performed on two pairs of models, M1a with M2a, and M7 with M8 (Yang, 2007).

The expression of ZP proteins was studied using proteomic analysis in Mastomys coucha, Mus mattheyi, and Mus pahari ZP. Ovaries were trimmed using small scissors and dissected to remove fat and connective tissue. The solubilized ZP was obtained according to the protocol previously described by our group (Izquierdo-Rico et al., 2009; Jiménez-Movilla et al., 2009). Solubilized ZP was also obtained from oocytes, for which oocyte ZP was solubilized at 65°C in PBS buffer for 30 min. The analysis was carried out on an HPLC-MS system consisting of an Agilent 1100 Series HPLC (Agilent Technologies, Santa Clara, CA) equipped with a μ-well-plate autosampler and a capillary pump, and connected to an Agilent Ion-Trap XCT Plus mass spectrometer (Agilent Technologies, Santa Clara, CA) equipped with an electrospray (ESI) interface.