Phelps WA et al. (2023), Hybridization led to a rewired pluripotency net...

XB-ART-60326

Elife 2023 Oct 03;12. doi: 10.7554/eLife.83952.

Show Gene links Show Anatomy links

Hybridization led to a rewired pluripotency network in the allotetraploid Xenopus laevis.

Phelps WA , Hurton MD , Ayers TN , Carlson AE , Rosenbaum JC , Lee MT .

???displayArticle.abstract???

???displayArticle.pubmedLink??? 37787392
???displayArticle.pmcLink??? PMC10569791
???displayArticle.link??? Elife
???displayArticle.grants??? [+]

Species referenced: Xenopus tropicalis Xenopus laevis
Genes referenced: hes3 homer1 lrrfip1 pou5f3 pou5f3.2 pou5f3.3 sox2 sox3
GO keywords: zygotic genome activation

???attribute.lit??? ???displayArticles.show???

	Figure 1 with 1 supplement Identifying the first wave of genome activation across the two subgenomes. (A) The allotetraploid X. laevis genome contains two distinct subgenomes “L” and “S” due to interspecific hybridization of ancestral diploids. (B) Triptolide inhibits genome activation, as measured in the late blastula, while cycloheximide inhibits only secondary activation, distinguishing genes directly activated by maternal factors. NF = Nieuwkoop and Faber. (C) Heatmap of RNA-seq coverage over exons (left) and introns (right) of activated genes.
	Figure 1—figure supplement 1 Measuring genome activation. (A) (Top) Animal and vegetal views of embryos treated with DMSO (vehicle) versus triptolide. Triptolide-treated embryos fail to gastrulate. (Bottom) Comparison of DMSO versus cycloheximide treated embryo. Treatment was at stage 8, which inhibits progression to stage 9. Scale bar = 0.5 mm. (B) Biplot of RNA-seq for untreated versus DMSO-treated embryos at stage 9, showing no effect on the transcriptome. (C) Biplot of RNA-seq for DMSO versus triptolide treated embryos, showing inhibited activation as detected by exonic (purple) and intronic (green) signal. The predicted mir-427 primary transcript is labeled and exhibits >95% expression inhibition. (D) Histogram of mRNA percent identity between homeolog pairs, as measured by Needleman-Wunsch alignment. Maximum is 0.975 (5 in every 200 bases, which should be generally distinguishable by RNA-seq using 2x100 sequencing reads).
	Figure 2 with 2 supplements Homeologous genes are differentially activated in the early embryo. (A) Proportion of genes encoded as homeologs on both subgenomes versus only one subgenome (singleton) (left), as compared to expression patterns in the early embryo. p Values are from χ-squared tests, 10 d.o.f., comparing genomic to expressed proportions, 16 d.o.f., comparing proportions between activated genes and the maternal contribution, 12 d.o.f., comparing proportions at subsequent stages of activation. (B) Browser tracks showing log2 reads-per-million RNA-seq coverage of equivalently activated homeologs (top) and differentially activated homeologs (L-specific, middle; S-specific, bottom). Trip = triptolide, CHX = cycloheximide. (C) Biplot comparing log2 fold primary activation over triptolide treated embryos for the S homeolog (x axis) versus the L homeolog. (D) Left, proportion of genes activated symmetrically or asymmetrically from the L or S subgenomes, stratified into whether there is a maternal contribution for either homeolog (MZ) or not (Z) (p=0.02, χ-squared test, 4 d.o.f.); and whether a gene is activated only in the stage 9 blastula or is additionally increased in only one or more than one differentiated lineage from stages 10–13 (p=1.3 × 10–22, χ-squared test, 8 d.o.f.). Right, homeolog proportions of later gene activation in epidermal (Epi), neural progenitor (Neur), ventral mesodermal (Meso), and endodermal (Endo) lineages from stages 10–13 (p=3.3 × 10–64, χ-squared test comparing stage 9 and the four lineages, 16 d.o.f.). Lineage-specific gene expression data are from Johnson et al., 2022. (E) Homeolog-specific stage 9 activation proportions, versus maternal contribution homeolog expression patterns, for maternal-zygotic genes. (F) Concordance of homeolog activation patterns across the differentiated lineages at stages 10–13, for genes initially activated at stage 9 and also increased in at least two differentiated lineages. (G) Browser track showing strand-separated reads-per-million RNA-seq coverage over the mir-427 encoding locus on the distal end of Chr1L (v10.1).
	Figure 2—figure supplement 1 Differential homeolog activation over early development. (A) Proportion of genes activated only from one homeolog or from both homeologs at stage 8, as compared to their homeolog activation patterns at stage 9. (B) Biplots showing L dominant (red) or S dominant (blue) activation at stage 8 (left), that resolves to more balanced homeolog activation patterns at stage 9. Each gene is the same color in both plots. (C) Biplot comparing log2 fold all stage 9 activation over triptolide treated embryos for the S homeolog (x axis) versus the L homeolog. (D) Proportion of genes activated symmetrically or asymmetrically from the L or S subgenomes, stratified into whether there is a maternal contribution for either homeolog (MZ) or not (Z), for primary activated genes, all stage 9 activated genes, and stage 10.5 activated genes. p Values are from χ-squared tests, 4 d.o.f. (E) Proportion of genes activated symmetrically or asymmetrically from the L or S subgenomes from stages 10–13 in differentiated lineages. Top plots represent all genes increased at each stage, bottom plots represent genes with maternal contribution <1 TPM. Lineage-specific gene expression data are from Johnson et al., 2022. (F) Proportion of genes encoded on both subgenomes that change homeolog activation patterns in a differentiated lineage in stages 10–13. (G) Concordance of homeolog activation patterns across the differentiated lineages at stages 10–13, for all genes up in at least two different lineages (top) or only genes with maternal contribution <1 TPM (bottom). (H) Boxplots showing log2 L versus S expression ratio for each differentiated lineage over stages 10–13, in the whole transcriptome. N = 8271 homeolog pairs. (I) Zoom of boxplots in (H), demonstrating a gradual bias toward L-dominant expression over time. (J) Significantly (FDR <0.05, Fisher’s exact test, two-sided) enriched Gene Ontology terms in genes activated from both homeologs, as compared to genes activated from only one subgenome. (K) Boxplots of non-synonymous to synonymous substitution rate ratio (dN/dS) shown on a log10 scale, for genes activated from both subgenomes or only one subgenome (p=0.15, Kruskal-Wallis test; median L=0.15, LS = 0.14, S=0.14). All gene groups trend toward stabilizing selection. N = 465, 585, 305 for L, L+S, and S groups, respectively. (L) Boxplots of CDS percent similarity for activation groups as in (K) (p=0.23, Kruskal-Wallis test). For all boxplots: center line, median; box limits, upper and lower quartiles; whiskers, 1.5x interquartile range; points, outliers.
	Figure 2—figure supplement 2 The mir-427 locus. (A) Browser tracks showing strand-separated log2 reads-per-million RNA-seq coverage over the predicted mir-427 primary transcript near the telomere of Chr1L on the v10.1 genome assembly. Xenbase gene isoforms are annotated at the top, aligned precursor mir-427 sequences are annotated in the middle according to strand orientation. (B) Browser track showing log2 reads-per-million RNA-seq coverage over the presumed mir-427 encoding region on the v9.2 genome assembly. The overlapping antisense transcript is not transcribed (all coverage shown is sense to the mir-427 transcript). (C) Dot matrix alignment plot showing BLAST local alignments between v10.1 Chr1L and Chr1S in the region flanking the mir-427 locus. Repetitive sequence alignments (Xenbase soft-masked genomic sequence) are shown in gray, non-repetitive alignments in black. Upstream (dnajb5) and downstream (npr2) homeologous genes are labeled. The L-specific mir-427 locus is highlighted in light red, showing no alignments to Chr1S. (D) Region of v10.1 Chr3S where two additional sequence matches to the mir-427 hairpin are found by BLAT. However, there is minimal RNA-seq coverage, suggesting the Chr1L locus is the only bona fide mir-427 encoding region in the v10.1 assembly. Log2 reads-per-million coverage is shown on the same scale as panel (A).

	Figure 3—figure supplement 1 Profiling homeologous regulatory elements. (A) CUT&RUN for X. laevis blastulae requires cell dissociation prior to nuclear extraction. (B) Comparison of different nuclear extraction techniques. Percent DNA recovered was estimated by NanoDrop quantification of phenol-chloroform extracted DNA after nuclear extraction, as a percentage of theoretical total nuclear DNA mass based on the length of the reference genome sequence. Three nuclear extraction methods were tested with and without cell dissociation: gentle washing by pipeting buffer on the surface of the cells, pipet mixing, vortexing at 1500 rpm. (C) Heatmap of pairwise sample correlation between CUT&RUN samples, as measured by log2 coverage in a 1 Kb window around the center of ATAC-seq open regions (N=58223). (D) Heatmaps showing CUT&RUN coverage over transcription start sites for all samples. (E) Boxplots as in Figure 3B showing CUT&RUN log2 coverage L/S on genes in which only one homeolog is activated at stage 9 (L, S) or both activated (LS), stratified into maternal-zygotic genes (MZ) with >1 TPM maternal contribution (N = 547, 889, 398 for L, LS, and S respectively), and strictly zygotic genes (Z) (N = 46, 126, 35 for L, LS, and S respectively). Two-way ANOVA for each of the four sets indicates there is no significant interaction between MZ/Z and homeolog activation pattern (p=0.37–0.91). (F) Biplots comparing log2 L versus S RNA-seq activation ratio (x axis) and log2 L versus S CUT&RUN read coverage. Linear regression lines are overlaid in red; each has positive slope. p Values are from two-sided Pearson’s correlation tests. (G) CUT&RUN coverage over paired homeologous gene regions around the TSS. Gene pairs are sorted according to L versus S RNA-seq activation ratio (right).


	Figure 4—figure supplement 1 Assessing Pou5f3 and Sox3 roles in genome activation. (A) Embryos injected with 40 ng of each pou5f3.2, pou5f3.3, and sox3 morpholinos have gastrulation defects compared to stage 10.5 control embryos, and fail to close the blastopore. Embryos injected with two of the three morpholinos plus 40 ng of gfp morpholino also had developmental delays, with one of two batches of pou5f3.2+pou5f3.3 morpholino embryos also failing to close the blastopore. Stage 12 phenotype counts are summed over two different rounds of injections on different days. Scale bar = 0.5 mm. (B) Heatmap showing RNA-seq intron log2 fold difference compared to control for exonic and intronic signal, with individual replicates shown. The left columns are for morpholino treatments without cycloheximide, right columns with cycloheximide comared against a cycloheximide control. Triptolide samples are shown for comparison. All samples are from embryos collected when untreated wild-type controls were at stage 9. (C) Heatmap for two morpholino samples that likely had problems with the cycloheximide treatment, based on the RNA-seq patterns. These samples are excluded from subsequent analyses. (D–F) Biplots comparing exonic expression levels in cycloheximide-treated control embryos versus embryos with deficient genome activation. Primary-activated genes with minimal maternal contribution (strictly zygotic) are purple circles, maternal-zygotic genes detected by exonic increases are orange triangles. TPM = transcripts per million. (G–I) Venn diagrams showing overlap of significantly different genes in the morpholino treatments. (J) Proportions of genes classified as significantly up-regulated in morpholino treatments detected by intronic levels. (K) log2 fold difference RNA-seq expression of pou5f3.2+pou5f3.3+sox3 morpholino treated embryos over control for different groups of genes, according to whether they were classified as significantly up-regulated in various morpholino combinations. Although some gene groups achieve significance in the double but not triple morpholino conditions, all gene groups significantly up in at least one treatment nonetheless have elevated expression in the triple morpholino conditions. Center line, median; box limits, upper and lower quartiles; whiskers, 1.5x interquartile range; points, outliers. N = 60, 778, 210, 102, 140, 215, 553, left to right. (L) Proportions of morpholino-affected genes (no cycloheximide) stratified into strictly zygotic genes (Z, maternal contribution <1 TPM) versus maternal-zygotic genes and whether wild-type activation is only in the stage 9 blastula or also observed in differentiated lineages (left, middle). Strictly zygotic genes overall are significantly more likely to be down-regulated by morpholino treatment (p=6.6 × 10–18, χ-squared test, 3 d.o.f.), and this is a stronger effect in blastula-specific genes compared to genes also up-expressed in differentiated lineages (p=2.0 × 10–34, χ-squared test, 6 d.o.f.). (Right) proportions are not significantly different when comparing in which lineages those genes are further expressed (p=0.11, χ-squared test, 6 d.o.f.). (M) Proportions of genes affected by morpholino treatment in the presence of cycloheximide, according to their homeolog expression patterns. Strictly zygotic genes encoded in both subgenomes are more likely to depend on Pou5f3 and Sox3 for their primary activation than singleton genes, regardless of the patterns of homeolog activation (p=0.0020, χ-squared test, 5 d.o.f.). Maternal-zygotic genes have a subtler trend. (M) Proportions of genes affected by morpholino treatment without cycloheximide, according to their homeolog expression patterns. Genes expressed only on the S homeolog, particularly S singletons, are somewhat less likely to be regulated by Pou5f3 and Sox3 (p=3.4 × 10–5, χ-squared test, 10 d.o.f.), though only among maternal-zygotic genes.

	Figure 5 with 1 supplement Regulatory divergence underlies dosage maintenance. (A) Biplots comparing relative expression levels of activated genes in X. laevis and X. tropicalis, treating L and S homeolog contributions separately (middle, right) or summed (left). Individual subgenome expression is scaled 2 x, since transcript per million (TPM) normalization is calculated relative to the entire X. laevis transcriptome. Individual labeled genes are color coded according to the dominant expressed homeolog (red = L, blue = S, purple = equivalent). (B) Barplots showing the proportion of X. laevis genes across homeolog activation categories whose orthologs are also activated in X. tropicalis or part of the maternal contribution. (C) Barplots showing the proportion of X. laevis enhancers across homeolog activity categories that are acetylated in X. tropicalis. (D) Barplots showing the proportion of Xenopus genes whose orthologs are regulated by Pou5f3/SoxB1 and Nanog in zebrafish. Xenopus genes are classified according to how many homeo/orthologs are regulated by Pou5f3/Sox3. Genes with conserved regulation in both X. laevis homeologs and X. tropicalis are more likely to be regulated by Pou5f3/SoxB1 in zebrafish, but also more likely to be regulated by Nanog.

	Figure 6 Model for pluripotency network evolution. X. laevis likely underwent extensive enhancer turnover between its two subgenomes, which nonetheless maintained stoichiometry of pluripotency reprogramming in the early embryo.

References [+] :

Adams, Polyploidy and genome evolution in plants. 2005, Pubmed