Zrsr2 and functional U12-dependent spliceosome are necessary for follicular development.

ZRSR2 is a splicing factor involved in recognition of 3'-intron splice sites that is frequently mutated in myeloid malignancies and several tumors; however, the role of mutations of Zrsr2 in other tissues has not been analyzed. To explore the biological role of ZRSR2, we generated three Zrsr2 mutant mouse lines. All Zrsr2 mutant lines exhibited blood cell anomalies, and in two lines, oogenesis was blocked at the secondary follicle stage. RNA-seq of Zrsr2 mu secondary follicles showed aberrations in gene expression and showed altered alternative splicing (AS) events involving enrichment of U12-type intron retention (IR), supporting the functional Zrsr2 action in minor spliceosomes. IR events were preferentially associated with centriole replication, protein phosphorylation, and DNA damage checkpoint. Notably, we found alterations in AS events of 50 meiotic genes. These results indicate that ZRSR2 mutations alter splicing mainly in U12-type introns, which may affect peripheral blood cells, and impede oogenesis and female fertility.

Introduction

Alternative splicing (AS) of messenger RNA is an essential eukaryotic mechanism whereby multiple transcripts are generated from a single gene, therefore contributing to protein diversity (Nilsen and Graveley, 2010). The splicing machinery is composed of five small nuclear ribonucleoproteins (snRNPs) that together with different splicing factors catalyze the splicing reaction. There are two different splicing machineries: the major class (named U2-dependent spliceosome), which removes the majority of introns, and the minor class or U12-dependent spliceosome, which removes U12-type introns (<0.4% of all introns) (Burge et al., 1998). U12-type introns are highly conserved across distant taxa (Bartschat and Samuelsson, 2010; Russell et al., 2006) and coexist with U2-type introns in essential genes involved in development, RNA processing, DNA replication, and cell cycle functions (Merico et al., 2015; Patel and Steitz, 2003; Turunen et al., 2013). The minor spliceosome is necessary during development in several plants and animal species (Bai et al., 2019; Gault et al., 2017; Markmiller et al., 2014; Otake et al., 2002; Xu et al., 2016) and also plays a role in the mouse central nervous system (Baumgartner et al., 2015, 2018), hypothalamus (Alen et al., 2019), early embryo development (Gomez-Redondo et al., 2020), and spermatogenesis (Horiuchi et al., 2018). Moreover, alterations of the minor spliceosome have been associated with several diseases (Edery et al., 2011; Elsaid et al., 2017; He et al., 2011; Horiuchi et al., 2018; Merico et al., 2015) and have a broader impact on human diseases than previously appreciated (Olthof et al., 2020).

The splicing factors ZRSR1 and ZRSR2 (encoded by Zrsr1 and Zrsr2 genes, respectively) are important components of both the major and minor spliceosome and have been implicated in 3′ splice site recognition of U12-type introns and in the second step of U2-type intron splicing (Shen et al., 2010). These proteins are highly homologous to the 35-kDa subunit of the splicing factor U2AF, which is implicated in the recognition of the 3′ splice site of major introns. U2AF35-related proteins include U2AF26 (encoded by U2afl4 gene), ZRSR1, and ZRSR2, and they all share the same basic domains: an RNA-recognition motif (RRM) flanked by two zinc-finger domains and one arginine-serine rich region (RS) (Kitagawa et al., 1995; Shepard et al., 2002; Tronchere et al., 1997). Redundant expression of these factors has been observed in many tissues, so in these cases, the lack of expression of one of them could be compensated by its paralog.

Various sequencing studies have identified somatic ZRSR2 mutations in hematological diseases, such as myelodysplastic syndrome (MDS), chronic lymphocytic leukemia (CLL), chronic myelomonocytic leukemia (CMML), or thyroid cancer, describing mis-splicing of U12 (Inoue et al., 2021) or both U12- and U2-type introns (Madan et al., 2015, 2020, 2021). In addition, ZRSR2 is frequently mutated in solid and non-solid tumors (Bejar, 2016). It has been suggested that ZRSR2 is affected by nonsense or frameshift mutations, which presumably results in loss of ZRSR2 (Escobar-Hoyos et al., 2019). However, it has been reported that the most common alterations in ZRSR2 associated with cancer are ZRSR2 missense and synonymous substitutions (∼50%), and only 30% are nonsense substitutions or frameshift alterations that produce ZRSR2 loss (Tate et al., 2019). Analysis of human ZRSR2 mutations revealed that 46% of them are located within the 10 nucleotides adjacent to a 3′ss or 5′ss, suggesting that those mutations could affect the splice site recognition by the spliceosome. Interestingly, a study using Zrsr2 knockout mice has reported that ZRSR2 is not essential for hematopoietic stem cell (HSC) development in mice, as the transgenic mice were phenotypically normal, and the authors suggested that ZRSR1 (which is a retro-transposed copy of ZRSR2) can compensate for the deficiency of ZRSR2 (Madan et al., 2021). Similarly, it has been observed that Zrsr1 knockout mice do not show an abnormal phenotype (Sunahara et al., 2000). However, a conditional Zrsr2 knockout mouse line has revealed certain roles of this splicing factor in hematopoiesis, which enhances HSC self-renewal (Inoue et al., 2021). Moreover, using Zrsr1 mutant mice we recently demonstrated that Zrsr1 has an essential role in muscle strength, hematopoiesis, and spermatogenesis, producing severe alterations at the spermatocyte stage, causing retention of U12-type introns and showing functional resemblance between Zrsr1 and Zrsr2 (Horiuchi et al., 2018). Furthermore, we have described that the double mutation of Zrsr1 and Zrsr2 led to the developmental arrest of embryos at the two-cell stage, indicating that their expression is essential for early embryo development (Gomez-Redondo et al., 2020).

To further determine the role and function of Zrsr2 in vivo, we generated three Zrsr2 mutant mouse lines containing nonsense mutations targeting the region immediately after the AUG initiation and in the two zinc finger domains. We hypothesize that these mutations in homozygosity could interact with the correct performance of minor splicing processing, thus allowing the identification of their specific functions. We found that mice containing mutations within the zinc finger domains produced truncated proteins that exert a recessive negative effect on ZRSR2 function. Zrsr2 homozygous mutant (Zrsr2) mice were viable but showed alterations in peripheral blood cells. Also, Zrsr2 females exhibited impaired oogenesis, with follicles blocked at the secondary stage and eventually leading to infertility, in the same way in which Zrsr1 mutants present spermatogenesis alterations (Horiuchi et al., 2018). RNA sequencing (RNA-seq) analysis of cumulus-oocyte complexes (COCs) of Zrsr2 secondary follicles showed an aberrant splicing pattern, causing retention mainly in U12-type introns. Altogether, this study confirms the function of ZRSR2 for splicing of U12-type introns and reveals its role in oogenesis and peripheral blood cells.

Results

Zrsr2 mutant female mice display severe defects in oogenesis

We used CRISPR-Cas9 technology to produce Zrsr2 mutant mice with different nonsense mutations, generating three different mutant lines named Zrsr2, Zrsr2, and Zrsr2. Mutant Zrsr2 line was generated by editing the first exon of Zrsr2 (Figure 1A, top panel), immediately after the AUG initiation codon (Figures S1A and B). Three different sublines were produced, with deletions of 23, 4, and 22 nucleotides (Zrsr2, Zrsr2, and Zrsr2, respectively) (Figure S1B). In Figure 1B we can see a western blot analysis of three selected transgenic sublines (Zrsr2, Zrsr2, and Zrsr2) using spleen from adult mice. Although all three mutant lines of Zrsr2 produced frameshift mutations and theoretically short, truncated proteins, in the western blot we can see that Zrsr2 produces a protein of 399 aa that could be originated from an alternative start codon located upstream of the first zinc finger domain of Zrsr2 (Figure 1A, second panel, Figure 1B, second lane). Similar relative mRNA expression of Zrsr2 was detected in WT and Zrsr2 (Figure 1C).

Figure 1

Diagrams of Zrsr2 gene structure and mutations generated in this study and their expression

(A) Diagram of the Zrsr2 gene with its functional domains: ZF1 and ZF2 (zinc finger domains 1 and 2), RRM (RNA recognition motif), and RS (arginine-serine rich). The red lightning indicates the target region for each different mutation. The protein sizes and structures produced by the different mutations (Zrsr2, Zrsr2, and Zrsr2) are represented in the lower panels, shaded in pink. The position of the primers used for Figure 1C (2-3F and 7R) and Figure 1C (9F and 11R) are indicated.

(B) Western blot analysis of ZRSR2 protein expression in spleen of WT and the three mutant mice, using ACTB as a loading control. The positions of molecular mass markers are indicated on the left. Blue arrows indicate ZRSR2 mutant protein produced by different alternative splicing events.

(C) Expression level of Zrsr2 mRNA in ovaries from 3-month-old females of the three mutant lines and WT mice determined by RT-qPCR. Biological triplicate data for qPCR are presented as mean ± SEM. Bars with different superscripts differ significantly (p < 0.05).

(D) Non-quantitative RT-PCR analysis of Zrsr2 in WT and Zrsr2, Zrsr2, and Zrsr2 mice. The right panel shows a diagram of the affected exons in Zrsr2 mice. Isoforms were confirmed by sequencing.

The second mutant line (Zrsr2) was produced by targeting the first zinc finger motif (ZF1) (Figure 1A, top panel) of Zrsr2, generating three sublines with 11-, 7-, and 8-nucleotide deletions (Zrsr2, Zrsr2, and Zrsr2, respectively) (Figure S1C). These alterations are located in a region that contains an exonic splicing enhancer (Figure S1C). These mutations appear to produce a truncated protein that is not recognized by the anti-ZRSR2 antibody in the western blot (Figure 1A third panel, Figure 1B third lane), and their Zrsr2 mRNA expression is significantly lower than the WT and the other two mutant lines (Figure 1C).

The Zrsr2 line was generated by targeting the second zinc finger motif (ZF2) (Figure 1A, top panel). The transgenic sublines generated comprise animals with 26-, 10-, and 17-nucleotide deletions in Zrsr2 (Zrsr2, Zrsr2, and Zrsr2, respectively) and produced truncated proteins (Figure S1D). From this line we expected to observe different size proteins: a 410-aa protein with a nonsense sequence after amino acid 307 and a 504-aa protein originated by exon 10 skipping because the deletions originated in the mutant lines cover the 3′ss region of exon 10 and the first nucleotide (G) of the 5′-splice site (Figure 1A, bottom panel; Figure S1D). However, in the western blot we observed three bands that could be produced by alternative splicing due to the deletions produced in the mutants (Figure 1B, blue arrows; Figure S1D). No differences were observed in the levels of mRNA expression between WT and Zrsr2 mice (Figure 1C). This observation was also consistent with the RT-PCR electrophoresis results, which clearly show the Zrsr2 spliced fragments, whereas no alterations in splicing are detected in Zrsr2 or Zrsr2 (Figure 1D).

The mutant lines did not show any phenotype in heterozygosis; for this reason, all analyses were performed in homozygous animals. Because all three sublines in each case showed similar phenotypes, the following phenotypic results refer to Zrsr2, Zrsr2, and Zrsr2, hereafter named Zrsr2, Zrsr2, and Zrsr2, respectively.

Although homozygous transgenic mice developed normally and males were fertile (with normal spermatogenesis and sperm parameters), females from the transgenic lines Zrsr2 and Zrsr2 showed smaller ovaries than wild-type (WT) females (Figures 2A–2C). Females from the Zrsr2muA line exhibited no differences in ovary size compared with WT (Figures 2A–2C). Histological observation of homozygous mutant ovaries from the three lines revealed arrested folliculogenesis at the secondary follicle stage without the appearance of preovulatory follicles in transgenic lines Zrsr2 and Zrsr2 (Figures 2B and 2D), whereas no deficiencies were observed between Zrsr2muA ovaries and WT. Because alternative splicing may contribute to gonadotropin expression (Das and Kumar, 2018), the ZRSR2 mutation could also affect oogenesis by modifying AS of genes that regulate pituitary gonadotrophins affecting follicle development. For this reason, we examined whether exogenous gonadotropins could induce ovarian stimulation using a classic superovulation protocol (Ramos-Ibeas et al., 2014) and whether in vitro maturation could produce follicular development. After pregnant mare’s serum gonadotropin (PMSG) injection, females from the Zrsr2 line showed no differences with respect to the WT, and both Zrsr2 and Zrsr2 showed increased numbers of secondary follicles and lower early antral and antral follicles than WT ovaries (Figures S2 and S2D) (p < 0.01). Sixteen hours after the second injection of luteinizing gonadotropin hCG, very few corpora lutea were detected in Zrsr2 and Zrsr2 compared with WT ovaries (Figure S2), indicating that exogenous gonadotrophins are unable to restore the Zrsr2 and Zrsr2 anovulatory phenotype.

Figure 2

Ovarian morphology and histology and follicle count in WT and Zrsr2 mutant mice

(A) WT and Zrsr2 mutant ovaries from the three lines at 5-, 8-, and 11-week-old mice. Scale bar, 200 μm.

(B) Hematoxylin and eosin-stained sections of WT, Zrsr2, Zrsr2, and Zrsr2 ovaries from 5- and 8-week-old mice. Scale bar, 200 μm.

(C) Bar graph representing the ovary length of WT, Zrsr2, Zrsr2, and Zrsr2 ovaries from 5-, 8- and 11-week-old mice (n = 5 per line). Biological triplicate data for qPCR are presented as mean ± SEM. Bars with different superscripts differ significantly (p < 0.05). (two-way ANOVA followed by Tukey’s post hoc test).

(D) Primary, secondary, and antral follicles count in dissected ovaries from 8- to 10-week-old females from the WT and the three Zrsr2 lines (n = 5 per line) injected with pregnant mare’s serum gonadotropin (PMSG). Biological triplicate data for qPCR are presented as mean ± SEM. Bars with different superscripts differ significantly (p < 0.05) (two-way ANOVA followed by Tukey’s post hoc test).

Moreover, COCs from the Zrsr2 line obtained 48 h after PMSG injection were unable to properly mature in vitro and showed several abnormalities, such as a few or absence of cumulus cells, large perivitelline space, abnormal zona pellucida, non-spherical oocytes, and whitish oocytes (Figure S3). This indicates that the blockage of follicular development is not mainly related to some hypothalamic effect over hormonal stimulation of oocyte development but to a local effect of mutant Zrsr2 expression.

Zrsr2 mutant mice exhibit blood cell anomalies

Hematological evaluation of WT and Zrsr2 mutant mice revealed alterations in peripheral blood cells in all mutant lines (Zrsr2, Zrsr2, and Zrsr2). The number of white blood cells was higher in mutant mice compared with WT (Figure 3), which goes in accordance with a previous study (Inoue et al., 2021). Specifically, more granulocytes and monocytes were present in all mutant lines, more lymphocytes were found in Zrsr2, and no differences were found between mutant lines in the number of monocytes or granulocytes (Table S1). Defects in red blood cell series were also observed, with the mutants showing differences in the erythrocyte volume and the hemoglobin content, as well as in the variation coefficient. Zrsr2 and Zrsr2 also showed an increased number of platelets similarly to the conditional Zrsr2 knockout mice (Inoue et al., 2021), displaying mild thrombocytosis, and all three mutants had a higher percentage of blood occupied by platelets than the WT (Figure 3, Table S1). These overall observations are consistent with the described role of Zrsr2 in hematopoiesis and could be attributed to CLL or CMML, although at the time point of evaluation, the animals appeared healthy, the life expectancy of the mice was normal, and the appearance of lymphadenopathy was not observed.

Figure 3

Peripheral blood cell counts

Boxplot showing the most representative parameters of each cell type series (white blood cells, red blood cells, and platelets). WBC, white blood cells; LYM, lymphocytes; MON, monocytes; GRAN, granulocytes; RBC, red blood cells; MCV, erythrocyte volume; MCH, hemoglobin content; RDW, variation coefficient; PLT, platelets; MPV, mean platelet volume; PDW, heterogeneity; PCT, percentage of blood occupied by platelets. Bars with different superscripts differ significantly (p < 0.05) (one-way ANOVA followed by Tukey’s post hoc test).

Transcriptomic analysis of Zrsr2 secondary follicles

To establish the molecular basis for the oogenesis defects observed in Zrsr2 mutant mice, we first evaluated the Zrsr2 mRNA expression profile in developing follicles to identify where the highest expression of Zrsr2 occurs (Figure 4A). Then, we undertook RNA-seq analyses on WT and Zrsr2muC secondary follicles, which was the cell type that had the highest expression of Zrsr2 and the developmental stage where follicular development was blocked in the mutant ovary. Recent investigations have shown that the expression of minor intron-containing genes (MIGs) across different mouse and human tissues is dynamic (Olthof et al., 2019), but no data from ovary or follicles were analyzed. To characterize the expression of these MIGs in secondary follicles, we followed the pipeline described in Olthof et al. (2019) to evaluate the expression of 666 MIGs in secondary follicles RNA-seq data obtained from WT mice. Then, we merged our data with the rest of the tissues analyzed and set the expression threshold at 1 transcript per million, finding that secondary follicles express a similar number of MIGs (557) than tissues with a higher general expression like testis (585) or thymus (553) (Figure S4A). Moreover, not only many MIGs are expressed in secondary follicles but they are expressed at a relatively high level (Figure S4B). We detected two groups of tissues with similar MIG expression patterns, being the secondary follicles clustered with the bone, lung, thymus, spleen, cerebrum, and testis in the group with higher expression (Figure S4C). This analysis allowed us to identify 22 MIGs that constituted a specific signature of secondary follicles, being significantly up-regulated (fold-change>2) compared with all other tissues (Figure S4D). These genes are involved in functions such as RNA binding (Myef2, Esrp1, Gar1, Lsm5), transcription (Myef2, E25f, Polr2k), and splicing (Esrp1). We did not identify any MIG significantly down-regulated in secondary follicles compared with the rest of the tissues analyzed. The expression of a high number of genes carrying minor introns in the secondary follicles suggests that these genes play an important role in this stage of development.

Figure 4

Zrsr2 expression profile and RNA-seq analysis of WT and Zrsr2 follicles

(A) Relative expression levels of Zrsr2 in developing follicles; biological triplicate results presented as the mean ± SEM. Bars with different superscripts differ significantly (p < 0.05 by two-way ANOVA followed by Tukey’s post hoc test).

(B) Pie charts show the proportion of up-regulated and down-regulated genes as well as the cellular components of the DEGs.

(C) Volcano plot of differentially expressed genes between Zrsr2muC and wild-type samples, in which genes with an false discovery rate value under 0.01 and a fold change over 2 are represented in red. The genes with the greatest fold change and/or significant difference are indicated. Four genes of calcium-binding protein are highlighted surrounded by blue.

(D) Bar plot of the number of genes belonging to the three main gene types in the set of DEGs.

(E) Number of DEGs that are expressed preferentially in the oocyte, in the granulosa cells or both (oocyte expression data: GEO-GSE111687, granulosa cells expression data: GEO-GSE158218). Missing genes in those datasets are indicated as “Non-classified.”

(F) Significantly enriched molecular functions terms in down-regulated transcripts in the Zrsr2muC follicles.

(G) Significantly enriched biological process terms in down-regulated transcripts in the Zrsr2muC follicles.

(H) RNA-seq coverage plot of Zrsr2 gene (top panel) representing the WT data (green) and the Zrsr2 data (red). The lower panel shows a zoom of the last exons of Zrsr2, showing the coverage and the splice junctions (shown in blue) of WT and Zrsr2 data.

We performed gene quantification of our RNA-seq data, keeping only those genes that were detected in at least two samples, analyzing a total of 31,952. In Zrsr2 follicles, 227 and 198 genes were up- and down-regulated, respectively, when compared with WT (Figure 4B, Table S2A and S3). We detected expression differences in 10 genes that contain U12-type introns. The volcano plot of DEGs shows similar distribution between up- and down-regulated genes (Figure 4C). Among the genes with the greatest fold change and/or significant difference there are several pseudogenes and four genes of calcium-binding protein (C1rb, C1s1, Duox2, and Efemp1). Among the DEGs there are mainly protein-coding genes and to a lesser extent lncRNAs and pseudogenes (Figure 4D). The majority of the DEGs detected correspond to genes that are preferentially expressed in the granulosa cells (Figure 4E). However, we observed preferential down-regulation of those genes that are expressed in the oocyte, detecting 67 down-regulated DEGs. In contrast, most up-regulated genes are expressed exclusively in the granulosa cells (193 DEGs). These observations indicate that the blockage in follicle development produced by the Zrsr2 mutation affects gene expression in both the oocyte and the granulosa cells.

DAVID-based gene ontology (GO) analysis of DEGs showed a marked dichotomy between up-regulated and down-regulated genes. The down-regulated genes in Zrsr2 follicles were mainly enriched in helicase activity, ATP and protein binding, DNA replication, cell cycle, and DNA unwinding GO terms (Figures 4F and 4G and Tables S4A and S4C), pointing to an effect of ZRSR2 in essential steps of cell division. This result is consistent with the reduced size of the mutants' ovaries and with the low number of granulosa cells present around the Zrsr2 oocytes (Figure S3). In contrast, genes up-regulated in Zrsr2 follicles were enriched in terms linked to immune system processes and innate immune responses (Tables S4B and S4D), indicating alteration of the immune mediators of follicular homeostasis (Herath et al., 2007) and suggesting that inflammation could be one of the mechanisms responsible for the follicular block of the mutant Zrsr2 mice. Separated GO analysis of differentially expressed genes was also performed after separating those genes that are specifically expressed in the oocyte and the granulosa cells. No relevant differences were found for those genes down-regulated in Zrsr2, but we observed that the enrichment of terms related to immune processes is exclusively associated with those up-regulated genes expressed in the granulosa cells (Table S4E), whereas meiotic nuclear division term was enriched in the oocytes.

RNA-seq results allowed us to accurately map the splicing event generated in our mutant Zrsr2 line, where 43% of the reads mapped against exon 10 and 57% of the reads indicated skipping of exon 10, (Figure 4H) confirming the results observed in RT-PCR (Figure 1D).

Intron retention is a frequent splicing alteration in Zrsr2muC secondary follicles

Alternative splicing events were categorized using vast-tools software (Tapial et al., 2017) into five different groups: 3′ AS site (3SS), 5′ AS site (5SS), exon skipping (ES), intron retention (IR), and micro-exon skipping/alternative micro-cassette exon ≤15 nucleotides (MIC) based on their inclusion levels. Overall, 2,564 differential splicing events were identified in Zrsr2muC follicles, 1,343 up-regulated and 1,221 down-regulated (Tables S2B and S5). All categories of alternative splicing (AS) events were affected (Figure 5A), and we detected a significant overrepresentation of IR, 3ss, and 5ss events when compared with the overall distribution of the event categories in the mouse transcriptome (Figure 5B). We were able to identify differentially spliced events in very few genes that also showed differential expression (Table S6), suggesting that the altered pattern of splicing does not affect directly the DEGs.

Figure 5

RNA-seq analysis of alternative splicing in Zrsr2muC follicles

(A) Distribution of categories of alternative splicing (AS) events differing in Zrsr2muC secondary follicles versus WT. Percentages of each class of event in vast-tools (VastDB annotation) and in Zrsr2muC follicles are indicated. 3SS, alternative 3′ splice sites; A5SS, alternative 5′ splice sites; ES, exon skipping; MIC, alternative micro cassette exon ≤15 nucleotides; IR, intron retention.

(B) Enrichment of differentially spliced events in Zrsr2muC secondary follicles versus WT with respect to the total number of events annotated in the mouse genome. MIC values are chopped for clarity purposes.

(C) Violin plots of AS events differing in Zrsr2muC secondary follicles versus WT.

(D) Differences in intron retention events detected in Zrsr2mu compared with WT follicles (measured as delta percent-spliced-in [DPSI] ratio of intron read counts in Zrsr2muC versus WT for different intron categories, as indicated).

(E) Distributions of the locations of affected U2-type introns in Zrsr2muC follicles relative to U12-type introns present in the same gene.

(F) Percentage of events altered in Zrsr2muC compared with WT follicles that were also detected in Zrsr1 testis (Horiuchi et al., 2018). DEG, differentially expressed genes; Other AS: 3SS, 5SS, ES, and MIC; IR (U2), U2-type introns that showed IR; IR (U12), U12- and U2-type introns located in genes with U12 that showed IR. The table on the right shows the genes with common IR (U12) events between the two experiments.

(G) Significantly enriched GO terms and KEGG pathways in MIGs showing retention in U12- or U2-type introns in genes containing U12-type introns in Zrsr2muC.

(H) Significantly enriched KEGG pathways in MIGs showing retention in up-regulated U2 introns in Zrsr2muC

The violin plot showed similar up- and down-regulated alternative splicing events, except for IR, which has an increase of up-regulated retentions (Figure 5C). Increased IR events have also been reported in MDS bone marrow samples harboring ZRSR2 mutations (Inoue et al., 2021; Madan et al., 2020), in ZRSR2 mutant cells with myelodysplastic syndrome (Madan et al., 2015), and in Zrsr1mu mice (Horiuchi et al., 2018). As Zrsr2 has been attributed roles in both U2- and U12-type intron splicing (Shen et al., 2010), we distinguished between these two intron classes. In fact, patients with MDS and diverse ZRSR2 mutations showed that one-third of U12-type introns exhibited significantly increased retention without affecting U2-type intron retention (Inoue et al., 2021). In addition, mice with HSC deficiency for Zrsr2 exhibited global increases in U12-type intron retention without significant changes in U2-type splicing (Inoue et al., 2021). Here we found that 22 of 417 (∼5.3%) IR events up-regulated in Zrsr2muC follicles corresponded to the U12 class of introns, which were enriched over the overall proportion of U12-type introns in the mouse (∼0.04%) (Figure 5D, left panel). A general view of the RNA-seq coverage of some representative U12-type introns is shown in Figure S5. Fourteen percent (58 of 417 events) of IR events up-regulated in Zrsr2muC follicles corresponded to U2-type introns within MIGs, following a similar distribution to that detected in U12-type introns (Figure 5D, central panel). Remarkably, ∼40% of these introns were located immediately upstream or downstream of the U12-type introns, suggesting mutual influence between neighboring U12- and U2-type introns (Figure 5E), in agreement with observations in Zrsr1 mutant mice (Horiuchi et al., 2018). Moreover, we also detected a tendency in intron retention of U2-type introns in Zrsr2muC follicles (Figure 5D, right panel, and Table S2D). The results indicate that the Zrsr1 paralog is not able to offset the splicing alterations in Zrsr2muC oocytes as the maternal imprinting alleles are already methylated at this developmental stage and it is not expressed (Joh et al., 2018; Obata, 2011).

Interestingly, the highest difference in intron retention was found in the GC-AG U2-type intron of the long-noncoding RNA (lncRNA) Rpph1 (Ribonuclease P RNA component H1) (52% more retained in Zrsr2) (Figure S6A). This intron of 21 bp is not identified in Ensembl annotations. Analysis of the RNA secondary structure of the different Rpph1 transcripts revealed that the free energy of the transcript generated from the RNA that retains the intron is significantly higher than the one from the canonical transcript (Figure S6B), suggesting that the canonical transcript is much more stable and thus potentially functional, but no differences in gene expression were observed between mutant and WT.

The highest difference in alternative 5′ splice sites was obtained for Ska1 (Spindle and Kinetochore Associated Complex Subunit 1) (Figure S7), a microtubule-binding subcomplex of the outer kinetochore that is essential for proper chromosome segregation, required for timely anaphase onset during mitosis (Hanisch et al., 2006). It is noteworthy that the alteration occurs in the first exon. The gene with the greatest difference in exon skipping is Hyi (Hydroxypyruvate isomerase [putative]). This AS event also takes place in the second exon, and this gene also shows a high IR of the second intron (Figure S7). Hyi encodes a putative hydroxypyruvate isomerase, which likely catalyzes the conversion of hydroxypyruvate to 2-hydroxy-3-oxopropanoate and may be involved in carbohydrate transport and metabolism.

Modifications in gene expression and AS observed in Zrsr2 follicles were also compared with our previous RNA-seq data for Zrsr1 in testis throughout development (Horiuchi et al., 2018). The percentage of common events was low except for IR of U12- and U2-type introns located in genes with U12, where 20% of the events were shared between Zrsr2 follicles and Zrsr1 testis, corresponding to 11 genes (Figure 5F). These genes are principally related to DNA damage response during genome replication, testis- and ovary-specific version of the TFIIS family of transcription factors, and centrosome and microtubule organization. Validation by quantitative reverse transcription PCR (rt-qPCR) of intron retention changes predicted by RNA-seq analyses from Zrsr2muC follicles (related to Figure 5) and from Zrsr2muB line follicles are indicated in Figure S8. Similar results were obtained by rt-qPCR and by RNA-seq analyses.

Through GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of MIGs showing increased U12- or U2-type intron retention in Zrsr2muC follicles (Figures 5G, 5H, and Tables S7A and S7B), we identified molecular pathways potentially involved in the infertile phenotype, like progesterone-mediated oocyte maturation and MAPK signaling pathways. We also observed enrichment in key cellular pathways such as FoxO and NOD-like receptor signaling pathway and regulation of cyclase activity, among others. GO and KEGG analyses of all genes showing increased U2-type IR in Zrsr2 follicles revealed enrichment in key cellular functions such as regulation of DNA replication, DNA repair, centriole replication and assembly, and centrosome duplication and cycle (Tables S7C and S7D). We also noted enrichment in genes related to myeloid leukemia and cancer (Braf, Ikbkb, Mapk14, Mapk8, etc.) (Baudot et al., 2010). Nevertheless, the most relevant affected KEGG pathways are similar to those observed in MIGs (Figures 5G and 5H, Tables S7B and S7D).

Discussion

In this study, we report a physiological role of ZRSR2 in oogenesis and peripheral blood cells by investigating the effect of three different Zrsr2 mutations in mice. Moreover, by systematic transcriptomic analyses, we identified alterations in gene expression and defective splicing of several U12-type introns in Zrsr2 mutants, as well as aberrations in several splicing of U2-type introns, which could explain the phenotypes described here, confirming the importance of Zrsr2 for white blood cells development, and indicating that Zrsr2 is necessary for oogenesis.

The fact that we only observe phenotypic effects of Zrsr2 mutations in certain tissues suggest that there is a compensation in the expression of Zrsr2 and its paralog Zrsr1; thus, the alterations are only evident when one of these genes is expressed preferentially in a certain tissue, in accordance with our previous studies (Gomez-Redondo et al., 2020; Horiuchi et al., 2018). Therefore, Zrsr1mu has a critical effect on spermatogenesis at the stage when the X chromosome is inactivated, and Zrsr2 is therefore not expressed (Horiuchi et al., 2018). In contrast, Zrsr2 mutations have a critical effect on oogenesis at around 20 days post-partum (dpp), when Zrsr1 is maternally imprinted in the ovary (Obata, 2011). In addition, Zrsr1 is expressed preferentially in red blood cells and Zrsr2 in HSC and multipotent/myeloid/granulocyte progenitors (Lara-Astiaso et al., 2014). For this reason, we will expect that mutations in Zrsr1 affect erythrocyte development, whereas Zrsr2 mutations give rise to MDSs, leukopenia, etc. (Horiuchi et al., 2018; Madan et al., 2015). However, in this work we have observed that mutations in Zrsr2 also affect red blood cells, indicating that Zrsr2mu could be acting as a dominant-negative mutation, and affected not by Zrsr2 absence but by Zrsr2mu interaction with the splicing machinery. Interestingly, although the mutant line A does not show any effect on the development of oocytes, it did produce effects in the blood that even presented more lymphocytes than the other two lines, suggesting that the deletion of 143 aa in the NH2-terminal region is important for the development of peripheral blood cells but it does not affect the development of the follicle. It should be remembered that Zrsr2 has been implicated in 3′ splice site recognition of U12-type introns and in the second step of U2-type intron splicing (Shen et al., 2010), so this deletion in the NH2-terminal region could affect these two functions in different ways. In addition, the truncated ZRSR2 proteins presented here appear to have a more severe phenotype than the Zrsr2 null allele reported previously (Madan et al., 2021), suggesting that these truncations could be recessive neomorphic alleles affecting splicing machinery.

According to our RNA-seq results, the phenotype observed in Zrsr2 and the differences in gene expression could be the outcome of not only the dysregulation and aberrant AS of several MIGs but also of several U2 intron-containing genes. This AS pattern seems to be severely affecting gene expression, causing a preferential down-regulation of those genes that are specifically expressed in the oocyte, thus potentially affecting folliculogenesis. Moreover, up-regulation of genes that are related to the immune system and inflammatory processes in the granulosa cells could indicate the presence of a defense mechanism to prevent the oocytes with an aberrant expression pattern from further developing. Also, the highest up-regulation of the IR events is observed in the GC-AG U2-type intron of Rpph1, a lncRNA that is part of the ribonucleoprotein RNAse P, involved in the cellular translation system, in rRNA and mRNA cleavage, as well as mRNA stability, senescence, inflammation, and chromatin remodeling (Jarrous, 2017; Lee et al., 2021). Up-regulation of Rpph1 has been related to breast and colorectal cancer (Liang et al., 2019); furthermore, it has been suggested that RPPH1 could enhance human acute myeloid leukemia cell proliferation, migration, and invasion (Lei et al., 2019). In zebrafish, there is one maternal-specific variant of Rpph1 present just in oogenesis and embryogenesis (Locati et al., 2018), and in Drosophila, mutations in RNase P induce complete sterility by triggering several DNA damage checkpoints in the oogenesis (Molla-Herman et al., 2015). LncRNAs show overall splicing inefficiency compared with protein-coding genes (Derrien et al., 2012); however, efficient splicing has been observed among lncRNAs with specific functions (Mele et al., 2017), and it has been suggested that GC-AC introns represent new regulatory elements mainly associated with lncRNAs (Abou Alezz et al., 2020). We have not identified in our results a preferential role of Zrsr2 in the splicing of GC-AG introns, or in specific splicing of some U12-type intron, and one question remains unresolved: of all the genes with U12-type introns that are expressed in the follicle, why only some of these introns are affected by the Zrsr2 mutants? One possible explanation is that, because Zrsr1 is maternally imprinted at late-growing stages of oocyte development (Obata, 2011) (after 20 dpp), only genes with U12-type intron expressed after this stage will be affected by the Zrsr2 mutation in transgenic mice.

A total of 50 genes that showed altered splicing of U2-type introns in Zrsr2 are involved in meiotic processes (Table S5, genes with red font), highlighting the importance of a correct splicing regulation in such mechanism. Also, 112 of the differentially spliced events take place in transcription factors or regulator genes (Table S5, genes with green font) (Lambert et al., 2018), indicating that the splicing defects could be producing a secondary effect on gene expression and thus explain at least partially the differences in gene expression taking place between the WT and Zrsr2 follicles. Remarkably, we identified enrichment of IR in 12 genes related to centriole replication and assembly and 15 genes related to centrosome duplication and cycle, essential for formation and relocation of the GV upon meiosis reinitiation (Miyazaki et al., 2005). We also detected abnormal intron retention of Mcm9, described as essential for gametogenesis for its implication in homologous recombination (Lutzmann et al., 2012).

Furthermore, we detected anomalous IR in MIGs that participate in both the FoxO signaling and in the progesterone-mediated oocyte maturation pathways (such as Mapk8, Mapk14, and Braf) (Figure 5G), further supporting the idea of the aberrant splicing pattern being the cause of the mutant phenotype. Recently, it has been observed in complete maturation-deficient human oocytes that alterations in alternative splicing events (primarily associated with metabolism and cell cycle) constitute a critical underlying mechanism governing oocyte maturation (Li et al., 2020).

Conclusion

The significance of minor splicing during gametogenesis is still poorly understood. This comprehensive analysis of Zrsr2 mutant mice reveals a physiological role of Zrsr2 in oogenesis, similar to the established role of Zrsr1 in spermatogenesis (Horiuchi et al., 2018), and describes its role in blood cells. Through transcriptomic analyses, we identified abnormalities in gene expression, and abnormalities in AS, principally in splicing of U12-type introns and U2-type introns present in MIGs, and we identified that the principal genes affected are related to cell cycle, DNA damage checkpoint and repair, oocyte meiosis, FoxO and MAPK signaling pathways, centriole replication/assembly, and centrosome duplication/cycle.

Limitations of the study

Our AS analysis was performed with a bulk of cells, not single cells (oocyte, different types of granulose cells). Our analysis implied that different types of AS might occur in a single cell. Future work at the single-cell resolution will reveal how different types of AS occur in transcripts from the same gene.

Our analysis has not determined the effect on bone marrow. Future work may make it possible to analyze the effects of the different mutations in the bone marrow.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Alfonso Gutierrez-Adan (alfonso.gutierrez@csic.es).

Materials availability

Mouse lines generated in this study are available from the lead contact upon request with a completed Materials Transfer Agreement.

Experimental model and subject details

Generation of Zrsr2 mutant mice

To target the mouse Zrsr2 gene, Cas9-10DA mRNA and sgRNAs (Table S8) were produced with mMESSAGE mMACHINE T7 Ultra Kit and GeneArtTM Precision gRNA Synthesis Kit (Thermo Fisher Scientific, USA), respectively, and injected into B6CBAF1 (C57BL/6xCBA) embryos, which were transferred to pseudo-pregnant females. Pups were genotyped by PCR in standard conditions with Zrsr2 primers (Table S8) and screened for mutations using the SURVEYOR® Mutation Detection Kit (Transgenomic, NE, USA). Founders were confirmed by Sanger sequencing. Three transgenic lines were successfully generated. Homozygotes male and female were used in all experiments. WT Zrsr2+/+ mice and both genotypes (heterozygotes and wt) were used as controls. All experiments were performed in accordance with general ethical principles with the recommendations of the EEC Council directive 2010/63/EU about the protection of animals, which are used for scientific purposes. All study protocols were approved by the Ethical Committee on Animal Experimentation of the INIA (Madrid, Spain) (21 September 2015) and was registered on the Direccion General de Agricultura y Ganadería de la Comunidad de Madrid (Spain) (PROEX 261/15, 4 November 2015; and PROEX 115.5-21, 28 April 2021). Throughout the manuscript we use “Zrsr2muA”, “Zrsr2” and “Zrsr2” to refer to results that are similar between the three generated lines, always carrying the mutation in homozygosity.

Method details

Western blotting

Proteins were extracted from the spleen of the three lines of transgenic and WT adult mice. Briefly, 100 μl of RIPA buffer were added to each spleen sample that was digested for 2 h at room temperature and centrifuged at 16000G for 20 minutes at 4ºC. Proteins were collected with the supernatant and stored at -20ºC for future use. Proteins were resolved by SDS-PAGE (10% acrylamide, loading 20 μg of total protein per well) and transferred onto a nitrocellulose membrane for immunoblotting following standard procedures. To check the amount of transferred protein, red ponceau staining was performed.

Immunodetection was performed by incubating 30 minutes in blocking media (3%BSA in PBS-T) and incubating overnight at 4ºC with 1:1000 mouse-antirabbit Zrsr2 polyclonal antibody (PA5-41797, Invitrogen). Monoclonal anti-mouse ß-actin (ACTB) antibody (A3854, Sigma) was used as the loading control. Incubation with the secondary antibody anti-rabbit IgG-HRP (Santa Cruz Biotechnology, US, sc-517576) was conducted for 2 h at room temperature.

Ovary histology and morphometry, and follicle count

For histological analysis, ovaries were extracted from 5 WT and 5 Zrsr2 (from each of the 3 lines) euthanized mice aged 5, 8 and 11 months, fixed sequentially in 4% PFA (neutral buffered, Sigma-Aldrich, MI, USA) and 70% ethanol, and then paraffin-embedded, cut into 7-μm sections, and stained with hematoxylin and eosin (H&E).

To classify mouse follicular morphology, 11-week-old WT (n=5) and Zrsr2 (n=5 per line) mice were sacrificed by cervical dislocation and their ovaries removed and segregated. First, COCs from pre-antral and antral follicles were recovered by puncture, and then ovaries were sliced to obtain the primary and secondary follicles. Primary and secondary follicles and cumulus-oocytes-complexes were quantified based on morphological classification of mouse follicles (Myers et al., 2004). Follicles were classified as primary follicles if they contained an oocyte surrounded by a single layer of cuboidal granulosa cells, and as secondary if they had more than one layer of granulosa cells with no visible antrum.

Complete blood counts

Peripheral blood was collected from the submandibular vein from 6 WT, 8 Zrsr2, 6 Zrsr2 and 5 Zrsr2 mice by cheek puncture and transported into 20% EDTA-coated tubes. Blood was then subjected to automated analysis using a V-Sight hematology analyzer (A. Menarini Diagnostics, Badalona, Spain) previously calibrated to measure mice blood.

Zrsr2 expression analysis in follicles, COCs and granulose cells

Ten 8- to 10-week-old WT (C57BL/6xCBA) female mice were superovulated by intraperitoneal injections of 5 IU PMSG, and the ovaries were collected 46 to 48 h later. The ovaries were cleaned of any connective tissue and placed in handling medium M2 supplemented with 4mg/ml bovine serum albumin fraction V. Antral follicles were punctured with 30-gauge needles, and immature COCs were collected in handling medium. For gene expression analysis, 5 groups of 20 COCs were frozen in liquid nitrogen. Granulose cells (GC) were isolated by gentle stroking of the inner follicular wall with a pair of fine forceps to release sheets of GCs. The remaining follicular tissue (mainly theca) was discarded. Sheets of GCs were collected, washed three times in PBS by centrifugation at 300 g. Five groups containing samples from 10 GC isolations were frozen in LN2 until RNA extraction. After slicing the ovaries with a scalpel to collect the remaining COCs and granulosa cells, preantral follicles (primary and secondary follicles) were mechanically isolated from the ovaries in M2 medium under a stereo-microscope using fine needles and collected with a micropipette. Five groups of 20 primary and secondary follicles were immediately frozen in LN2 until mRNA purification. The mRNA extraction, reverse transcription (RT) and quantitative PCR procedures used have been previously described (Lopez-Cardona et al., 2017). Briefly, mRNA was extracted using the Dynabeads mRNA Direct Extraction Kit (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer’s instructions with minor modifications. Immediately after extraction, the RT reaction was performed using MMLV Reverse Transcriptase First-Strand cDNA Synthesis (Epicentre, Madison, WI, USA) according to the manufacturer’s instructions. Five groups of cDNA were used to determine the presence of Zrsr2 mRNA, and H2afz and Gapdh were used as the housekeeping genes.

Quantitative PCR

Quantitative real time-PCR (qPCR) was conducted following standard protocols as previously described (Horiuchi et al., 2018). RNA was converted to cDNA using a reverse transcription kit (Applied Biosystems, Foster City, CA). To prime the RT reaction and synthesize cDNA, poly(T) 35 primer, random primers, and Moloney murine leukemia virus reverse transcriptase enzyme were used in a total volume of 40 μL. All qPCR reactions were carried out in duplicate in a Rotorgene 6000 Real-Time CyclerTM (Corbett Research, Sydney, Australia) using a PCR mix (GoTaq® qPCR Master Mix, Promega Corporation, Boston, MA) containing the specific primers selected for Zrsr2 (Table S8) and a 2-μL aliquot of each sample. Two technical replicates were run for all genes of interest. Expression levels were normalized against that of the endogenous controls H2AFZ and/or GAPDH as described previously (Horiuchi et al., 2018). Expression levels were quantified by the comparative cycle threshold (CT) method. Fluorescence was acquired in each cycle to determine the threshold cycle or the cycle during the log-linear phase of the reaction at which fluorescence increased above the background for each sample. Within this region of the amplification curve, a difference of one cycle is equivalent to the doubling of the amplified PCR product. According to the comparative CT method, the ΔCT value was determined by subtracting the CT value obtained for the control gene from the CT value for each gene of interest in each sample. To calculate ΔΔCT, the highest sample ΔCT value (i.e., the sample showing the lowest target expression) was used as an arbitrary constant to subtract from all other ΔCT sample values. Fold changes in relative gene expression levels of target genes were determined using the formula 2−ΔΔCT (Pericuesta et al., 2021).

Fertility and superovulation of female Zrsr2mu mice

Nine- to 11-week-old females from the three transgenic Zrsr2mu lines in homozygosis, heterozygosis, and WT mice (n=10 for each genotype) were included in a continuous mating study. Two female mice were housed with one 12- to 14-week-old male of known fertility, and female mice were rotated weekly. Cages were monitored daily, and the number of pups and litters was recorded. For the ovarian stimulation experiment, female mice were superovulated at the ages indicated in the figures. Five animals per genotype were i.p. injected with 5 units of pregnant mare serum gonadotropin (PMSG) (Sigma) for 48 h, followed by 5 units of equine chorionic gonadotropin (hCG) (Sigma) for further 18h.

In vitro maturation (IVM) of COCs

Six 8- to 10-week-old WT (C57BL/6xCBA) or Zrsr2 female mice (n=5 per line) were superovulated by intraperitoneal injections of 5 IU PMSG, and ovaries were collected 46 to 48 h later. COCs were recovered as previously indicated and only COCs with 3 compact cumulus cells were used. For the IVM experiment, COCs were matured for 17 h in TCM-199 supplemented with 10% (v/v) fetal calf serum (FCS) and 10 ng/ml epidermal growth factor at 37°C under an atmosphere of 5% CO2 in air with maximum humidity.

RNA extraction and RNA-seq analysis of follicles

Total RNA was extracted from 4 pools of 200 secondary follicles of Zrsr2muC (n=17) and 3 pools of 200 secondary follicles of WT (n=28) 2-month-old females using the Arcturus Pico Pure RNA Isolation Kit (Molecular Devices, USA). The purified total RNA was stored in nuclease-free water, and then used for first-strand synthesis. RNA concentration was measured using a Qubit® RNA Assay Kit in a Qubit® 2.0 Fluorometer (Life Technologies, CA, USA). RNA-seq libraries were prepared from the follicle pools described above using the TruSeq RNA kit from Illumina® (NEB, USA) according to the manufacturer’s recommendations, and cDNA libraries were used for sequencing with an Illumina® HiSeq2500 sequencer, generating an average of 54 million 125 bp paired-end reads per sample. The quality of the raw sequences was assessed with FastQC v0.11.9 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Then, the sequences were trimmed using Trimmomatic v0.39 (Bolger et al., 2014), removing adapters and bases with a quality below 30 from the start and end of the reads. Full genomic and transcriptomic reference sequences of Mus musculus (version GRCm39/mm39) were downloaded from Ensembl v103 (Aken et al., 2016). In an initial approach, we followed the pipeline described by Olthof et al. (2019), to evaluate the expression of minor intron-containing genes (MIGs) to compare our results with the data shown in that study. For the rest of the RNA-seq analyses, each sample was aligned against the mouse reference genome and transcriptome using STAR v2.7 (Dobin and Gingeras, 2015), with the following options: --quantMode TranscriptomeSAM –outSAMstrandField intronMotif –outSAMtype BAM Unsorted SortedByCoordinate –outFilterMismatchNmax 999 –outFilterMismatchNoverReadLmax 0.04 –outFilterMultimapNmax 1. The alignments obtained with STAR were sorted with samtools software v1.10 (Li et al., 2009). On average, ∼43 million stranded 125-bp paired-end sequencing reads per sample were aligned uniquely.

Differential gene expression analysis

Raw read counts were calculated from the alignment files with HTSeq-count software v0.11.1 (Anders et al., 2015) setting ‘s’ parameter as yes, as appropriate for stranded reads. Differential gene expression analyses were performed independently using the R packages edgeR v3.32.1 (Robinson et al., 2010) and DESeq2 v1.30.1 (Love et al., 2014). To improve the reliability of our results, only genes identified by both softwares with an adjusted p-value < 0.05 and a FC > 1.5 were considered as differentially expressed. Gene ontology enrichment analysis was performed with the David Gene Functional Classification Tool (Huang da et al., 2009), considering ontology terms as enriched when their p-value was above 0.05.

Classification of oocyte- and granulosa cell-gene expression

To create a classification of the patterns of gene expression in oocyte and granulosa cells, we used public datasets including full genomic expression in oocytes (GEO-GSE111687) and granulosa cells (GEO-GSE158218). We then normalized the expression converting their values into a range of 0-100 to be able to compare them. We considered a gene as tissue-specific when their expression was >2 fold with respect to the other tissue. Genes showing no expression in any of the tissues were marked as “Non-classified”, and genes with expression in both tissues but showing no differences were marked as “Both”.

Differential splicing analysis

To identify differentially spliced events in WT and Zrsr2 mice, we used vast-tools to calculate the levels of inclusion of each transcript in the mRNA in each sample. First, each sample was aligned with default parameters against the mouse reference transcriptome. The samples were combined in a single table with the ‘vast-tools combine’ command. The distribution of each AS event was normalized to the overall number of that event in the mouse transcriptome (VastDB annotation). AS events differentially spliced in the two groups were then identified by calculating differences in their average inclusion levels (delta percent-spliced-in–DPSI), removing those events with low coverage, classified by the software as ‘VLOW’ and ‘LOW’. We considered an event as differentially spliced when its DPSIs was higher than 10, and events were classified as follows: exon skipping (ES), alternative 3’splice site (3’ss), alternative 5’ splice site (5’ss), intron retention (IR) and microexon skipping (MIC). The annotation of U12 type introns used was as described in our previous paper (Horiuchi et al., 2018), updating the U12-type intron database [U12db (Alioto, 2007)] as previously described and including those events annotated by (Olthof et al., 2019). Intron retention events corresponding to U12-type-introns were determined using custom scripts to check for U12 enrichment in the mutant mice.

RNA secondary structure prediction was performed in RNAfold WebServer (Gruber et al., 2008), calculating the structure with minimum free energy (MFE).

Validation of intron retention events

Intron retention events were validated by qRT-PCR using specific primers (Table S8). The qPCR reactions were always carried out in duplicate with three cDNA samples obtained from secondary follicles from WT and Zrsr2 mice (n=10 per sample), different from those employed in the RNAseq experiments. The expression was normalized against that of the endogenous control H2afz and/or Gapdh.

Quantification and statistical analysis

All data were compiled from experiments run in triplicate and reported as the mean ± SEM. All statistical tests were performed using R (R Development Core Team, 2008). Significant differences were determined based on one-way ANOVA or two-way ANOVA followed by Tukey’s post hoc test unless otherwise stated. Significance was set at p<0.05.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Mouse anti-rabbit Zrsr2 polyclonal antibody	Invitrogen	Cat# PA5-41797; RRID: AB_2607733
Anti-mouse β-actin monoclonal antibody	Sigma	Cat# A3854; RRID: AB_262011
Anti-rabbit IgG-HRP	Santa Cruz Biotechnology	Cat# sc-517576; RRID:AB_2848194
Biological samples
Zrsr2 WT mouse secondary follicles	This paper	N/A
Zrsr2 muA mouse secondary follicles	This paper	N/A
Zrsr2 muB mouse secondary follicles	This paper	N/A
Zrsr2 muC mouse secondary follicles	This paper	N/A
Zrsr2 WT mouse ovaries	This paper	N/A
Zrsr2 muA mouse ovaries	This paper	N/A
Zrsr2 muB mouse ovaries	This paper	N/A
Zrsr2 muC mouse ovaries	This paper	N/A
Zrsr2 WT mouse cumulus-oocyte complexes	This paper	N/A
Zrsr2 muA mouse cumulus-oocyte complexes	This paper	N/A
Zrsr2 muB mouse cumulus-oocyte complexes	This paper	N/A
Zrsr2 muC mouse cumulus-oocyte complexes	This paper	N/A
Critical commercial assays
mMESSAGE mMACHINE T7 Ultra Kit	Thermo Fisher Scientific	AM1345
GeneArtTM Precision gRNA Synthesis Kit	Thermo Fisher Scientific	A29377
SURVEYOR® Mutation Detection Kit	Transgenomic	706020
Dynabeads mRNA Direct Extraction Kit	Thermo Fisher Scientific	61011
MMLV Reverse Transcriptase First-Strand cDNA Synthesis	Epicentre	MM070150
Reverse transcription kit	Applied Biosystems	4368814
Arcturus Pico Pure RNA Isolation Kit	Molecular Devices	KIT0204
Qubit® RNA Assay Kit	Life Technologies	Q32852
TruSeq RNA kit	Illumina	RS-122-2001
Deposited data
Mouse reference assembly, release 103, version GRCm39/mm39	Ensembl	http://ftp.ensembl.org/pub/release-103/gtf/mus_musculus
Mouse reference genome, release 103, version GRCm39/mm39	Ensembl	http://ftp.ensembl.org/pub/release-103/fasta/mus_musculus/dna/
Public RNA-seq data of oocytes	GEO	GSE111687
Public RNA-seq data of granulosa cells	GEO	GSE158218
U12-type intron database	Alioto (2007)	https://genome.crg.es/cgi-bin/u12db/u12db.cgi
VastDB	Tapial et al. (2017)	https://vastdb.crg.eu/wiki/Main_Page
RNA-seq data of follicles from Zrsr2 WT and Zrsr2 mutant females	This paper	https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-7806
Experimental models: Organisms/strains
Mouse: B6CBAF1 (C57BL/6xCBA)	This paper	N/A
Oligonucleotides
Primers for Zrsr2muA and WT genotyping, see Table S10	This paper	N/A
Primers for Zrsr2muB and WT genotyping, see Table S10	This paper	N/A
Primers for Zrsr2muC and WT genotyping, see Table S10	This paper	N/A
Primers for amplification of alternative splicing events observed in Zrsr2muC, see Table S10	This paper	N/A
Primers for Zrsr2 gene expression, see Table S10	This paper	N/A
Primers for quantification of intron retention in Zrsr2mu for RT-qPCR, see Table S10	This paper	N/A
Software and algorithms
FastQC v0.11.9	http://www.bioinformatics.babraham.ac.uk/projects/fastqc/	http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
Trimmomatic v0.39	Bolger et al. (2014)	http://www.usadellab.org/cms/?page=trimmomatic
STAR v2.7	Dobin and Gingeras (2015)	https://github.com/alexdobin/STAR
Samtools v1.10	Li et al. (2009)	https://github.com/samtools/
HTSeq-count v0.11.1	Anders et al. (2015)	https://htseq.readthedocs.io/
edgeR v3.32.1	Robinson et al. (2010)	https://bioconductor.org/packages/release/bioc/html/edgeR.html
DESeq2 v1.30.1	Love et al. (2014)	https://bioconductor.org/packages/release/bioc/html/DESeq2.html
David Gene Functional Classification Tool	Huang da et al., 2009	https://david.ncifcrf.gov/
vast-tools	Tapial et al. (2017)	https://github.com/vastgroup/vast-tools
RNAfold WebServer	Gruber et al. (2008)	http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi