The X-linked splicing regulator MBNL3 has been co-opted to restrict placental growth in eutherians

In order to identify splicing regulators with a potential placenta-specific function, we began by screening a combination of publicly available and in-house transcriptomic datasets (S1 Table) for splicing factors that are differentially expressed between trophoblastic and embryonic tissues in mouse. The following expression comparisons were made: (i) placenta versus nonplacental tissues; (ii) trophoblast versus inner cell mass (ICM) at the blastocyst stage; and (iii) trophoblast stem (TS) cells versus embryonic and extraembryonic endoderm (XEN) stem cells, for 197 previously defined splicing regulators (Fig 1A–1C, S2 Table; see Methods [25,26]). These analyses revealed the X-linked gene Mbnl3 to be substantially more enriched than any of the other regulators examined in both placenta (analysis (i)) and E3.5 trophoblast (analysis (ii)) (Fig 1A and 1B). Six other genes were also identified as being enriched (log2FC > 1) in trophoblastic tissues in at least one of these comparisons: Mbnl2 (an autosomal Mbnl3 paralog), Esrp1, Esrp2, Rbm20, Rbm38, and Rbms1. For the top candidates, trophoblast expression was confirmed by in situ hybridization analysis (Fig 1D, S1A Fig). This showed Rbm20 to be specifically expressed in the extraembryonic ectoderm adjacent to the epiblast at E6.5 and in the chorionic hinge at late E7.5, regions that have been described as containing TS cells [27–29]. In contrast, expression of Mbnl3 and Esrp2 appeared to be restricted to differentiated trophectoderm (TE) tissue and absent from these regions. Finally, Mbnl2, Rbm38, and Rbms1 were expressed in both TS cell containing regions and more differentiated TE tissues.

Fig 1. Analysis of splicing factor expression enrichment in trophoblastic versus nontrophoblastic tissues.

(A–C) Scatter plots showing trophoblastic expression level and enrichment for 197 splicing regulators in 3 separate comparisons of trophoblastic versus nontrophoblastic tissues. (D) Whole-mount in situ hybridisation analysis of the expression of splicing regulators found to be enriched in trophoblastic tissues (A–C) in mouse embryos at the indicated stages. (E, F) Box plots showing expression levels of the indicated splicing factors in placenta and nonplacental tissues in mouse (E) and in human (Hs), cow (Bt), and opossum (Md) (F). The numerical data underlying this figure can be found in S1 Data. ESC, embryonic stem cell; ICM, inner cell mass; TE, trophectoderm; TPM, transcript per million; TS, trophoblast stem; XEN, extraembryonic endoderm.

https://doi.org/10.1371/journal.pbio.3001615.g001

Importantly, of the splicing factors examined, only Mbnl3 was highly placenta specific, with no or low expression in other tissues (Fig 1E, S1E Fig). Moreover, analysis of RNA sequencing (RNA-seq) data from human and cow showed placental expression of these factors to be highly conserved, with Mbnl3 also showing the strongest placenta-specific enrichment alongside very low levels of expression in nonplacental tissues (Fig 1F, S1C–S1E Fig, S1, S3, and S4 Tables). In contrast, Mbnl3 was strongly expressed across a range of adult tissues in noneutherians (S1E Fig, S1 Table), including the marsupial opossum, for which Mbnl3 showed no placenta-enriched expression (Fig 1F, S1E Fig, S5 Table).

Comparison of the Mbnl3 gene locus across a number of vertebrate species revealed it harbors a novel eutherian-specific transcription start site (TSS), adjacent to a eutherian-specific hAT family DNA transposon MER53, in addition to canonical TSSs, conserved across vertebrates, and more lineage-restricted TSSs (Fig 2A, S2A–S2C Fig). The use of the eutherian-specific TSS, or of rodent or primate specific TSSs found within <5 kbp of it, results in a transcript that codes for an amino-terminally truncated protein (Fig 2B). A similarly truncated MBNL3 protein isoform can also be generated from alternatively spliced transcripts that initiate from the ancestral TSSs but skip exon 2, which encodes the canonical start codon (Fig 2A). Consistent with previous studies, we found that both full-length and truncated MBNL3 protein isoforms were present in mouse placenta by western blot analysis (S2D Fig). To investigate the use of the different TSSs and amino-terminal isoforms, we quantified RNA-seq reads spanning competing exon–exon junctions (a–d in Fig 2A) across a range of tissues from mouse, human, and cow (Fig 2C, S2E and S2F Fig, S6 Table). This analysis revealed the eutherian-specific TSS, alongside the nearby lineage-specific ones, to be dominant in blastocyst stage trophoblast tissue in all 3 species. Furthermore, in placenta, the eutherian-specific TSS was found to account for around half of the transcripts in mouse and almost all of the transcripts in human and cow. Moreover, placenta-specific ChIP-seq data revealed an enrichment for H3K27ac and H3K4me3 active chromatin marks around the eutherian-specific TSS in mouse (S2G Fig). In line with these observations, transgenic reporter assays revealed that this region contains regulatory elements capable of driving trophoblast-specific expression in mouse embryos (Fig 2D).

Fig 2. Molecular evolution of Mbnl3 in eutherians.

(A) Schematic showing alternative splicing events and transcriptional (arrow) and translational (ATG) start sites leading to long and short MBNL3 isoform production. Splice junctions corresponding to transcripts that code for the long (b) and short (a, c, d) isoforms are indicated. Ancestral (black/gray), eutherian-specific (red) and mouse-specific (orange) TSSs and the location of the DNA transposon MER53 (blue) are also shown. (B) Diagram showing long and short Mbnl3 isoforms and the location of eutherian-specific aa changes. Zinc finger domains are shown in blue. Sequence alignment is shown in S3 Fig. (C) Quantification of the relative usage of the splice junctions highlighted in (A) across a range of mouse tissues. (D) Transgenic analysis of a putative Mbnl3 promoter region. A approximately 8-kb region surrounding the eutherian-specific TSS indicated by the black bar in (A) produced lacZ staining in the trophoblast of 4/10 injected embryos at E7.5. (E) Phylogenetic tree showing the degree of MBNL3 (red) and MBNL2 (blue) protein sequence divergence in eutherian (circles) and noneutherian (triangles) species. The tree was produced by the neighbor-joining method using Kimura corrected distances. (F) RNAcompete-derived sequence logos for the indicated human (Hs), mouse (Mm), opossum (Op), and chicken (Gg) MBNL proteins. (G) Analysis of the Z-score contribution of sequences containing GCUU or GCAGC motifs to the total cumulative Z-score of the top 100 RNAcompete-derived 7-mers for the indicated MBNL proteins. (H) Bar plot showing the Z-score of the top RNAcompete-derived 7-mer for each of the indicated MBNL proteins. (I) Analysis of the Z-score contribution of sequences containing GCG, GCC, GCA, or GCU motifs to the total cumulative Z-score of the top 100 RNAcompete-derived 7-mers containing a single GC dinucleotide for the indicated MBNL proteins. (J) Ratio of the z-score contribution of 7-mer’s containing a pair of GC dinucleotides with a 1-bp spacer (GCxGC) versus a 2-bp spacer (GCxxGC). * in panels G–J indicate chimeric proteins. The numerical data underlying this figure can be found in S1 Data. aa, amino acid; MBNL, Muscleblind-like; TSS, transcription start site.

https://doi.org/10.1371/journal.pbio.3001615.g002

MBNL proteins contain 2 pairs of zinc fingers (ZF1-2 and ZF3-4), which mediate RNA binding, separated by a linker region, which is required for splicing regulation and that may play a role in RNA binding (Fig 2B [30–34]). The first pair of zinc fingers is missing from the truncated eutherian MBNL3 isoform (Fig 2B), and equivalent truncated protein isoforms have not been observed for MBNL1 or MBNL2 [35–37]. Using phylogenetic analysis, we found that Mbnl3 has undergone accelerated sequence evolution in eutherians compared to noneutherian Mbnl3 and all Mbnl2 genes (Fig 2E). In particular, eutherian MBNL3 proteins show specific changes in conserved amino acids in ZF1, the spacers separating ZFs 1 and 2 and ZFs 3 and 4, as well as the linker region (Fig 2B, S3A Fig). In contrast, the protein sequence of MBNL2 is highly conserved across vertebrates. Close examination of Mbnl3 gene structure identified an additional feature to be associated with eutherians: loss of a 54-bp alternative exon [38], which triggers nuclear localization upon inclusion [33] and that is found in all Mbnl1 and Mbnl2 paralogs and in most noneutherian Mbnl3 genes (S3B Fig). Altogether, these eutherian-specific changes to Mbnl3, alongside the loss of its expression from nonplacental tissues and its location in the X chromosome, suggest that it may have undergone a process of evolutionary specialization in this clade.

To investigate how these eutherian-specific novelties may have affected MBNL3 functionality, we first examined the RNA-binding preference of full-length (long) eutherian (mouse) and noneutherian (chicken and opossum) MBNL3 proteins using RNAcompete [39–41]. We found that full-length opossum and chicken MBNL3 show a strong GCUU binding preference, similar to the binding preferences that have previously been described for mammalian MBNL1 and MBNL2 and nonvertebrate MBNL proteins using a variety of in vitro and in vivo techniques (Fig 2F and 2G) [38,39,42–44]. In contrast, we found full-length mouse MBNL3 to have overall weaker binding preferences, as measured by the top Z-score values (Fig 2H), and to preferentially bind GCA/U (Fig 2F, 2G and 2I). This motif is similar to that identified by previous CLIP-seq–based studies [38,45,46], confirming the validity of using this in vitro technique to explore MBNL binding preferences. Furthermore, we observed that among 7-mer motifs containing a pair of GC di-nucleotides, opossum and chicken MBNL3 show a strong preference for a pair of uracil’s between the GCs, while mouse MBNL3 appears to favor a single adenosine residue in that position (GCUUGC versus GCAGC; Fig 2F, 2G, 2I and 2J). Altogether, these results demonstrate that the binding specificity of MBNL3 has specifically diverged in eutherians.

Given the high prevalence of the short MBNL3 protein isoform in trophoblastic tissues, we also assessed its binding preferences alongside those of artificially created short isoforms for chicken and opossum MBNL3 and for mouse MBNL1 and MBNL2. The short isoform of mouse MBNL3 was found to have similar GCU/A binding preferences to the long isoform suggesting both may bind the same target sites (Fig 2F–2J). Intriguingly, similar binding preferences were also seen for all the artificially created short MBNL isoforms, with the partial exception of opossum MBNL3, which, although to a lesser extent than its long isoform, maintained a GCU preference (Fig 2F–2J). These results are in agreement with previous work suggesting the binding preferences of canonical MBNL proteins are largely driven by ZF2 [37,47,48] and suggest that the changed binding properties of the full-length mouse MBNL3 isoform may be due to a reduction in ZF2 binding dominance likely driven by the eutherian-specific mutations we identified (S3A Fig).

To investigate the molecular basis of the changes in binding preferences seen for the long isoform of mouse MBNL3, we next generated a number of chimeric mouse–opossum proteins. Introducing the 2 opossum zinc finger pairs into mouse MBNL3 resulted in the mouse protein reverting to an ancestral like binding preference (Fig 2G–2J, S3C and S3D Fig). However, in contrast, introducing the mouse zinc finger pairs into opossum MBNL3 only resulted in a modest change in binding preference. Specifically, although the consensus sequence remained strongly GCU-biased (Fig 2I, S3C Fig), motifs containing a GCA sequence contributed to 20.0% of the accumulated Z-score of the top 100 motifs, while this contribution was only 8.1% for the wild-type (WT) opossum MBNL3 protein (S3D Fig). Given that the linker region may also play a role in RNA binding [30,33,34], we made an additional chimeric protein consisting of the opossum protein with the mouse linker. Similar to the opossum MBNL3 with the mouse zinc fingers, this chimeric protein showed a modest increase in GCA preference at the expense of GCU when compared to the WT protein (21.7% versus 8.1%; Fig 2I, S3D Fig). Together, these data suggest that the novel binding preference of eutherian MBNL3 is a result of sequence changes in both the zinc finger and linker regions.

To explore novel and ancestral physiological roles of Mbnl3 in the developing placenta, we then made use of a mouse line in which both long and short Mbnl3 isoforms had been knocked out [49]. In addition to analyzing Mbnl3 single KOs, we also examined Mbnl2 KOs [37] and the effect of Mbnl2/3 double knockout (DKO), since, despite the divergence of Mbnl3 in eutherians, both proteins share 64% of sequence identity in mouse and show a large degree of expression overlap in differentiated trophoblastic tissues, as assessed by in situ hybridization (Fig 1D) and using a previously published single-cell RNA-seq dataset (S4 Fig [50]). Strikingly, Mbnl3 KO resulted in a significant increase in placental weight at embryonic day 13.5 (E13.5) and E18.5 (Fig 3A, S5A Fig; P < 0.001, Wilcoxon rank-sum tests). Analysis of placentas harvested at earlier stages revealed this weight increase first becomes apparent at E11.5 (S5B Fig). In contrast, no significant changes in placenta weight were seen in DKOs, while placentas from Mbnl2 KOs were found to be slightly smaller than those from their WT littermates at E13.5 (Fig 3A; P = 0.036, Wilcoxon rank-sum test). To further characterize the effects of the KOs on placental development, hematoxylin–eosin staining was performed at E18.5 (Fig 3B), and a number of phenotypic abnormalities of varying severity were observed. In DKOs, these included placental disorganization, with substantial quantities of junctional tissue (spongiotrophoblast and trophoblast glycogen cells) found within the labyrinth zone (Fig 3C) and an increased amount of eosinophilic necrotic tissue within the junctional zone and decidua (Fig 3D). Similar, but much less severe histological abnormalities were observed in Mbnl2 and Mbnl3 single KOs (Fig 3B–3D). Furthermore, an increase in the amount of both labyrinth and junctional zone tissue was seen in Mbnl3 KO placentas, in line with their increased placental weight (S5C–S5E Fig).

Fig 3. Mbnl3 restricts placental growth.

(A, E) Analysis of the effects of M2KO, M3KO, or Mbnl2:Mbnl3 DKO on placental (A) and embryo (E) weight at E13.5 and E18.5. (B) Hematoxylin–eosin staining of WT and KO placenta sections obtained at E18.5. In the zoomed images, note the GlyT and SpT cells within the labyrinth zone of the KO placentas (iv, vi, viii). Also note the nuclei undergoing karyorrhexis (black arrow head) and karyolysis (empty arrow head) in the DKO junctional zone (vii). Scale bar for full placental sections = 1 mm. Scale bar for zoomed images = 100 μm. (C, D) Quantification of the proportion of total junctional zone tissue found within the labyrinth (C) and total area of esophilic necrotic tissue (D) at E18.5. A single medial section was analyzed per placenta. (F) Analysis of the effects on embryo weight of placenta specific M3KO with (M2KOM3^m/+) or without (M3^m/+) universal M2KO. (G) Analysis of the effects on embryo weight of epiblast/fetus specific M2KO with (epiDKO) or without (M2epi) universal M3KO. For (A) and (C–G), the number of placentas/embryos analyzed for each genotype are indicated in brackets, and significance levels are calculated by Wilcoxon rank-sum tests. The numerical data underlying this figure can be found in S1 Data. DKO, double knockout; GlyT, glycogen trophoblast; JZ, junctional zone; LZ, labyrinth zone; M2KO, Mbnl2 knockout; M3KO, Mbnl3 knockout; MBNL, Muscleblind-like; SpT, spongiotrophoblast; WT, wild-type.

https://doi.org/10.1371/journal.pbio.3001615.g003

To assess how these differences in placenta development may affect placenta function and thus embryo development, we next analyzed embryo size. Mbnl3 KO alone had no effect on embryo weight, even at late gestational stages (Fig 3E), indicating that the increase in placenta size that occurs on Mbnl3 KO is accompanied by a reduction in placental efficiency (S5F Fig). In contrast, DKO embryos suffered from intrauterine growth restriction; this was severe at E18.5 (P < 0.001, Wilcoxon rank-sum tests), when placental efficiency was seen to be significantly lower than in Mbnl3 single KOs (S5F Fig; P = 0.015, Wilcoxon rank-sum tests), but could be observed from as early as E13.5 (Fig 3E; P = 0.041, Wilcoxon rank-sum tests). A modest but significant growth restriction was also seen in Mbnl2 single KO embryos at E18.5 (Fig 3E; P = 0.010, Wilcoxon rank-sum tests). To test if the decreased embryo size seen in the DKOs was a direct result of loss of placental and not embryonic Mbnl2/3, we used 2 conditional KO strategies. First, we made use of the X-linked nature of Mbnl3, which means that placental expression occurs almost exclusively from the maternal allele in mouse [8]. Therefore, female embryos in which the maternal allele is KO and the paternal allele is WT are effectively KO for Mbnl3 in the placenta and have mosaic Mbnl3 expression in embryonic tissues (M3^m/+; S5G Fig). Mbnl2 KO embryos with maternal Mbnl3 KO (M2KOM3^m/+) showed a similar growth restriction to full DKO embryos (Fig 3F). Additionally, placentas from offspring with maternal Mbnl3 KO alone showed a similar size increase to those from full Mbnl3 KOs (S5H Fig). Secondly, we made use of a paternal Sox2-Cre mouse line [51] to specifically knock out Mbnl2 from the epiblast (M2epi; S5I Fig). This resulted in a similar fetal and placental size reduction to that seen in full Mbnl2 KOs, suggesting that these phenotypes are at least in part due to a role for this gene in the epiblast and its derivatives and not just in trophoblastic tissues (Fig 3G, S5J Fig). Importantly, unlike with the full Mbnl2 KOs, Mbnl3 KO did not accentuate the fetal growth phenotype of epiblast-specific Mbnl2 KOs (epiDKOs, Fig 3G). Altogether, the data from the conditional KOs show that the strong embryo growth restriction seen in full DKOs is specifically due to a role for Mbnl2 and Mbnl3 in the placenta.

To investigate the molecular bases of these phenotypes as well as to gain insight into the mechanism of action of Mbnl2 and Mbnl3 in placenta, we next performed RNA-seq of placental samples harvested at E13.5 and E18.5. MBNL proteins have been extensively characterized as regulators of alternative splicing; therefore, we began by using this data to assess the effect of each KO on exon inclusion. At both time points, substantially more and larger splicing changes were seen in DKO than single KO placentas (Fig 4A–4F, S6A–S6C Fig), indicating that Mbnl2 and Mbnl3 regulate an overlapping set of events and partially compensate for one another’s loss. In agreement with previous reports [38,46,52], the number and extent of splicing changes in Mbnl2 KOs were greater than those in Mbnl3 KOs (Fig 4A–4C, S6A Fig), and no evidence for substantial Mbnl3-specific splicing effects was observed. We also found no evidence of Mbnl3 acting in placenta as a competitive inhibitor of splicing regulation by Mbnl2 [45,49], as almost all splicing events significantly altered in DKOs were misregulated in the same direction in Mbnl2 and Mbnl3 single KOs (Fig 4C). In line with the known function of MBNL proteins in regulating differentiation-associated alternative splicing [53], we found that misregulated exons acquired more immature-like inclusion levels, as shown both by analyzing splicing changes occurring between E13.5 and E18.5 in WT placentas (Fig 4D) and between embryonic stem cells (ESCs) and differentiated tissues (S6B Fig). Functional enrichment analysis for genes containing misregulated exons identified functions related to cell signaling, migration, and adhesion that could relate to the disrupted placental morphology seen in the DKOs (S7 and S8 Tables). Among the most misregulated events were alternative exons in Numb (Fig 4E), Gpr137 (Fig 4F), and Numa1 (S6C Fig), whose immature-like inclusion patterns have been shown to inhibit differentiation, decrease epithelial tightness, and associate with increased cell proliferation, respectively [54–57].

Fig 4. Effect of Mbnl2, Mbnl3, and double Mbnl2/3 KO on the placental transcriptome.

(A) Bar charts showing the number of differentially spliced exons in M2KO, M3KO, and DKO placentas at the indicated embryonic stages. Exons with increased inclusion in KOs are tallied above the x-axis, while exons with reduced inclusion are tallied below. (B) Venn diagrams showing the overlap in differential spliced exons between M2KO, M3KO, and DKO placentas at E13.5. (C) Boxplots showing the difference in exon inclusion (PSI) in the indicated KOs versus WT placentas for all exons differentially spliced in DKOs versus WTs at E13.5. (D) Boxplots showing the change in PSI in the indicated KO versus WT placentas at E18.5 for all exons with differential inclusion between E13.5 and E18.5 WT placenta. (E, F) Sashimi plots showing usage of Numb exon 9 (E) and Gpr137 exon 3 (F) in placentas with the indicated genotypes at E13.5. The PSI values for the exons in the different genotypes are shown to the right of the plots. (G) Bar charts showing the number of differentially expressed genes in M2KO, M3KO, and DKO placentas at E13.5. Genes with increased expression in KOs are tallied above the x-axis, while genes with reduced expression are tallied below. (H, I) Bar charts showing significance levels for enrichment of TF binding sites among genes up-regulated in M3KO placentas at E11.5. The top 5 most significantly enriched TFs are shown for ChEA (H) and ENCODE (I) data sets. (J, K) Plots showing Myc expression levels in WT and M3KO placentas at E11.5 based on RNA-seq (J) and qPCR (K). In (J), matched placenta pools derived from the same litters are indicated by dashed lines. For (K), values correspond to the fold change between each tested placenta and the average of WT placenta; N = 13 placentas for each genotype. Significance levels were calculated by Wilcoxon rank-sum tests. The numerical data underlying this figure can be found in S1 Data. DKO, double knockout; KO, knockout; MBNL, Muscleblind-like; M2KO, Mbnl2 knockout; M3KO, Mbnl3 knockout; PSI, percentage spliced in; qPCR, quantitative real-time PCR; RNA-seq, RNA sequencing; TF, transcription factor; WT, wild-type.

https://doi.org/10.1371/journal.pbio.3001615.g004

In addition to alternative splicing regulation, MBNL proteins have also been reported to modulate alternative polyadenylation site usage [58], and MBNL3, in particular, has been found to preferentially bind the 3′ end of transcripts [46]. Therefore, we next generated 3′ sequencing data from the same set of placenta samples. As with alternative splicing, we found MBNL2 and MBNL3 to redundantly regulate polyadenylation site usage, and their loss to lead to more immature-like usage patterns (S7A–S7D Fig). However, in contrast to alternative splicing, relatively few genes showed significant alterations in polyadenylation site usage at E13.5 (S7A Fig). Among the misregulated targets (S9 Table) was the transcriptional repressor Mbd2, which showed an increase in the use of a proximal polyadenylation site from approximately 0% of transcripts in the WTs to approximately 50% in the DKOs (S7E Fig). Use of this site results in the production of a shortened protein isoform that promotes pluripotency in ESCs, whereas the long isoform enhances cell differentiation [59,60]. This result suggests that a similar mechanism likely operates in TE tissue.

To further characterize the molecular changes occurring in single and DKO placentas, we next used the RNA-seq data to examine differences in steady-state mRNA levels (hereafter gene expression). Overall, the results from this analysis mirrored those from splicing and polyadenylation site usage, with substantial changes in gene expression being seen in DKOs and only moderate changes in Mbnl2 and Mbnl3 single KOs (Fig 4G, S8A–S8C Fig). Maturation defects were also observed at the level of gene expression in E18.5 KO placentas (S8D Fig). Indeed, among the genes up-regulated in DKO placentas at this stage were a number of transcription factors (TFs) associated with TS cell maintenance, including Esrrb [61], Eomes [62], and Cdx2 [63] (S10 Table). This was accompanied by a reduction in expression of a number of genes involved in labyrinth formation and maturation, including the transcriptional regulators Gcm1 [64] and Cebpa [65], and the Gcm1 target SynB, which plays a role in labyrinth cell fusion [66,67] (S10 Table). Gene Ontology (GO) analyses of differentially expressed genes revealed multiple functions altered in DKOs, including terms related to cell adhesion, migration, and extracellular matrix organization (S11 Table).

In contrast to alternative splicing and polyadenylation usage, analysis of gene expression at E13.5 revealed a subset of genes that appeared to robustly respond to Mbnl3 but not Mbnl2 loss (S8B and S8E Fig). These genes were not significantly misregulated in Mbnl2 KOs and had similar levels of misexpression in Mbnl3 KOs and DKOs (S8E Fig, S10 Table). This indicates that eutherian Mbnl3 has some regulatory activities not possessed by Mbnl2, in line with its unique morphological phenotype. Among these genes were a number with a potential role in nutrient uptake, Slc2a3, LipG, Atp1b1, Steap3, and Slc19a2, which where all up-regulated following Mbnl3 loss (S10 Table). This would be consistent with the X-linked nature of Mbnl3 and the expected maternal bias in its regulatory activities resulting from this; however, it could also be a result of the placenta compensating for functional deficiencies resulting from Mbnl3 loss.

To further investigate the molecular bases behind the increase in placenta growth specifically seen in Mbnl3 KOs, we performed RNA-seq of Mbnl3 KO placentas at E11.5. This is the stage when the increase in placental size was first observed (S5B Fig) and a time of proliferative growth for the organ [68,69]. Differential gene expression analysis at this stage revealed 10 genes to be significantly up-regulated and 23 genes to be significantly down-regulated, over half of which came from a single genomic cluster and were members of the carcinoembryonic antigen (CEA) gene family (S10 Table). Extensive examination of the genes showing significantly altered expression did not reveal any obvious candidates that may be responsible for the observed increase in placenta size. Therefore, we next performed a series of functional and regulatory enrichment analyses using looser cutoffs for differential gene expression (P < 0.1, log2FC > 0.1). These revealed enrichment for GO terms consistent with increased growth among up-regulated genes, including ribosome, translation, and metabolic processes (S12 Table). In addition, they showed that the pro-proliferative TF Myc has the most significant level of binding site enrichment in the up-regulated genes according to both ChEA and ENCODE ChIP datasets (Fig 4H and 4I, S12 Table). Similarly, MycN was found to be the top candidate in a TF perturbation dataset (S12 Table). In line with these results, expression of Myc was up-regulated approximately 1.2 fold at E11.5, both based on RNA-seq data and qPCR assays using independent placenta samples (Fig 4J and 4K). Up-regulation of Myc and its associated regulatory network is consistent with an increase in proliferation contributing to the larger placentas seen in Mbnl3 KO mice.

The increase in placental growth seen on Mbnl3 KO is consistent with the theory that the eutherian X chromosome is enriched for placentally expressed genes that favor increased maternal fitness [9–11]. In this context, Mbnl3 would act by limiting the amount of resources allocated to placenta growth, and thus fetal development. However, while there is an obvious maternal cost in terms of energy expenditure for increasing placenta size in Mbnl3 KO animals, we did not observe any fetal advantage, such as increased embryo size. This could be either due to a suboptimal placenta function in Mbnl3 KOs or because embryo size is already optimal under standard laboratory conditions, where mothers have ad libitum access to food; indeed, in normally gown near-term litters from WT mice, there is no correlation between placenta and fetus size [70]. The second possibility led us to hypothesize that the increased placental size seen in Mbnl3 KOs may provide additional placental reserve capacity to Mbnl3 KO embryos when resources are limited. To test this, we examined the effects of Mbnl3 KO on fetal growth following midgestational maternal calorie restriction, which results in a reduction in embryo size of approximately 50% at E18.5 (S9A Fig [71]). Under these conditions, placentas from Mbnl3 KOs were still found to be larger than those of their WT littermates (Fig 5A and 5B, S9B Fig). However, in contrast to noncaloric restricted conditions and in agreement with our hypothesis, analysis of embryo weights revealed Mbnl3 KO fetuses from calorie restricted pregnancies to be significantly larger than their WT littermates at E18.5 (Fig 5C and 5D, S9C Fig).

Fig 5. M3KO reduces the effect of midgestational maternal calorie restriction on fetal growth.

(A, C) Dot plots showing the weights of individual WT (gray) and M3KO (red, M3KO) placentas (A) and embryos (C) harvested from 6 calorie restricted mothers and the corresponding empirical cumulative distribution plot for each genotype. Significance levels were calculated by a permutation test with 1,000 iterations swapping the genotype labels within each litter. (B, D) Analysis of the effects of maternal calorie restriction on the weight of WT and M3KO placentas (B) and embryos (D) at E18.5. For each embryo/placenta, the weight was normalized with respect to the average of WT and KO weights within each litter. Significance levels are calculated by Wilcoxon rank-sum tests. The numerical data underlying this figure can be found in S1 Data. KO, knockout; M3KO, Mbnl3 knockout; WT, wild-type.

https://doi.org/10.1371/journal.pbio.3001615.g005

We used a list of previously defined RNA-binding proteins and other genes known to be involved in splicing and/or alternative splicing regulation [25,26]. We restricted the analysis to 197 splicing factors with one-to-one orthologs between mouse and human based on Ensembl BioMart. One-to-one orthologs versus mouse were also retrieved for cow and opossum, when available. All studied splicing factors for the 4 species are listed in S13 Table.

To investigate the enrichment of these splicing factors in trophoblastic tissues in mouse, we focused on 3 comparisons: (i) placenta versus nonplacental tissues; (ii) trophoblast versus ICM at the blastocyst stage; and (iii) TS cells versus embryonic and XEN stem cells. Gene expression quantifications, using the transcript per million (TPM) metric, for public RNA-seq data for different placenta and nonplacenta tissues as well as embryonic blastomeres and ESCs were obtained from VastDB [80], assembly mm10. Short Read Archive (SRA) identifiers and associated information are listed in S1 Table. In addition, for comparison (iii), we generated de novo RNA-seq of different lines of cultured TS cells and XEN stem cells. Culture details and sequencing and mapping information are provided below and in S1 Table. TPM values for TS and XEN samples were obtained using vast-tools v2.5.1 [80].

For comparison (i), for each gene, we first calculated the average TPM expression for each tissue type (groupings are provided in the column “Group” of S1 Table). Next, we calculated the log2 fold change of expression between placental tissues and the average of nonplacental tissues, which are plotted in Fig 1A and shown in S2 Table. Only genes with TPM > 10 in placental tissues were considered. For the other 2 comparisons (Fig 1B and 1C), we followed a similar logic and averaged the expression of each gene in (ii) TE and ICM blastomeres at E3.5 and (iii) TS cells and the average of ESCs and XEN stem cells (ESC and XEN).

To compare placental and nonplacental tissues in human, cow, and opossum, we used a similar approach to the one described for mouse, relying on public RNA-seq data. For human (hg38) and cow (bosTau6), we retrieved TPM values from VastDB [80]. For opossum and nonmammalian vertebrates, we collected RNA-seq data from the SRA (S1 Table) and calculated TPMs using vast-tools v2.5.1 for each species.

Single-nuclei RNA-seq data were obtained from [50]. The Seurat object containing processed and integrated data from all trophoblast nuclei sequenced at different time points (E9.5, E10.5, E12.5, and E14.5) was downloaded from figshare (https://figshare.com/projects/Single_nuclei_RNA-seq_of_mouse_placental_labyrinth_development/92354). Data from individual time points were extracted from the joint Seurat object using the subset function from the Seurat package [81]. For each time point, we generated UMAP plots showing the distribution of the cells of each sample across clusters. Additionally, we generated Mbnl2 and Mbnl3 co-expression plots using the FeaturePlot function with the parameters blend = TRUE, order = T, blend.threshold = 0, and split violin plots using a custom extension of ggplot2 GeomViolin function [82] to generate split violins of the density distribution of Mbnl2 and Mbnl3 expression in each individual cluster.

To identify potential TSSs and alternative splicing events leading to different amino-terminal isoforms of MBNL3, we looked for competing splicing donors (5′ splice sites) that are spliced to the splicing acceptor (3′ splice site) of the second coding exon, which contains the in-frame ATG giving rise to the short protein isoform and that is present in all known Mbnl3 transcripts. For this purpose, we took advantage of vast-tools align, which generates an output file for each RNA-seq sample (sample.eej2) containing the read counts for all exon–exon junctions annotated in the vast-tools junction-based library [80]. The reference splicing acceptor for Mbnl3 for each species were the following: 11 from ENSMUSG00000036109 (mm10), 16 from ENSG00000076770 (hg38), and 0 from ENSBTAG00000014088 (bosTau6). The corresponding donors are numbered starting from 0 (the most upstream) for each species, and the coordinate on the X chromosome and its location in the genome are provided in S2A–S2C Fig and S6 Table. Each donor was grouped and colored based on the evolutionary origin of the associated TSS region and whether or not it would lead to a short MBNL3 protein isoform: black, long isoform from the ancestral TSS; gray, short isoform from the ancestral TSS by skipping of the first coding exon; red, short isoform from the eutherian-specific TSS; and orange, short isoform from the rodent/primate-specific TSSs. To provide the percent of reads coming from each group of donors (Fig 2C, S2E and S2F Fig), we calculated the percent of reads for each of those groups out of the total number of exon–exon junction reads to the reference acceptor as provided in the sample.eej2 files, and, then, we grouped similar samples (e.g., neural tissues, muscle, etc.) as indicated in the column “GROUP” of S6 Table.

Western blot analysis was conducted according to standard laboratory procedures. Lysates were prepared from staged placental samples by disruption with a sterile pestle in Laemmli buffer (0.05M Tris-HCl at pH 6.8, 1% SDS, 10% glycerol, 0.1% b-mercaptoethanol) on ice. Samples were then boiled at 95°C for 10 minutes, sonicated using a bioruptor (10 cycles, 30 seconds on, 30 seconds off), and spun at 20,000 g for 90 seconds. Proteins were separated by 4–15% SDS-PAGE gel (Bis-Tris Criterion TGX pre-cast gels, Bio-Rad) and transferred to PVDF membranes (Thermo Scientific). Membranes were incubated overnight with MBNL3 antisera [46] at 1:500 or anti Lamin B1 (1:1,000, Cell Signaling Technology) in 5% nonfat dairy milk in PBS plus Igepal (0.05%, Merck).

Three fragments of the putative enhancer/promoter region surrounding the eutherian-specific TSS of Mbnl3 were amplified from mouse genomic DNA using the primers shown in S14 Table. The 3 fragments were then cloned together into a single modified pBluescript vector [83] containing a lacZ reporter gene. Prior to microinjection, this construct was linearized and plasmid sequences removed. The construct was then microinjected at 3 to 6 ng/μl into the pronucleus of fertilized mouse oocytes at E0.5. Embryos were collected from E6.5 to E7.5, fixed, and stained for β-galactosidase activity.

Protein sequence alignments and tree production were made using the CLC main work bench (ver. 5.5.1, QIAGEN) software package. MBNL isoforms used for alignment only contained the coding exons common to all Mbnl2 and Mbnl3 orthologs (coding exons 5 [54 bp] and 9 [approximately 74 bp] were excluded from all aligned sequences). Deeply conserved residues (red box, S3A Fig) were defined as those conserved across MBNL2 and noneutherian MBNL3 proteins, whereas all other residues were considered nondeeply conserved (green boxes, S3A Fig). The protein tree was constructed by the neighbor-joining method with Kimura protein corrections. To investigate loss of coding exon 5 (54 bp), Mbnl2 and Mbnl3 orthologs were first examined on the Ensembl (version 102) and UCSC (http://genome.ucsc.edu) genome browsers. For Mbnl genes where exon 5 was not annotated, exon loss was confirmed by manual inspection of the intronic region.

For all Mbnl constructs tested, the exon usage pattern followed that of the predominant MBNL3 isoform in mouse placenta: All MBNL proteins tested included alternative coding exon 8 (95/98 bp) and lacked alternative coding exons 5 (54 bp), 7 (36 bp), and 9 (approximately 74 bp). Mouse Mbnl isoforms were cloned from placenta/limb cDNA using the primers detailed in S14 Table. DNA sequences corresponding to opossum and chicken MBNL isoforms alongside the chimeric MBNLs were synthesized in vitro by Bio Basic and IDT, respectively. Corresponding short isoforms were generated either from genomic DNA or cloned long isoforms using the primers detailed in S14 Table.

The RNA pool generation, RNAcompete pull-down assays, and microarray hybridizations were performed as previously described [39–41]. For the pull down assays, GST-tagged RNA-binding proteins (20 pmoles) and RNA pool (1.5 nmoles) were incubated in 1 mL of Binding Buffer (20 mM Hepes pH 7.8, 80 mM KCl, 20 mM NaCl, 10% glycerol, 2 mM DTT, 0.1 mg/mL BSA) containing 20 mL glutathione sepharose 4B (GE Healthcare) beads (pre-washed 3 times in Binding Buffer) for 30 minutes at 4°C and subsequently washed 4 times for 2 minutes with Binding Buffer at 4°C.

One-sided Z-scores were calculated for all 7-mer motifs as described previously [39]. Z-scores approximate the affinity of the RNA-binding protein for that specific 7-mer. To obtain a measure of the relative contribution of a specific motif (e.g., GCUU or GCAGC) to the overall affinity, we summed the Z-scores of all 7-mers containing the motif among the top 100 7-mers and calculated the fraction of that sum with respect to that of all 100 7-mers (Fig 2G). Similarly, to assess the relative contribution of the nucleotide after a GC, for each 7-mer containing a GC among the top 100 7-mers, we summed the Z-scores for 7-mers containing GCA, GCC, GCG, or GCU and calculated the fraction corresponding to each of them (Fig 2I). Finally, for those 7-mers containing GCxGC or GCxxGC among the top 100 7-mers, we calculated the relative contribution of those with GCxGC or GCxxGC motifs as the log2 ratio between the sums of their Z-scores (Fig 2J).

All protocols were carried out in accordance to the European Community Council Directive 2010/63/EU and approved by the local Ethics Committee for Animal Experiments (Comitè Ètic d’Experimentació Animal-Parc de Recerca Biomèdica de Barcelona, CEEA-PRBB, CEEA 9086 and MDS 0035P2). Mbnl3^tm2.1Sws/+ [49], Mbnl2^tm1.1Sw/+ [37], Mbnl2^tm1Sws/+ [37], and Edil3^{Tg(Sox2-cre)1Amc} [51] mice alongside the required intercrosses were maintained on mixed backgrounds. Embryo genotyping was conducted on tail biopsies as described [84] using published primers. Embryo sex determination was based on the presence or absence of the Y-linked gene Sry as determined by PCR using primers detailed in S14 Table.

In all analyses, WT littermates were used as controls. Furthermore, female Mbnl3 heterozygous embryos carrying a WT maternal allele were considered to be WT for Mbnl3 as it is almost exclusively expressed from the maternal allele in placenta (S5F Fig), and Mbnl2 heterozygotes were considered WT for Mbnl2 as neither placenta nor embryo size was significantly altered in Mbnl2 heterozygotes versus full WT embryos regardless of the KO status of Mbnl3. For the analysis of placenta weights from full KO embryos at E18.5, data were combined from crosses generating just constitutive KO embryos and crosses generating embryos with placenta specific Mbnl3 KO in addition to constitutive KO embryos; therefore, some of the data from these crosses was used in both Fig 3A and S4H Fig. For the analysis of placenta specific Mbnl3 KOs (a genotype only obtainable in female embryos) all results from male siblings of all genotypes were excluded.

For maternal calorie restriction experiments, pregnant dams were fed 50% of their normal daily chow intake from E11.5 until embryo harvesting [71]. Only litters with at least 3 WTs and 3 KOs were used for the analyses. Raw values for placenta and embryo weight for each genotype within each litter are provided in Fig 5A and 5C (calorie restricted diet) and S9B and S9C Fig (normal diet). To normalize for the large variability observed among litters, we performed 2 complementary approaches. First, we plotted the values as empirical cumulative distribution function (ECDFs) to visualize the effects. To test the hypothesis that KO placentas/embryos are heavier than their littermates, we ranked each placenta/embryo within each litter and obtained a total sum of ranks per genotype within each litter. We then randomly swapped the genotype labels within each litter 1,000 times and computed in how many iterations the sum of ranks was the same or smaller than in the test case. Second, we normalized each placenta/embryo measure by dividing its weight by the average of the mean WT and mean KO weights within each litter (Fig 5B and 5D). Statistical significance of the difference between the normalized values was assessed using Wilcoxon rank-sum tests.

Whole-mount in situ hybridization was carried out following standard procedures [85]. Probes were synthesized from cDNAs cloned into the pGEM-T Easy (Promega) plasmid. The cloning primers used for probe generation are shown in S14 Table. Hematoxylin–eosin staining was conducted according to standard laboratory practice on 6-μm sagittal midsections from placentas harvested at E18.5, fixed overnight in 4% PFA and embedded in paraffin wax. Sections from 10 or more placentas were examined per genotype.

TS^JR [27], TS^FxL4 [86], XEN^JR [87], and XEN^DIZ6 [86] cell lines were cultured as described [27,87]. For E11.5 placental samples, the whole embryonic portion of the placentas was separated from the maternal decidua by blunt dissection. For E13.5 and E18.5 placental samples, broad medial sections of placenta were cut and tissue layers separated using forceps. For E13.5 samples, the embryonic portion of the placenta was collected. For E18.5 samples, only the labyrinth zone was collected to avoid signal from the necrotic cells found in the junctional zone of DKO placentas (Fig 3B and 3D) overwhelming the analysis. For sequencing of E11.5 placentas, 3 samples, each containing 3 placentas, were collected per genotype, and samples were paired so that WTs and mutants came from the same litter(s). For sequencing of E13.5 and E18.5 placentas, 2 samples, each containing tissue from a minimum of 5 placentas, were collected for each genotype. In all cases, RNA was extracted using the RNeasy MiniKit (QIAGEN).

For qPCR, cDNA was prepared from total RNA using Superscript III reverse transcriptase (Invitrogen) and random nonamer primers (Invitrogen). qPCR was performed using the NZYSpeedy qPCR Green Master Mix (NZYTech). Primers used are detailed in S14 Table. For analysis of Mbnl2 and Mbnl3 expression, the housekeeping genes Ywhaz and Hmbs were used for normalization. For analysis of Myc expression, the housekeeping genes Ywhaz, Srp14, Gapdh, and Eef2K were used for normalization, and expression was calculated relative to average WT expression.

For gene expression and alternative splicing analyses, standard polyA-selected Illumina libraries were generated using standard protocols by the CRG Genomics Unit. An average of approximately 100 million (XEN and TS cells), approximately 47 million (E13.5 and E18.5) 125-nt paired end, or approximately 47 million 50-nt single end (E11.5) reads were generated using a HiSeq2500 machine. For alternative polyadenylation analysis, the same E13.5 and E18.5 RNA samples were processed by the CRG Genomics Unit with the QuantSeq 3′ mRNA-Seq Library Prep Kit Reverse (REV) from Lexogen following the manufacturer’s instructions. An average of approximately 11 million 50-nt single end reads were generated using a HiSeq2500 machine. General and mapping statistics for all samples are provided in S1 Table.

To obtain exon inclusion levels, we processed RNA-seq samples with vast-tools v2.5.1 [80] following default options. In brief, each sample was mapped using vast-tools align and the mm10 VASTDB library (vastdb.mm2.23.06.20.tar.gz). Next, tables with inclusion levels using the percentage spliced in (PSI) metric for all exons were generated for E13.5 and E18.5 samples separately using vast-tools combine with the default options. Events with sufficient read coverage across all samples within a given time point were then extracted using vast-tools tidy with the following options:—min_N 8—min_SD 0—noVLOW (i.e., only exons with LOW or higher coverage scores were allowed). From these exons, differentially spliced exons between each KO and the WT for each time point were defined as those with an average absolute change in PSI (ΔPSI) between KO and WT ≥15 and at least ΔPSI ≥5 between each of the 2 replicates in the same direction. GO enrichment analyses for genes with differentially spliced exons for each KO versus WT comparison at each time point were performed using DAVID v6.8 [88] for GOTERM_BP_DIRECT, GOTERM_MF_DIRECT, GOTERM_CC_DIRECT, and KEGG_PATHWAY. As a background set for each time point, we used genes with at least 1 exon event annotated in vast-tools with sufficient read coverage across all samples as defined above (9,211 genes for E13.5 and 8,263 for E18.5). Sashimi plots for selected events were generated using ggsashimi [89] with BAM files produced by mapping RNA-seq reads to mm10 genome using STAR [90].

Exons that change during placental development were defined as those with a significant change between WT samples at E13.5 and E18.5 based on the same cutoffs (|average ΔPSI| ≥ 15 and ΔPSI between each replicate ≥5 in the same direction). ESC-differential exons (i.e., exons that are differentially spliced between ESCs and differentiated cell and tissue types) were defined as previously reported [53]. For this purpose, we obtained PSI values from VastDB for mouse (mm10 release 21/12/21) and obtained the average ΔPSI between ESCs and differentiated cell and tissue types using the script Get_Tissue_Specific_AS.pl [91] (https://github.com/vastdb-pastdb/pastdb) with the config file provided in S15 Table and the option—test_tis ESC.

We aligned the Lexogen Quantseq Reverse sequences to the mm11 assembly (Ensembl 104 gene annotation) using STAR with parameters—outFilterMultimapNmax 1—outFilterMatchNminOverLread 0.2—outFilterScoreMinOverLread 0.2—sjdbGTFfile. We then clustered 5-end of alignments with inverted strand (with polyA signal filtering to avoid most internal priming events) to construct the polyA database from our data. We then used this database to compute polyA site counts (clustering of alignments around detected polyA sites with 50 down/up stream nt tolerance) and used DEXSeq [92] to identify regulated polyA sites in genes. We report regulated genes with polyA site pairs that undergo fold changes in opposite directions and for which changes are significant (false discovery rate [FDR] < 0.1).

For further analyses, we selected a single polyA pair per gene and stage (S9 Table). For this purpose, among the pairs (if more than 1) that were significantly misregulated in the largest number of KO versus WT comparisons (from 0 to 3), we performed a ranked selection: (i) we took the pair with the largest absolute fold change if the total number of reads across the 4 conditions for the pair were at least 15% of the pair with the largest number of reads; (ii) else, we took the pair with the largest number of reads of all (not only significant) pairs; (iii) else, if the pair had more than 30% of the number of reads of the pair with the largest number of reads, we took the one with the highest fold change; and (iv) else, we discarded the gene.

RNA-seq samples were used to estimate transcript abundances using Salmon quasi-mapping approach with default parameters plus the—seqBias and the—gcBias options to perform sequence-specific and GC bias correction, respectively, and the—validateMappings option to ensure that mappings are validated using an alignment-based verification [93]. The reference fasta file used for the quantification was the set of cDNA sequences for Ensembl and ab initio genes predicted in the mm10 genome assembly in Ensembl [94] (Version 88). Gene expression estimates were obtained from transcript abundance estimates using the txtimport package [95] with the option countsFromAbundance = lengthScaledTPM, to obtain estimated gene counts scaled using the average transcript length over samples and the library size.

Differential gene expression analysis was performed using these gene expression estimates with the DESeq2 package [96] with default settings. In WT versus KO differential expression analyses, genes with an FDR-adjusted P value (padj) smaller than 0.05 (E13.5 and E18.5 samples) or 0.1 (E11.5 samples) were considered differentially expressed. For display purposes, shrunken log2 fold changes were used. To investigate the effect of the different KOs on placenta maturation, log2 fold changes for each KO at E18.5 were plotted for genes that were differentially expressed between E13.5 and E18.5 in WT samples using strict cutoffs (padj < 0.01 and |FC| > 5) (S8D Fig).

GO analyses for differentially expressed genes for each KO at E13.5 and E18.5 were performed using DAVID v6.8 [88] for GOTERM_BP_DIRECT, GOTERM_MF_DIRECT, GOTERM_CC_DIRECT, and KEGG_PATHWAY. Genes quantified by DESeq2 were provided as background set (17,488 for E13.5 and 17,221 for E18.5). For E11.5, to increase the list of potentially differentially expressed genes, we relaxed the cutoffs to select genes with P < 0.1 and |log2FC)| > 0.1 (409 and 331 up- and down-regulated genes, respectively). We conducted GO enrichment analyses for these subsets of genes using DAVID as described above with all genes for which a valid P value (non-NA) is proved by DESeq2 as background set (16,972 genes). In addition, we used the gene set enrichment analyses for TF binding (ChEA_2016 and ENCODE_TF_ChIP-seq_2015) and the “TF Perturbations Followed by Expression” test from the web tool EnrichR [97].