Different Genes are Recruited During Convergent Evolution of Pregnancy and the Placenta.

The repeated evolution of the same traits in distantly related groups (convergent evolution) raises a key question in evolutionary biology: do the same genes underpin convergent phenotypes? Here, we explore one such trait, viviparity (live birth), which, qualitative studies suggest, may indeed have evolved via genetic convergence. There are >150 independent origins of live birth in vertebrates, providing a uniquely powerful system to test the mechanisms underpinning convergence in morphology, physiology, and/or gene recruitment during pregnancy. We compared transcriptomic data from eight vertebrates (lizards, mammals, sharks) that gestate embryos within the uterus. Since many previous studies detected qualitative similarities in gene use during independent origins of pregnancy, we expected to find significant overlap in gene use in viviparous taxa. However, we found no more overlap in uterine gene expression associated with viviparity than we would expect by chance alone. Each viviparous lineage exhibits the same core set of uterine physiological functions. Yet, contrary to prevailing assumptions about this trait, we find that none of the same genes are differentially expressed in all viviparous lineages, or even in all viviparous amniote lineages. Therefore, across distantly related vertebrates, different genes have been recruited to support the morphological and physiological changes required for successful pregnancy. We conclude that redundancies in gene function have enabled the repeated evolution of viviparity through recruitment of different genes from genomic "toolboxes", which are uniquely constrained by the ancestries of each lineage.

Introduction

Viviparity (live birth) is one of the most striking examples of convergent evolution in animals: The trait has evolved repeatedly from the ancestral state of oviparity >150 times in vertebrates (Blackburn 2015), and many more times in invertebrates (Ostrovsky et al. 2016). In some pregnant animals, placentas formed by the “intimate apposition or fusion of the fetal organs to the maternal (or paternal) tissues for physiological exchange” (Mossman 1937) have also developed repeatedly. In fact, complex nutritive placentas, which transport large quantities of nutrients to the embryo, have evolved independently at least 16 times across vertebrates (Blackburn 2015; Whittington et al. Forthcoming). Both viviparity and the placenta are among the only complex traits in vertebrates that have evolved repeatedly so many times (Griffith and Wagner 2017; Whittington 2021).

The evolution of viviparity is constrained by the need to meet the biophysical requirements of internally incubated embryos until they complete development. Extensive maternal tissue remodeling is necessary to convert the viviparous uterus from a quiescent, nonreproductive state to a reproductive state capable of supporting an embryo throughout gestation. Therian mammals, viviparous lizards, and viviparous sharks exhibit similar changes in uterine remodeling and expansion of the vascular bed during gestation (Murphy et al. 2000; Biazik et al. 2010; Parker et al. 2010; Wooding et al. 2010; Murphy and Thompson 2011; Buddle et al. 2019). These morphological changes regulate a set of core critical processes for embryonic development: the exchange of respiratory gases, water transport, nutrient provision, and parental immune regulation (e.g., Murphy and Thompson 2011; Van Dyke et al. 2014). Collectively, the morphology and physiology of viviparity and the placenta have evolved convergently across vertebrate lineages, and diverse taxa exhibit remarkably similar uterine and placental structures and functions (Blackburn 2015).

The natural evolutionary replication of viviparity and the placenta enables us to test whether the same genetic changes have driven their convergent evolution in different lineages, or whether each origin of these traits has arisen through novel genetic mechanisms. Both scenarios are possible: although studies of some convergent phenotypes (e.g., eusociality, electric organs, bioluminescent organs) have revealed convergent genetic changes (e.g., Woodard et al. 2011; Gallant et al. 2014; Pankey et al. 2014, respectively), others have determined that convergent molecular changes producing the same phenotype are rare (e.g., Foote et al. 2015; Zou and Zhang 2015). We can make predictions for our traits of interest (viviparity and the placenta) based on qualitative studies that have identified shared genes underpinning independent origins of pregnancy (e.g., Murphy and Thompson 2011; Brandley et al. 2012; Griffith et al. 2016; Carter 2018; Foster et al. 2020; Recknagel et al. 2021; van Kruistum et al. 2021). Strikingly, similarities are evident across very divergent taxa separated by 450 million years of evolution (Irisarri et al. 2017), including amniotes and anamniotes (e.g., Schartl et al. 2013; Whittington, Griffith, et al. 2015; Guernsey et al. 2020; Whittington and Friesen 2020; Parker et al. 2021; Recknagel et al. 2021; van Kruistum et al. 2021; Du et al. 2022). However, these studies do not quantitatively test whether more genes are shared across pregnant taxa than would be expected by chance alone, an issue that we address here.

Here, we evaluate the changes in uterine gene expression that occur during independent origins of pregnancy across the vertebrate phylogeny. We predict that there should be significant overlap of gene use across independent origins of vertebrate viviparity and the placenta, particularly given that many relatively phylogenetically distant species use the same gestational tissue (the uterus) (Lombardi 1998). We compare the gestation response in five lizard species that differ in both parity mode (oviparous vs. viviparous) and nutritional mode (lecithotrophy: embryos predominantly rely on yolk for nutrition; placentotrophy: embryos receive substantial nutrition via a placenta), two viviparous mammals (one metatherian, one eutherian), and one viviparous shark. Each of these species was carefully chosen because all gestate embryos within homologous tissues, at least temporarily, in a muscular section of the female reproductive tract (derived from the Müllerian duct, Lombardi 1998, termed the uterus [also known as the oviduct]). In animals with complex nutritive placentas (placentotrophic lineages), uterine tissues form the maternal portion of the placenta. Our comparisons reveal that, as expected, many of the same processes supporting gestation occur across viviparous lineages. Additional shared functions occur across placentotrophic lineages. Strikingly, though, we find that there are no orthologous genes upregulated in common across viviparous lineages. Instead, convergent gestation-related functions that are shared across viviparous lineages are produced by different genes. Our results show that, when considering viviparity, convergent evolution results in similarities in morphology and physiology that are predictable due to their shared functions, yet are produced by unpredictable divergences in gene expression.

Results and Discussion

The genetic basis of convergent evolution of viviparity and the placenta can be determined by identifying the changes in uterine gene expression of different species across gestation, and then testing for overlaps between species. To do this, comparable datasets are required, generated from the same tissues (uterus) and reproductive stages (see below) across species. To this end, we assembled high-quality uterine transcriptomes de novo for two oviparous lizards (Bassiana duperreyi, Lampropholis guichenoti), one viviparous lecithotrophic lizard (Niveoscincus coventryi), two viviparous placentotrophic lizards (Niveoscincus ocellatus, Pseudemoia entrecasteauxii), and a viviparous placentotrophic shark (Rhizoprionodon taylori) (supplementary table S1, Supplementary Material online), and leveraged uterine transcriptomes and reference genomes for two placentotrophic mammals (Monodelphis domestica, Rattus norvegicus). These species represent four independent origins of viviparity and the placenta. All of our de novo transcriptome assemblies are of high quality, with BUSCO completeness scores ranging from 89.6% to 94.8% (median: 93.3%), and read alignment rates ranging from 84.3% to 94.2% (median: 93.2%) (supplementary table S1, Supplementary Material online). High-quality assemblies are important for our study, ensuring that any lack of signal for convergent gene expression evolution across lineages is real, rather than an artifact of poor data quality.

Convergence in Function, Not Gene Recruitment

We tested for changes in uterine gene expression between the nongestating and gestating states (here, uterus late in pregnancy, when placental species are transporting significant quantities of nutrients) for each species separately, using differential expression and gene ontology (GO) term enrichment analyses (supplementary tables S2 and S3, Supplementary Material online). GO term enrichment analyses provide an approximation of the functional role of genes that are expressed differently across gestation within each species. We cross-referenced the functional hypotheses derived from GO term enrichment analyses with functional changes known to occur in our study species based on morphological and physiological experiments (e.g., Murphy et al. 2000; Carter 2012; Van Dyke et al. 2014; Blackburn 2015; Buddle et al. 2021). Most species exhibited a strong shift in gene expression profile between nongestation and gestation, with hundreds to thousands of up/downregulated genes, enriched for a broad suite of GO terms (supplementary table S3, Supplementary Material online). The only outlier species was the oviparous lizard L. guichenoti, where only 269 genes were differentially expressed, many of which were not annotated despite an overall robust, well-annotated transcriptome assembly (supplementary table S1, Supplementary Material online). The differential expression results for L. guichenoti were not enriched for any meaningful functional categories (supplementary tables S2 and S3, Supplementary Material online), demonstrating a lack of consistent gene expression response to bearing an embryo among biological replicates.

We next tested whether the GO terms enriched in differentially expressed gestation-related genes correlate with reproductive mode and the presence of a placenta. We consider overlaps in enriched GO terms to represent convergence in physiological functions occurring in the uterus during gestation. Using this approach, we found no enriched GO terms shared across all species (fig. 1, supplementary table S4, Supplementary Material online), suggesting that gestation-related functions are different between oviparous and viviparous species. Importantly, we found 23 enriched GO terms shared by all placentotrophic viviparous species, and 16 enriched GO terms shared by all viviparous species (fig. 1, supplementary table S4, Supplementary Material online). Enriched GO terms shared by viviparous species largely relate to transport functions, metabolism, chemical homeostasis, and the apical plasma membrane (supplementary table S5, Supplementary Material online). The enriched GO terms shared by placentotrophic species were similar, but included an additional seven enriched biological processes, including the transport of lipids, anions, organic acids, and sodium ions. A generalized Fisher’s exact test revealed that the enriched GO term overlaps between all placentotrophic species, and all viviparous species, are significantly greater than expected by chance alone (fig. 1, supplementary table S4, Supplementary Material online). These results provide gene expression evidence that many of the same broad physiological processes occur in all viviparous species during gestation, which supports results from morphological and physiological studies (e.g., Murphy et al. 2000; Thompson et al. 2000; Murphy and Thompson 2011; Carter 2012; Van Dyke et al. 2014; Blackburn 2015; Dudley et al. 2017; Buddle et al. 2019, 2021). We explore the functional significance of these enriched GO terms, including why the biological processes are likely to be important for viviparity, in our investigation into the genes driving the GO enrichment in each species (see below).

Fig. 1.

UpSet plots visualizing the degree of overlap in enriched GO terms (a and b) and differentially expressed genes (c and d) across eight vertebrate species, taking into account parity mode and nutritional mode (placentotrophic vs. lecithotrophic embryonic nutrition). Vertical bars indicate overlaps that are greater than expected by chance alone, as determined using generalized Fisher’s exact tests. Horizontal bars represent the total number of enriched GO terms or differentially expressed orthogroups per species.

However, similarity in function does not necessarily imply similarity in gene recruitment, so we next tested whether our species achieve the same physiological outcomes by differentially expressing the same genes. We inferred groups of orthologous genes (orthogroups) to allow meaningful cross-species comparisons to identify the same genes. An orthogroup is “the set of genes that are descended from a single gene in the last common ancestor of all the species being considered”, so contains both orthologues across species and paralogues within a species (Emms and Kelly 2015). We considered an orthogroup to be differentially expressed in a species if any of its constituent genes were differentially expressed.

We tested for overlaps in differentially expressed orthogroups within all placentotrophic viviparous species, all viviparous species, and all species (both oviparous and viviparous). Henceforth, “upregulated”/“downregulated” refers to genes with a higher/lower expression in the gestating uterus compared with the nongestating uterus. Surprisingly, only one orthogroup was upregulated in common across all placentotrophic viviparous lineages (fig. 1, supplementary tables S6 and S7, Supplementary Material online). This orthogroup contains a single cathepsin gene, putatively functioning in intracellular protein catabolism, which is necessary for tissue remodeling (Coulombe et al. 1996). No orthogroups were upregulated or downregulated in common across all viviparous species, or across all species. Nine orthogroups were downregulated in common across both oviparous species; these orthogroups predominantly comprised immunity-related genes such as those coding for major histocompatibility complex proteins and immunoglobulin-binding proteins (fig. 1, supplementary tables S6 and S7, Supplementary Material online). Generalized Fisher’s exact tests revealed that the few overlaps in differentially expressed orthogroups were not greater than expected by chance alone (fig. 1, supplementary table S6, Supplementary Material online). Therefore, contrary to prevailing assumptions about the genetic basis of viviparity, our results demonstrate no statistically significant convergence at the gene expression-level across all viviparous lineages in our study, despite the convergence of their uterine morphology and physiology.

To test whether any particular species may be disproportionately driving these results, we iteratively removed one species at a time from the analysis and reexamined overlaps in differentially expressed genes across viviparous species. By this less conservative approach, we still found very little overlap in orthogroup recruitment. For example, when R. norvegicus was removed from the comparison, there was only a maximum of eight upregulated orthogroups and one downregulated orthogroup overlapping across all remaining viviparous species (supplementary table S8, Supplementary Material online). Even then, the differentially expressed genes that overlap between viviparous species in each permutation do not fulfill all physiological functions that we previously identified as being critical for viviparous pregnancy (supplementary tables S3 and S9, Supplementary Material online); thus, some functions must be controlled by different differentially expressed genes across species. Therefore, our tests reveal an overall striking lack of convergence in the genes that are differentially expressed during gestation in viviparous species or viviparous placentotrophic species, despite a strong similarity in their uterine functions.

Pregnancy Relies on Redundancy in Gene Functions

Since the viviparous uterus exhibits convergence in function and morphology across species, which is controlled by different genes, we predict that there must be redundancy in uterine gene function across species. Our overlap comparison of enriched GO terms is a higher-order approximation of the functional changes occurring within the uterus of these species. To further determine how gene redundancy might influence patterns of differential expression, we also explored the differential expression results to identify orthogroups involved in the core set of critical gestational functions (discussed in Van Dyke et al. 2014): nutrient transport, respiratory gas exchange, water transport, and immune regulation (supplementary table S10, Supplementary Material online).

Embryonic development requires a supply of amino acids (AA) (Avagliano et al. 2012), some of which are transported across the placenta in placentotrophic species. The solute carrier (SLC) gene family fulfills essential transport functions, including the transport of AA. These genes are broadly expressed and are highly redundant in the specific substrates they transport (Schumann et al. 2020). They, thus, represent an ideal gene family to identify convergent gene recruitment across independent origins of the placenta. Therefore, we identified all AA-transporting SLCs (henceforth: AA-SLCs) upregulated in the gestating stage relative to the nongestating one. We found a remarkable divergence in AA SLC recruitment during gestation among species.

We detected 84 AA SLCs in the uterine transcriptomes of our species, of which 31 were detected in all species, and 44 were upregulated in at least one lineage (fig. 2). All species, except for L. guichenoti, upregulated at least one AA SLC. Additionally, at least one upregulated AA SLC was species-specific, in every viviparous species, except R. norvegicus. For example, out of a total of 22 AA SLCs upregulated in gestating M. domestica, 13 were upregulated only in this species (fig. 2). There were no AA SLCs upregulated in common across all species, all viviparous species, or all oviparous species. We predicted that there would be a strong upregulation of AA SLCs in animals with a complex nutritive placenta, because AA transport is an important function of this organ. Although we did observe this pattern, we found no overlap in AA SLCs across our placentotrophic species. Thus, although AA SLC recruitment must be important for pregnancy in placentotrophic viviparous species, which individual AA SLC family member is recruited differs between lineages.

Fig. 2.

The phylogenetic relationships among the eight vertebrate species in our study, as well as their recruitment of AA-transporting SLC genes (AA SLCs) during gestation. All orthologous AA SLCs detected in the transcriptomes of our species are listed, and the recruitment of AA SLCs from this “toolbox” to different lineages is demonstrated using colored boxes. AA SLCs with orange boxes are upregulated by more than one lineage, and those with red boxes are specifically upregulated in only one lineage. Image credit for opossum silhouette: Sarah Werning (CC BY 3.0). V, viviparous; P, placentotrophic; L, lecithotrophic; O, oviparous.

One AA SLC was shared by all viviparous species, except R. norvegicus: SLC38A3. SLC38A3 encodes a membrane-bound protein that mediates sodium-coupled cotransport of glutamine, histidine, asparagine, and alanine (Fei et al. 2000), and is implicated in energy metabolism, ammonia detoxification, and gestation (Rubio-Aliaga and Wagner 2016). Presumably, SLC38A3 contributes to important gestation-related functions in most of our viviparous species. This result raises the question: is R. norvegicus able to carry out the same physiological functions during gestation as the other viviparous species, despite not upregulating SLC38A3? To address this question, we compared all functions assigned to SLC38A3 based on annotated GO terms, to the functions for all AA SLCs upregulated in R. norvegicus. Crucially, we found that most functions normally carried out by SLC38A3 are also carried out by (1) other AA SLCs and (2) genes in other families, which are upregulated in R. norvegicus (supplementary table S11, Supplementary Material online). The same is true of SLC3A2: all viviparous species, except R. norvegicus, upregulate SLC3A2, but other genes upregulated by R. norvegicus are implicated in most of the same physiological functions as SLC3A2 (supplementary table S11, Supplementary Material online). Collectively, our AA SLC results demonstrate that different species predominantly recruit different genes that drive the same function of AA transport, with no correlation between parity mode or nutritional mode (placentotrophic vs. lecithotrophic embryonic nutrition) and AA SLC recruitment. Instead, redundancy in gene function accounts for differences in differentially expressed genes among lineages.

All species differentially express genes implicated in respiratory gas exchange (supplementary table S10, Supplementary Material online), which is critical for embryonic development (Andrews and Mathies 2000). No orthologous genes related to respiratory gas exchange are, however, differentially expressed in common by all species, all viviparous species, or all oviparous species. Morphological evidence suggests that parental regulation of respiratory gas exchange during gestation is often enabled by increased gestational vasculature via angiogenesis (Murphy et al. 2010; Parker et al. 2010; Dudley et al. 2021). All eight species differentially express different genes implicated in the angiogenesis process, implying that dynamic upregulation and downregulation of a suite of angiogenic genes is necessary in vertebrate gestation (supplementary table S10, Supplementary Material online). During gestation, oxygen is transported to the fetus to allow fetal metabolism and growth (Carter 2000). In viviparous animals that develop a placenta, fetal excretion of CO2 occurs via the placenta, with morphological changes to the maternal uterus reflecting maternal participation in the excretion process (Van Dyke et al. 2015). The dynamic role of the maternal uterus in providing gas exchange during gestation is clear, with all eight species examined here differentially expressing different genes involved in the biological response to hypoxia. Additionally, all species, except L. guichenoti, upregulate different genes that contribute to a biological response to carbon dioxide, mostly via carbonic anhydrase orthologues. Thus, the same requirements for respiratory gas exchange during gestation can be met through recruitment of different genes.

Another critical requirement for embryonic development is the provision of adequate water. All species in our study differentially express many genes that presumably contribute to water provision to embryos, through roles including the maintenance of water homeostasis (supplementary table S10, Supplementary Material online). Some genes differentially expressed by our species have previously been linked to the active and passive transport of water during gestation (aquaporins, SLCs) (Lindsay and Murphy 2007). The suite of differentially expressed genes controlling these functions is variable among species. From the same background set of orthologous genes that contribute to water provision, there is no single route to controlling this process during gestation; different lineages use different combinations of genes that contribute to these physiological functions (e.g., aquaporins, SLCs, hyaluronidases).

Many immunological studies of the evolution of viviparity have focused on why the mother’s immune system does not reject the developing embryo, which is semiforeign genetically (Medawar 1953). All species in our study differentially express genes that are implicated in the adaptive and innate immune responses, including genes coding for zinc finger proteins and cell surface ligands (supplementary table S10, Supplementary Material online). All viviparous lineages, except for N. ocellatus, upregulate an immune-related orthogroup containing C–X–C motif chemokine ligands, also known as interleukins or inflammatory cytokines (CXCL1, CXCL2, CXCL3, CXCL8, CXCL13, PPBP). Interleukin-family cytokines are important for the prevention of inflammatory damage to pregnant human mothers (Southcombe et al. 2015) and to regulate neutrophil trafficking to successfully resolve tissue injury (Sawant et al. 2016). Extensive differential expression of interleukins during gestation also occurs in squamate lizards (Biazik et al. 2007; Brandley et al. 2012), amphibians (Jantra et al. 2007), syngnathid fishes (Whittington, Griffith, et al. 2015), and chondrichthyans (Cateni et al. 2003). Instead of CXCL genes, N. ocellatus upregulates a suite of other genes that function in adaptive and innate immune responses (supplementary table S9, Supplementary Material online). Collectively, both the present study and previous research suggest a highly conserved role for inflammatory cytokines in viviparous gestation, but some lineages use alternative genes to control immune regulation during pregnancy.

Gene Expression Largely Reflects Phylogeny, Not Reproductive Phenotype

We estimated a species tree based on individual trees for all orthogroups (supplementary fig. S1, Supplementary Material online), which agrees with recent phylogenomic estimates of the relationships among our species (Brandley et al. 2015; Irisarri et al. 2017). Trees can also be estimated based on distance matrices formed from gene expression data. In some cases of convergent evolution, these expression-based trees reflect shared phenotypes rather than evolutionary history (Stern and Crandall 2018). We searched for convergence at the gene expression-level across lineages with the same parity mode or nutritional mode by inferring the relationships among our species based on orthogroup expression estimates (see Materials and Methods). Using the same approach, we also observed whether lineages clustered based on the stage of gestation, which would suggest a convergent evolution of a “gestational expression phenotype.”

The relationships we inferred in our phylogeny based on gene expression (fig. 3) largely agree with our phylogenetic trees inferred from sequence data (supplementary fig. S1, Supplementary Material online), with the only differences being between closely related lizard species. Principal components analysis of normalized expression estimates also recovers strong clustering based on phylogenetic relationships rather than reproductive phenotype (fig. 3). We do not observe clustering of lineages based on parity mode or nutritional mode (e.g., all viviparous lineages forming a clade, or all placentotrophic lineages forming a clade), and we do not recover interspecific clustering based on the stage of gestation (nongestation vs. late-gestation). These results reinforce our conclusions that there is no clear convergence in gene expression in lineages with the same parity/nutritional mode.

Fig. 3.

Assessing the extent of convergence in gestation-related differential expression of orthologous genes across eight vertebrate species. (a) Phylogenetic tree for all species in our study, inferred using neighbor-joining analysis of a distance matrix of orthologous gene expression patterns. This “expression tree” is largely congruent with estimates of the species tree based on sequence data. (b) PCA of cross-sample normalized expression values for the eight species. The color of shapes is specific to each species, and the type of shape (circle vs. triangle) demonstrates the stage of gestation for each replicate sample (nongestating vs. gestating).

The percentage similarity of overlap in differential expression results, as assessed through pairwise comparisons of species, does not strictly follow phylogeny (supplementary fig. S2, Supplementary Material online). In general, however, overlaps in differential expression patterns are the closest in species that are the most closely related. For example, lizard species tend to differentially express more of the same genes in common than they do with mammals or the shark, and the highest percentage overlap in upregulated orthogroups occurs between two lizard species in the same genus, despite their differences in the mode of nutritional provision to embryos (supplementary fig. S2, Supplementary Material online). Overall, our results presented in a phylogenetic context reinforce our assertion that the similarities in gestation-related physiological functions between viviparous species are not driven by differential expression of the same genes.

Different Genes Underpin Convergent Origins of Pregnancy

Convergent evolution of phenotypes is common in nature, and the genetic basis of this convergence is a key question in evolutionary biology. Our study uses the repeated origins of viviparity and the placenta as a powerful, quantitative test of genetic convergence of these traits, across major vertebrate groups separated by hundreds of millions of years of evolution. We found similarities in the physiological functions of the genes differentially expressed during gestation across viviparous lineages. However, contrary to our expectations, these functional similarities are driven by recruitment of different genes. These results are surprising because previous qualitative studies have identified similar genes being expressed in the pregnant uteri of relatively unrelated taxa (Brandley et al. 2012; Schartl et al. 2013; Griffith et al. 2015, 2016; Gao et al. 2019; Foster et al. 2020; Guernsey et al. 2020; Parker et al. 2021; Recknagel et al. 2021; Du et al. 2022), but we quantitatively demonstrate here that such similarities are not reflective of a broader convergent gene recruitment across the diversity of viviparous animals. Our study shows that successful gestation is not driven by the evolution of the same genes in evolutionarily distant species. Instead, the expression of different genes from ancestral orthologous “toolkits” has been uniquely modified in each lineage to achieve gestation-related functions that are convergent in both morphology and physiology.

The remarkable convergence in uterine remodeling and functional changes during vertebrate gestation, despite different genetic starting points, is probably a result of biophysical constraints on internal embryonic development. The strong phylogenetic signal of the gene expression response during gestation observed in this study suggests that gene recruitment for this trait is tightly constrained by evolutionary history. Overall, our results demonstrate that different genes can be recruited to drive independent origins of viviparity and the placenta across vertebrates.

Materials and Methods

All animal work in this study was conducted according to protocols approved by the University of Sydney Animal Ethics Committee (approval numbers 668, 2018/1405) and the James Cook University Animal Ethics Committee (approval numbers A1933, A2310). We chose eight species from across the vertebrate phylogeny to obtain independent origins of viviparity, as well as independent origins of nutritional mode. The two nutritional modes were lecithotrophy (embryos predominantly rely on yolk for nutrition) and placentrotrophy (embryos receive substantial nutrition via a placenta) (fig. 2, supplementary table S1, Supplementary Material online). These species are the Australian sharpnose shark (R. taylori), the southern grass skink (P. entrecasteauxii), the eastern three-lined skink (B. duperreyi), the common garden skink (L. guichenoti), the ocellated cool-skink (N. ocellatus), the southern forest cool-skink (N. coventryi), the gray short-tailed opossum (M. domestica), and the brown rat (R. norvegicus). Across all species, we sampled tissue from individuals at a late stage of gestation (during which a placenta is present in placentotrophic lineages) and at a nongestating stage to investigate gestation-specific changes in gene expression. We ensured that all gestating animals were at functionally equivalent stages of gestation to avoid confounding of our results, that is, late in gestation for all species, and in all placentotrophic species, at a stage where large quantities of nutrients are being transported to the embryos (Dufaure and Hubert 1961; Simpfendorfer 1992; Abrahamsohn and Zorn 1993; Mate et al. 1994). We specifically sampled tissues of the uterus (a muscular tube derived from the Müllerian duct, proximal to the cloaca/vagina) in all animals in our study, avoiding the Therian specialized region termed the vagina (Wagner and Lynch 2005), to ensure that the tissues were homologous in all species.

Data for P. entrecasteauxii, L. guichenoti, and M. domestica (Griffith et al. 2016, 2017) were downloaded from the Sequencing Reads Archive (accessions: SRR1200828, SRR1201775, SRR1201776, SRR1201777, SRR1201787, SRR1201788, SRR1201789, SRR1201790, SRR4293353, SRR4293354, SRR4293359, SRR4293360, SRR4293361, SRR4293362, SRR5822131, SRR5822132, SRR5822133, SRR5822134, SRR5822135, SRR5822136). For all other species, we sampled animals from the wild populations or from lab-bred populations. Collection permits include the following: 185560 (Queensland government), G15/37987.1 (Great Barrier Reef Marine Park Authority), SL100401 (NSW government), and FA18234 (Tasmanian government).

Tissue Sampling

The placenta comprises closely apposed embryonic and maternal tissues (Mossman 1937). In all pregnant animals, we sampled the maternal side of the placenta (uterine tissue). In all nonpregnant animals, we sampled the equivalent tissue (empty uterus). We specifically sampled oviductal regions excluding the vagina/cloaca, which in pregnant individuals are closely apposed to the embryonic tissue, forming a morphologically distinct placenta.

We collected females of R. taylori from populations ca. 1 km offshore from Townsville, QLD, Australia. The breeding season of R. taylori is highly cyclical (Simpfendorfer 1992), allowing a precise estimation of the reproductive stage. Prior to tissue sampling, all sharks were euthanized by severing the spinal cord with a sharp knife immediately after capture. We sampled tissue from the uterus of nonpregnant sharks and from the uterine portion of the placenta of late-pregnant sharks.

Gestating females of B. duperreyi were collected from field sites in Brindabella National Park, NSW, Australia, and we collected gestating females of N. coventryi from field sites in Kanangra-Boyd National Park, NSW, Australia, and gestating females of N. ocellatus from field sites near Orford, Tasmania, Australia. We maintained these lizards in captivity according to standard protocols (Linville et al. 2010) until they reached the late gestating stage of the reproductive cycle, according to an established embryonic staging system (50). Lizards were then euthanized by injection with 0.1 mL of sodium pentobarbital (6 mg/mL) as previously described (Whittington, Grau, et al. 2015). In each case, we excised tissue from the chorioallantoic placenta and omphaloplacenta into a sample of total placenta tissue. We also obtained nonreproductive replicates by processing lizards 4–6 weeks after parturition/oviposition.

Wistar laboratory rats were sampled from a captive breeding colony at the University of Sydney. Once at the desired stage of the reproductive cycle, rats were euthanized via lethal injection of pentobarbitone, and uterine tissues were sampled. Using this strategy, we obtained nonpregnant and late-pregnant replicates.

For a summary of sampling numbers and RNA quality across all species, see supplementary table S13, Supplementary Material online. In all cases, tissues were excised and fixed in RNA later for 24 h at 4 °C and then stored at −80 °C. We then extracted total RNA from tissues using an RNeasy Plus Mini Kit (Qiagen, Hilden, Germany), including an on-column DNAse digestion step (RNase-free DNase set, Qiagen). We also extracted RNA from some tissues using the Direct-zol RNA MiniPrep Plus kit (Zymo Research, CA, USA). RNA integrity was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), and library preparation was conducted using a TruSeq mRNA stranded library preparation kit (Illumina, San Diego, CA, USA). Finally, samples were sequenced either on the NextSeq 500 rapid run platform (llumina, San Diego, CA, USA), with 75 bp paired-end sequencing, or on the NovaSeq 6000 platform (llumina, San Diego, CA, USA), with 100 bp paired-end sequencing, at the Ramaciotti Centre for Genomics, Sydney, Australia.

Transcriptome Assembly, Annotation, and Abundance Estimation

We ensured that sequencing reads from all samples passed standard quality control tests implemented in FastQC (Andrews 2010) and removed any residual sequencing adaptors using BBDuk and error-corrected reads using Tadpole (both available from https://sourceforge.net/projects/bbmap/). With the exception of M. domestica and R. norvegicus, all species in our study lack a reference genome; therefore, we assembled transcriptomes de novo. Reads for P. entrecasteauxii were originally separated into chorioallantoic placenta and omphaloplacenta upon download from the Short Reads Archive. However, we combined all reads from both types of placenta into a single placenta sample for each P. entrecasteauxii individual to match sampling of the other lizard species. Each of our samples was sequenced to a sufficient depth to require in silico normalization of reads prior to transcriptome assembly. Therefore, we normalized the reads for each sample to an average depth of 40× coverage using BBNorm (available from https://sourceforge.net/projects/bbmap/).

In all species lacking a reference genome, we carried out de novo assembly of normalized reads using four different assemblers: (1) Trinity v2.8.6 (Grabherr et al. 2011), with all parameters set to default values, except for the –no_normalize_reads flag; (2) Trans-ABySS v2.2.3 (Robertson et al. 2010), assembled with four different k-mer sizes (k = 25, 35, 45, and 55), followed by merging with transabyss-merge; (3) rnaSPAdes in SPAdes v3.13.0 (Bushmanova et al. 2019), assembled with two different k-mer sizes (k = 33, 49), followed by automatic merging; and (4) SPAdes single-cell mode in SPAdes v3.13.0, with three different k-mer sizes (k = 21, 33, and 55), followed by an automatic merging of scaffolds into an assembly. For each species, we then combined all four transcriptome assemblies into a single overassembly, then reduced the redundancy of transcripts using the EvidentialGene pipeline (Gilbert 2019). The EvidentialGene pipeline also provides accurate protein reconstructions for all isoforms of assembled genes. We assessed the completeness of our transcriptome by searching against the Vertebrata (ODB9) database to find highly conserved orthologues using BUSCO v3 (Waterhouse et al. 2018). The input sequences for BUSCO were predicted peptide sequences obtained using the EvidentialGene pipeline. For additional assembly quality assessment, we calculated summary statistics, including the percentage of annotated transcripts, the read alignment rate using Salmon v1.2.1 (Patro et al. 2017), and the median transcript length (supplementary table S1, Supplementary Material online).

We annotated all de novo transcriptome assemblies by conducting homology searches against the SwissProt peptide database (supplementary table S12, Supplementary Material online). For all homology searches, we used MMseqs2 Search (Steinegger and Söding 2017) with the sensitivity flag (-s) set to 6, which has similar sensitivity to NCBI BLAST+ but is much faster. Searches with an E-value ≤ 10−5 were considered to be significant. We leveraged GO term annotations for each gene based on SwissProt annotations.

We estimated the abundance of transcripts using the alignment-free approach of Salmon v1.2.1, with automatic library-type detection, correction of potential fragment-level GC bias, and mapping validation to account for potential spurious mapping artifacts. For our species that lack a reference genome, we quasi-mapped reads against the indices of our de novo assemblies. For the two species that have reference genomes available (R. norvegicus: Rnor_6.0.99 and M. domestica: ASM229v1.100), we concatenated the reference genome and transcriptome into a single file, built a decoy-aware index with Salmon, and then quasi-mapped mapped processed reads against their respective indices.

Differential Expression and Identification of Orthologous Genes

We tested for differential expression between nongestating and gestating individuals within each species using DESeq2 (Love et al. 2014). We chose not to directly test for differential expression between oviparous and viviparous species, such as by directly contrasting gene expression in a gestating viviparous species with gene expression of a gestating oviparous species. Such an approach might identify genes that have a different expression across species, but would not allow any differences or similarities to be linked to gestation. Additionally, by focusing on how the expression of genes changes during gestation within a species, we are able to identify genes that are important for successful gestation in different lineages, rather than those genes that play a “housekeeping” role regardless of whether an animal is gestating or not.

For each species, we summarized counts for each transcript to the gene level using the tximport package in R (Soneson et al. 2015). We carried out minimal prefiltering of genes to remove those genes for which reliable expression estimates could not be estimated (raw counts < 10). We accounted for any latent, unwanted sources of variation unrelated to our predictor variable (nongestating vs. gestating) using surrogate variable analysis in the sva R package (Leek et al. 2012). Any genes with an absolute log fold change >1 and adjusted P-value (false discovery rate) <0.05 were considered to be differentially expressed. We also tested for functional enrichment of GO terms in differentially expressed genes using the goseq R package (Young et al. 2010), which employs overlap enrichment analyses while accounting for potential gene-length biases. Each species had a custom background gene set against which enrichment was tested, consisting of all genes that passed prefiltering of raw gene counts and independent filtering by DESeq2.

We aimed to compare the expression of the same genes in all species, but comparing species separated by such a great evolutionary distance raises issues of orthology. For example, the presence of unique genes in any lineage, or duplication of genes in any lineage(s), complicates such a comparison. Therefore, we searched all transcriptomes for groups of orthologous genes (“orthogroups”) using OrthoFinder v2.3.12 (Emms and Kelly 2015). For the purpose of orthology searching, we selected the transcript with the longest predicted peptide sequence for each gene. In total, we identified 48,684 orthogroups, including 5,591 orthologues that were single-copy in all eight species. This number of single-copy orthologues is fewer than would be expected from a genomic dataset because transcriptomic datasets only contain genes that are expressed in the target tissue at the time of tissue sampling.

Using the output of the OrthoFinder analysis, we mapped differentially expressed to their corresponding orthogroups. We then classed an orthogroup as differentially expressed within a species if any of its constituent genes were differentially expressed. This approach does not assume that all paralogues of a gene within a species will be differentially expressed, reflecting knowledge that expression profiles vary between gene copies (Kegel and Ryan 2019). We visualized patterns in our expression data using principal components analysis of variance-stabilization normalized expression data for all orthogroups for which all species have at least one gene represented.

We then compared the differential expression results across species by searching for overlaps in differentially expressed orthogroups. We focused on testing four specific comparisons: (1) all species, (2) all oviparous species, (3) all viviparous species, and (4) all placentotrophic viviparous species. For each comparison, we determined all orthogroups that were differentially expressed in all constituent species, as well as those orthogroups that were specifically differentially expressed in all constituent species. We followed the same procedure looking for overlaps in enriched GO terms. Finally, we investigated the upregulation of genes from the SLCs family as a case study of redundancy in gene function.

An important consideration when searching for overlaps is the possibility that we might find differentially expressed orthogroups to be overlapping between species based on chance alone. Therefore, we tested for significance of the overlaps within our groups of interest using a generalized Fisher’s exact test using the SuperExactTest R package (Wang et al. 2015). For a particular comparison (e.g., all viviparous species), given (1) the sets of differentially expressed orthogroups for each species and (2) the number of background orthogroups, the supertest function returns the probability that the orthogroup sets were independent random samples from the background population of orthogroups. If the null hypothesis is rejected (P-value < 0.05), overlaps between orthogroup sets within a comparison can be considered to be nonrandom recruitment of orthogroup(s) by each lineage. The supertest method assumes that the orthogroup sets being compared are unbiased independent random samples from a population, but the method is robust to violations of this assumption (Wang et al. 2015).

Phylogenetic and Expression-Based Clustering of Species

To determine the assignments of orthologues to orthogroups, OrthoFinder generates an estimate of the species tree for all input species (supplementary fig. S1, Supplementary Material online). This estimate agrees with recent phylogenomic estimates, demonstrating that (1) our assembled transcriptomes allow accurate phylogenetic inference and (2) orthogroup assignments in our study can be considered to be accurate. As a form of comparison, we carried out additional maximum-likelihood analysis of a concatenated alignment of all single-copy orthogroups using IQTREE v2.0.6 (Minh et al. 2020) with a JTT AA substitution model. This analysis recovered an alternative topology for the lizard species within our study (supplementary fig. S3, Supplementary Material online), but this alternative topology was also recovered in a recent phylogenomic study of lizards when a concatenation method was used instead of a multispecies coalescent approach (Brandley et al. 2015). Overall, we can be confident in the orthology assignment procedures.

We tested for convergence in gene expression across lineages by estimating a phylogeny based on gene expression values (fig. 3). Specifically, we tested whether any clustering occurs based on parity mode, nutritional mode, or the stage of gestation. First, we estimated expression at the orthogroup level by summing the expression of all genes within each orthogroup and then applied a variance-stabilization transformation to the data matrix. We then generated a distance matrix based on sample–sample correlations in gene expression. Each species was represented in the analysis by a mean expression value for all gestating replicates and a mean value for all nongestating replicates. We then used the boot.phylo function of the ape R package (Paradis and Schliep 2019) to estimate 10,000 bootstrap replicates of the phylogenetic relationships between species using a neighbor-joining approach, rooted using midpoint rooting. We visualized the resulting trees using the ggtree R package (Yu et al. 2017). We also used principal components analysis to observe whether any clustering occurs based on parity mode, nutritional mode, or stage of gestation (fig. 3).

Supplementary Material

Supplementary data are available at Molecular Biology and Evolution online.

Click here for additional data file.