Inflammation and convergent placenta gene co-option contributed to a novel reproductive tissue.

The evolution of pregnancy exposes parental tissues to new, potentially stressful conditions, which can trigger inflammation.1 Inflammation is costly2,3 and can induce embryo rejection, which constrains the evolution of pregnancy.1 In contrast, inflammation can also promote morphological innovation at the maternal-embryonic interface as exemplified by co-option of pro-inflammatory signaling for eutherian embryo implantation.1,4,5 Given its dual function, inflammation could be a key process explaining how innovations such as pregnancy and placentation evolved many times convergently. Pelvic brooding ricefishes evolved a novel "plug" tissue,6,7 which forms inside the female gonoduct after spawning, anchors egg-attaching filaments, and enables pelvic brooders to carry eggs externally until hatching.6,8 Compared to pregnancy, i.e., internal bearing of embryos, external bearing should alleviate constraints on inflammation in the reproductive tract. We thus hypothesized that an ancestral inflammation triggered by the retention of attaching filaments gave rise to pathways orchestrating plug formation. In line with our hypothesis, histological sections of the developing plug revealed signs of gonoduct injuries by egg-attaching filaments in the pelvic brooding ricefish Oryzias eversi. Tissue-specific transcriptomes showed that inflammatory signaling dominates the plug transcriptome and inflammation-induced genes controlling vital processes for plug development such as tissue growth and angiogenesis were overexpressed in the plug. Finally, mammalian placenta genes were enriched in the plug transcriptome, indicating convergent gene co-option for building, attaching, and sustaining a transient tissue in the female reproductive tract. This study highlights the role of gene co-option and suggests that recruiting inflammatory signaling into physiological processes provides a fast-track to evolutionary innovation.

Results

Histology reveals signs of microscopic injuries and indicates presence of immune cells in the plug

Most ricefishes deposit their eggs within hours after spawning. In contrast, female pelvic brooders carry eggs until hatching, (Figure 1). Eggs remain connected to the gonoduct by attaching filaments (AFs), cell-free strings on the egg surface that are built during oocyte maturation. While AFs of egg-depositing species fix eggs to a substrate, AFs of the pelvic brooding Oryzias eversi and Oryzias sarasinorum are retained and eventually anchored in the plug, (Figure 1). To date, tissue origin and cell types in the plug remain largely unknown. To gain insight into the microanatomy of the developing plug and investigate whether retaining AFs injures the gonoduct, we generated histological sections of the pelvic region.

Figure 1

Anatomy of pelvic brooding in Oryzias eversi

(A) A brooding female of O. eversi.

(B) The plug (pl) is connected to the eggs (e) by attaching filaments (af).

(C) A schematic illustration of the pelvic region during pelvic brooding. Eggs are situated next to the genital papilla (gp), covered by elongated pelvic fins (pf) and attached to the mother via attaching filaments (af) that terminate in the plug (pl). The plug is anchored inside the mother’s gonoduct (gd) (o, ovary; sk, skin).

(D) Microanatomy of the developing plug 7 days after spawning (longitudinal section, trichromatic Masson-Goldner staining). The upper left inlay provides an overview of the developing plug and shows the location of the main picture (black square). Interstitial cells (ic) surround the attaching filaments. Blood vessels (bv) enter the posterior part of the plug at a potentially injured site where attaching filaments (af) reach into the gonoduct epithelium (ge) (marked by an asterisk). Lower right inlay shows a multinucleated giant cell (for further details, see Figure S1).

At 1 day after spawning, individual small interstitial cells adhered to attaching filaments throughout the incipient plug (Figures S1A and S1B). Blood capillaries were absent. Seven days after spawning, at least two types of interstitial cells surround AFs in the developing plug. While these cells stained similar to epithelial cells, they differed in shape and size from cells in neighboring tissues including the gonoduct epithelium (Figures 1D, S1, and S2). Similar to observations in O. sarasinorum, the plug was still largely separated from the gonoduct. Notably, we found potentially injured sites where AFs reached into the gonoduct epithelium (Figures 1, S1C, and S1D). At these sites blood capillaries were highly abundant and entered the plug through protrusions of gonoduct tissue that embedded the tips of AFs (Figures 1D and S2B). Cells with similar morphologies to plasma cells, granulocytes, and macrophages forming multinucleated cells were visible in the plug (Figures 1D and S1). Multinucleated giant cells, which are a hallmark of chronic inflammation, were highly abundant in the plug at the end of brooding, but were never observed in gonoduct tissue (Figures 1D and S1Q).

Gene expression is organ-specific and impacted by brooding status

To gain insight into the molecular basis of pelvic brooding, we analyzed a total of 42 tissue-specific transcriptomes. These included the plug tissues of six brooding females, as well as ovaries, genital papillae, and skin samples from six brooding and six non-brooding O. eversi.

From a de novo transcriptome assembly of 1.1 billion paired end reads, we retained a final dataset of 30,643 genes (N50 of 2,917 bp. BUSCO scores of 85.4% complete, 11% missing, 3.6% fragmented, 2.1% duplicated; Tables S1 and S2). Gene expression within organs was more similar than across organs, and differences between brooding states were small compared to variation in gene expression across organs (Figure 2A). The skin and the papilla, a dermal protrusion from which the posterior part of the gonoduct develops, formed sister clusters in hierarchical clustering and principal component analysis (PCA) of gene expression (Figures 2A and S2C–S2F). Similarly, plug and ovary formed sister clusters, albeit with less similar gene expression profiles compared to papilla and skin. More genes were exclusively expressed in the plug (n = 1,498) compared to the papilla (n = 1,108), but less than in ovary (n = 2,012) or skin (n = 2,905) (Figure 2B).

Figure 2

Gene expression across organs and brooding stages in Oryzias eversi

(A) Hierarchically clustered Spearman correlation matrix of gene expression (log2 transformed FPKM). More similar gene expression is indicated with increasing yellow color in the heatmap. Samples from brooding (light blue) and non-brooding (brown) individuals are indicated in the hierarchical clustering tree on the left; organs are indicated in the hierarchical clustering tree on the top (papilla, green; skin, blue; plug, red; ovary, yellow).

(B) The number of genes expressed (FPKM ≥ 2) uniquely in each organ and co-expressed across organs is depicted by a Venn diagram.

(C) Volcano plots show organ-specific differences in gene expression between brooding and non-brooding O. eversi. Log2 fold change (FC) of gene expression on the x axis and the −log10 false discovery rate (FDR) on the y axis illustrate differences in gene expression. Dots represent individual genes with colors according to their FC and FDR (dark gray, FC < 2 and FDR ≥ 0.05; blue, FC ≥ 2 and FDR ≥ 0.05; light gray, FC ≥ 2 and FDR < 0.05). Significantly differentially expressed genes (FC ≥ 2 and FDR < 0.01) are colored in red. Numbers at the bottom indicate the number of significantly overexpressed genes in brooding and non-brooding O. eversi, respectively.

See also Figure S2.

Gene expression was most impacted by the brooding status in the ovary and least impacted in the skin, which is not involved in pelvic brooding (Figure 2C). In line with suppressed ovulation and oocyte maturation, brooding downregulated gene expression in the ovary (275 down- versus 16 upregulated differentially expressed [DE] genes). The opposite was true for the papilla (24 down and 52 up), which appears to be enlarged in pelvic brooders., While downregulated genes in the ovary during brooding are involved in tissue remodeling (mrc2, mmp19, tnfaip6, and angptl2), follicle maturation (fstl1, esyt1, and pcsk9), ovulation (tnfaip6), and follicular steroidogenesis (fstl1 and pcsk9), upregulated genes in the papilla included fish skin mucus components (ZPP, krt8, and C-C motive chemokine) and regulators of mucus secretion (abca12 and noxo1).

The plug transcriptome is dominated by inflammatory signaling

To investigate whether the plug exhibits signs of inflammation, we identified characteristic gene sets for each organ and carried out gene ontology (GO) enrichment analyses with organ-associated genes (overexpressed compared to all other organs with FDR ≤ 10−10, expression fold change [FC] ≥ 2) and organ correlated WGCNA modules (Figure S3; Tables S3 and S4). Plug-associated genes (pl-a) had the highest number of enriched biological processes (BPs) (plug, 595; ovary, 198; skin, 72; papilla, 0), indicative of highly distinct gene functions in the plug. Enriched BPs of the plug module (pl-m) and pl-a genes were largely overlapping (Figure 3B; Tables S3 and S4). Unless stated otherwise, plug genes (PLGs) were part of both the pl-a and the pl-m and enriched BPs were enriched in both gene sets.

Figure 3

Inflammatory immune response dominates gene expression in the plug

(A) Treeplot illustrates enriched biological processes (BPs) in plug-associated genes (pl-a) with area sizes representing significance (large areas: low p values). GO terms were summarized under one representative term with the smallest p value. Redundant terms were removed using REVIGO.

(B) Overlap of enriched BPs in the plug-associated gene set (pl-a) and the brown plug module (pl-m).

(C) Bar plot shows enrichment of BPs linked to inflammation in the pl-a genes (red) and the pl-m (brown). The dotted line indicates FDR = 0.05.

(D) Centered log2 transformed gene expression of myeloid cell marker genes (n = 28). Overexpression is shown in red and under-expression in blue in the heatmap. Similarity in gene expression among samples is shown in the hierarchical clustering tree on the top, while similarity in expression among genes is shown in the hierarchical clustering tree on the left.

(E) Odds ratio of myeloid cell genes is shown for tissue-associated (-a) gene sets and tissue-correlated WGCNA modules (-m). Whiskers show 95% confidence intervals and color indicates whether odds ratios are significantly different from one, i.e., whether focal genes are significantly more or less frequent than would be expected by chance based on the entire transcriptome using Fisher’s exact test.

For further details, see Figures S3 and S4 and Tables S3 and S4.

The immune system dominated enriched BPs in the plug (Figure 3; Tables S3 and S4). “Inflammatory response” and “regulation of acute inflammatory response” were enriched in the plug (Figure 3C). Inflammation protects organisms from infections and restores tissue homeostasis. The inflammatory pathway consists of inducers, molecules that indicate adverse conditions; sensors, receptors of the innate immune system; and mediators that alter functions in effectors, the target tissues.

PLGs included sensors like scarb1 and Toll-like receptors (tlr1 and tlr2), which mediate innate immunity16, 17, 18 and activate NF-κB and mitogen-activated protein kinases (MAPKs). Correspondingly, “positive regulation of MAPK cascade,” which induces secretion of inflammatory cytokines,, and “positive regulation of cytokine secretion” were enriched in the plug. Cytokines are messengers in cellular crosstalk of immune cells that initiate immune cell migration. They are secreted by monocytes and macrophages and signal transduced via the Jak-STAT pathway., “Macrophage activation” and “positive regulation of JAK-STAT cascade” were enriched in the plug. PLGs activate monocytes after trafficking into inflamed tissues (tnfrsf1a, cd40, tlr1, tlr2, and itam) and regulate differentiation of macrophages and monocytes (hck and cxcr3.2)., Proinflammatory signaling, via il1b (pl-m) and the PLGs il12b, il1r1, and il6rb, triggers an “acute inflammatory response” (enriched in the plug). IL-12 promotes the production of pro-inflammatory cytokines and induces PLG il18r expression,, which activates immune cells and induces “positive regulation of IFN-γ production” (enriched in the plug). Further, PLG and enriched BPs are related to activation of diverse immune cells and the complement system (crp, c1qa, c1qb, c1qc, and c1s), which links the innate and adaptive immune system.27, 28, 29 Inflammation ends in a resolution and repair phase that is controlled by macrophages. Correspondingly, PLGs included the anti-inflammatory cd163 and tgfb1, which are expressed by macrophages and mediate tissue repair.,

Myeloid cell genes are overrepresented in the plug

To validate their presence and identify immune cell types in the plug, we analyzed the expression of marker gene sets of T cells (adaptive immunity), myeloid cells, and NK cells (both innate immune response). Myeloid cell genes comprised markers for granulocytes (neutrophils and eosinophils) and monocytes (dendritic cells and macrophages). Myeloid cell marker genes were significantly enriched among pl-a genes (n = 6, odds ratio [OR] = 6.21, p < 0.001; Fisher’s exact test) and the pl-m (n = 11, OR = 3.25, p = 0.005) (Figures 3D and 3E) compared to the rest of the transcriptome. 10,000 randomly sampled genes with one-to-one orthologs confirmed that similar or more extreme results than for myeloid cell genes are expected in less than 0.3% (pl-a) and less than 1% (pl-m) of the cases (Figure S4B). Myeloid cell markers were not enriched in any other organ (Figure 3E). Neither NK cell markers nor T cell markers were significantly enriched in any gene set (Figure S4A). Analyses including all homologs instead of one-to-one orthologs revealed similar results (Figure S4B).

Immune response genes in the plug contribute to tissue remodeling and growth

Regulatory circuits link the immune system and the extracellular matrix (ECM), which controls tissue functions, cell migration, proliferation, and differentiation.33, 34, 35, 36 Immune cells regulate ECM composition and the ECM modulates immune responses.33, 34, 35 In line with collagenous matrix formation, “extracellular matrix organization” and “collagen metabolic process” were enriched in the plug. ECM constituents like collagens (col12a1, col8a1, col5a2, and col5a1) and non-collagenous glycoproteins (fn1a, emil1, lamb2, and lamb4) were PLGs. Emilins, laminins, and collagens are secreted by macrophages after wounding, and “regulation of wound healing” was enriched in the plug. Macrophages also secrete molecules initiating tissue repair including the PLGs platelet-derived growth factor (pdgfc) and tgfb1. The latter regulates expression of ECM molecules and ECM remodeling matrix metalloproteinases (MMPs).,, PLGs also included important matrix remodeling enzymes such as collagenases (like col), decorin (dcn), and several MMPs like ADAMs (mmp2, adam9, adam10, and adam33) and ADAMTS (adamt3) secreted by macrophages following infection or injury. Finally, TNF superfamily members (tnfsf11, tnfsf13, and tnfsf14) and receptors (tnfrsf1a and cd40) controlling cell proliferation and cell death, were PLGs.

Mammalian placenta genes are overrepresented in the plug

A solution for the paradox that novelties evolve although selection can only act on pre-existing traits is co-option: employing existing traits in novel contexts. Accordingly, genes have been recruited non-randomly for certain functions in convergent innovations,40, 41, 42, 43, 44, 45, 46 including pregnancy and placentation in fishes, reptiles, and mammals.47, 48, 49, 50, 51, 52, 53 Similar to mammalian placentation, plug development depends on well-coordinated cell proliferation, ECM formation, and angiogenesis. To explore whether genes with important roles in mammalian placentas were also recruited into the plug, we analyzed expression of 41 one-to-one orthologs of the mammalian core placenta transcriptome. Our results show that these genes are overexpressed in the plug and significantly overrepresented in the pl-m (n = 12, OR = 2.42, p = 0.013; Fisher’s exact test) (Figure 4). 10,000 random samples of 41 one-to-one orthologs between human and O. latipes found in our transcriptome confirmed that similar results occur in less than 3% of the cases (Figure S4C). When analyses were expanded to all homologs (59 genes), enrichment was also significant in pl-a genes (n = 9, OR = 4.11, p < 0.001) (Figure S4C), but not for any other tissue. Thus, genes with important functions in mammalian placentas contribute to plug formation.

Figure 4

Mammalian placenta genes are overrepresented in the plug

(A) Centered log2 transformed gene expression of mammalian core placenta transcriptome genes (n = 41). Overexpression is shown in red while under-expression is shown in blue in the heatmap. Similarity in gene expression among samples is shown in the hierarchical clustering tree on the top, while similarity in expression among genes is shown on the left.

(B) Frequency of mammalian core placenta transcriptome genes compared to the entire transcriptome is shown as odds ratio for tissue-associated (-a) gene sets and tissue-correlated WGCNA modules (-m). Whiskers show 95% confidence intervals and color indicates whether odds ratios are significantly different from one, i.e., whether these genes are significantly more or less abundant than would be expected by chance based on the entire transcriptome using Fisher’s exact test (further details are provided in Figure S4).

Genes involved in pregnancy and placentation expressed in the plug contribute to ECM remodeling, signaling, and angiogenesis

PLGs play central roles in vertebrate pregnancy, and “female pregnancy” was an enriched BP in the pl-m. PLGs that are also part of the mammalian core placenta transcriptome build the ECM (fn1a), modulate ECM-cell interactions (dcn), and contribute to signal transduction (grb2) and tissue remodeling (adam9, serpine2, and p4 ha1b). Serpine2 is upregulated during brood pouch remodeling in seahorses and grb2 controls trophoblast migration during embryo attachment in mammals. Additionally, PLGs included itga4, which is essential in mammalian placenta development, and pcdh10, which is highly expressed in mammalian placentas and likely contributed to the evolution of placentas in poeciliid fishes, seahorses, and pipefishes.

PLGs regulate vascularization in mammalian placentas (tfeb, angptl2, myo1c, and esr1),, and “positive regulation of vasculature development” was enriched in the plug. Estrogen and progesterone levels vary across the reproductive cycle and control processes in reproductive tissues including oogenesis in fishes and endometrial cell proliferation in mammals., Thus, estrogen (esr1) and progesterone receptors (paqr9) may allow the synchronization of processes in the plug with the reproductive cycle.

Additionally, PLGs involved in vertebrate pregnancy are related to the immune system. These include key inflammatory signaling genes that are vital for mammalian embryo implantation (il1r, il1b [pl-m], and il6r) and expressed at the maternal-fetal interface of placental reptiles (il1b and il1r), or during brood pouch development in pipefishes (il6r). Furthermore, tgfb1, MAPK signaling, and pik3cg regulate proliferation and invasion of rat placenta cells. prf1, a cytolytic enzyme, is upregulated in uterine NK cells when these aid embryo implantation in mammals., Moreover, xbp1 is activated by TLR and VEGF signaling, part of the mammalian core placenta transcriptome, and promotes angiogenesis. acod1, a PLG expressed in activated macrophages, contributes to mammalian embryo implantation and downregulates inflammatory responses. Finally, PLGs included the alpha-2-macroglobulin a2ml1. Alpha-2-macroglobulis are acute phase proteins expressed in placentae of live bearing fishes, upregulated in brood pouches during early pregnancy in seahorses, and important for mammalian gestation.,

Discussion

Illuminating the genetic bases of evolutionary innovations is fundamental to understanding how complex life histories and body plans evolve. Here, we analyzed tissue-specific transcriptomes of the pelvic brooding ricefish Oryzias eversi and present evidence that, similar to mammalian embryo implantation, the novel egg-anchoring plug tissue evolved from an acute inflammatory reaction and involved co-option of genes that play important roles in mammalian placentas.

Pelvic brooders rely on a novel transient tissue (the plug) in the female reproductive tract to carry eggs until hatching. Phylogenetic data indicate that ancestors of pelvic brooders deposited fertilized eggs shortly after spawning. Although ecological drivers of pelvic brooding remain unknown, a lack of spawning substrates in pelagic habitats may have favored prolonged egg-carrying., As no “non-self” antigens are required to elicit inflammation, prolonged brooding could trigger inflammation, if retaining attaching filaments injures/irritates the gonoduct. Accordingly, inflammatory signaling dominated the plug transcriptome and blood capillaries entered the developing plug through protrusions of tissue that embedded the tips of AFs (Figures 1D, S1C, and S1D). Our observations indicate that AF can injure the gonoduct, causing inflammation that induces cell migration, cell proliferation, and angiogenesis, i.e., sprouting of new blood vessels from neighboring tissues. Different immune cells, including multinucleated giant cells, were present, and expression of marker genes supports overabundance of myeloid cells in the plug. Multinucleated giant cells are characteristic for granulomas, macrophage-dominated tissues that encapsulate persisting irritants., Myeloid cells further mediate inflammation and subsequent tissue repair, and crucial processes for plug development such as formation of ECM, cell proliferation, and angiogenesis rely on genes that are known to be induced by inflammation. Hence, we propose that an inflammatory response to retention of AFs provided the molecular toolkit for the plug.

This scenario is in line with the recently introduced model of “stress-induced evolutionary innovation” (SIEI). SIEI highlights the innovative potential of recurrent stress responses and proposes that co-opted stress-response pathways can give rise to novel body parts. If the ancestral inflammation already anchored attaching filaments, e.g., by granuloma-like encapsulation of filaments, and thereby produced an approximation of the adaptive phenotype, evolution of the plug would fall somewhere between SIEI and the plasticity-first hypothesis., Since we predicted gene functions based on sequence similarity, demonstrating how inflammatory signaling contributes to plug formation will ultimately need functional validation. While the presence of individual cells on attaching filaments throughout the incipient plug at the beginning of brooding shows that the plug is not built by thickening gonoduct epithelium, the origin of these cells remains unknown. In the future, single-cell transcriptomics will help to characterize the tissue origin and cellular composition of the plug and knockdown of key inflammatory genes will shed light on how inflammation contributes to plug formation. Closely related egg-depositing species like O. dopingdopingensis provide promising opportunities to identify changes in inflammatory signaling as well as genomic regions linked to pelvic brooding.

Our data add to emerging evidence indicating that stress responses, including inflammation, can serve as sources of innovation.,,,,, In contrast to the common infection-centered perspective, inflammation likely evolved as a broad response restoring tissue homeostasis. Although inflammation is not commonly associated with innovation, eutherian embryo implantation is an important exception. Embryo implantation evolved from an acute inflammatory response to embryo attachment, and still relies on pro-inflammatory signaling that induces tissue remodeling and increases vascular permeability., Although inflammation is always costly,, immune cells infiltrating the reproductive tract of ricefishes cannot reach embryonic tissues. Tinkering with inflammation is thus expected to be less constrained compared to pregnancy, where inflammation can cause embryo rejection. Inflammatory signaling pathways and downstream modules control ECM composition, cell migration, cell proliferation, and angiogenesis, making them prime candidates driving the evolution of novel tissues. During the evolution of pregnancy, the involved tissues face new, potentially stressful conditions, so co-option of inflammatory signaling may provide a “jump-start” to complex adaptations. This idea is supported by key pro-inflammatory signaling genes that are expressed during plug formation, eutherian embryo implantation, brood pouch remodeling in pipefishes, and pregnancy in reptiles., Additionally, pro-inflammatory signaling is known to activate PLGs that are also important for pregnancy in both mammals and fishes,,, suggesting that their co-option is more common than previously thought.

In addition to embryo implantation and plug formation inducing overlapping sets of inflammatory genes, our analyses indicate that genes with important roles in placentas were independently co-opted into the plug. In contrast to placentas, which are defined as fusions of fetal and maternal tissues for physiological exchange, the plug is only formed maternally and likely does not serve physiological exchange. Nonetheless, both plug formation in ricefishes and placentation in mammals require well-coordinated cell proliferation, ECM formation, and angiogenesis to build, attach, and sustain a transient tissue in the female reproductive tract. Accordingly, mammalian placenta genes expressed in the plug contribute to ECM formation, cell-cell signaling, and angiogenesis, and some play important roles in independently evolved pregnancy in fishes and mammals.,, Our data thus indicate that these genes were independently co-opted for shared functions in vertebrate placentas and the plug. This finding is in accordance with previous studies that showed independent co-option of certain genetic toolkits into convergently evolved innovations.,,,80, 81, 82

Our study gives first insights into the molecular basis of pelvic brooding, identifies striking parallels to the evolution of pregnancy, and paves the way toward a deeper understanding of how this evolutionary innovation revolutionized parental investment in ricefishes. In addition to highlighting the importance of gene co-option for evolutionary innovation, it provides a new perspective on inflammation, indicating that tinkering with inflammatory signaling provides a fast-track to evolutionary innovation.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Leon Hilgers (leon.hilgers@seckenberg.de).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

Oryzias eversi used in this study are derived from a captive population of ZFMK in Bonn, Germany, established from wild caught individuals from Sulawesi, Indonesia in 2012. Adult fish were kept with 11.5/12.5 h light/dark cycle at temperatures between 25 and 27°C. Tissues were sampled as permitted by the Landesamt für Natur, Umwelt und Verbraucherschutz (§11 Abs. 1Nr 1b, 8a and 8d TierSchG). Each tank was checked for brooding females at least twice per day. To investigate the microanatomy of the developing plug, females were sampled one day (n = 3) and seven days (n = 3) after spawning. For gene expression analyses, twelve brooding females were separated and transferred to smaller net inlets. The ovary, the genital papilla, the skin and the plug were harvested from six egg-carrying females on the 7th day after oviposition (Figure S2B). At this time of sampling, the plug is well developed, but has not yet reached its final size and vascularization (Figure S2A). Therefore, we expect the cells that form the majority of the plug to be present and the molecular machinery that is needed to build the plug still to be active. Organs of non-brooding individuals were harvested from six females three days after the loss of the attaching filaments after a successfully completed brooding cycle. At this time, females of O. eversi (A.S. and T.S., unpublished data) and its sister species O. sarasinorum can spawn again, so we expected that changes in gene expression related to the preparation of the next spawning are already initiated. Brooding and non-brooding individuals were kept under the same conditions. Fishes were deeply anaesthetized with an overdose of MS222 and death was assured by cutting the spinal cord.

Method details

Histology

For histological analyses, specimens were fixed in paraformaldehyde or Bouin-Hollande solution for three to seven days and subsequently decalcified in 10%-ethylenediaminetetraacetic acid for seven days. Abdomens were infiltrated with paraffin, embedded in a paraffin block and cut into 5 μm thick longitudinal sections with a rotary microtome. Sections were mounted on glass slides and stained with a trichromatic Masson–Goldner staining (light green). Photographs were taken using a ZEISS AxioCam “HRC” coupled to a ZEISS Axio Imager.Z2m microscope. All observations reported in this study were verified across multiple slides.

RNA extraction and transcriptome sequencing

Organs were harvested and immediately stored in RNAlater to ensure RNA preservation. After incubation at 8°C for 2 days, samples were transferred to and stored at −20°C. Tissues were homogenized with mortar and pestle (in 1.5ml tube) in 15μl RLT buffer with DTT and total RNA was extracted using a customized protocol of the RNeasy Plus Micro Kit (QIAGEN). Briefly, lysis buffer was diluted and proteinase K was used to digest tissue fragments following mechanical homogenization. Subsequently, additional lysis buffer was added for efficient DNA removal with gDNA spin columns. Amount and quality of total RNA was inspected using Agilent’s 4200 TapeStation. All sequenced samples showed no signs of degradation or DNA contamination. Libraries were prepared using the Illumina TruSeq mRNA stranded sample preparation Kit with 1μg total RNA as input. After poly-A selection, mRNA was purified, fragmented and reverse transcribed using random primers, followed by second strand cDNA synthesis and indexing adaptor ligation. The products were purified and amplified to create the final cDNA libraries. Following library validation and quantification using Agilent’s 4200 TapeStation, equimolar amounts of each library were pooled. The pool was quantified by using the Peqlab KAPA Library Quantification Kit and the Applied Biosystems 7900HT Sequence Detection System. Sequencing was carried out at the Cologne Center for Genomics (CCG) on an Illumina HiSeq4000 and an Illumina NovaSeq6000 sequencer using PE100 and PE75 protocols, respectively (Table S2).

Transcriptome assembly & filtering

Raw reads were trimmed to remove terminal Ns and base calls with a Phred quality score below 30 using sickle v1.33 discarding reads shorter than 25 bp. The resulting 1.1 billion quality filtered reads from all individuals and organs were used for de novo transcriptome assembly (Table S2). Trinity v2.8.4, was run in strand-specific mode with a minimal contig length of 250 bp, in silico read normalization (max. read coverage = 50) and five-fold minimal kmer coverage to generate a single assembly of all organs. De novo assembly resulted in 140,754 contigs that represented isoforms of 66,967 “trinity genes.”

Since ribosomal RNA (rRNA) abundance primarily reflects polyA capture success, rRNA genes were removed following identification with a BLAST search using Oryzias latipes 28S rRNA (AY655694.1) and 18S rRNA (XR002874070) as query sequences. Furthermore, MCSC decontamination was used to remove contaminant transcripts from non-chordates. Finally, the remaining genes and isoforms, as assigned by Trinity, were filtered by expression to exclude spurious transcripts (genes: TMM normalized FPKM ≥ 2, i.e., at least two mapped fragment per kilobase of transcript per million mapped reads, isoforms ≥ 5% of the corresponding gene’s expression), using the filter_low_expr_transcripts.pl script provided with the Trinity v2.8.4 pipeline. Filtering increased the average contig size from 1,674 bp to 1,886 bp and the transcriptome N50 from 2,875 bp to 2,917 bp (Table S1). BUSCO v3.0.2 was used to generate estimates of transcriptome completeness, redundancy and fragmentation by searching for 4,584 Actinopterygii single copy orthologs (odb9).

Transcriptome annotation

The final assembly was annotated with Trinotate v3.0.1. Briefly, the longest open reading frame (ORF) of each transcript was extracted using TransDecoder v5.5.0 (http://transdecoder.github.io). Transcripts and predicted proteins were compared to the UniProtKB/Swiss-Prot database (March, 2019 version) using BLASTP and BLASTX v2.7.1, respectively. Significant BLAST hits (E ≤ 10−3) were extracted and HMMER v3.1b2 was employed to scan transcripts for known protein families based on hidden Markov models stored at Pfam (March, 2019 version). RNAs, signal peptides and transmembrane domains were identified using RNAmmer, SignalP-5.0 and TmHmm 2.0c, respectively. This pipeline provided annotations for transcripts of 19,926 genes, which were imported into the Trinotate-SQLite database, and the annotation report was generated using default parameters. The identity of all transcripts that are mentioned by name in this manuscript were further verified by BLASTN searches against the NCBI reference RNA sequences database (refseq_rna) and manual inspection of alignments with significant matches. It was further manually checked that annotations of organ-associated genes mentioned in this study did not occur in more than one organ-associated gene set.

Gene expression

Mapping and differential gene expression

Bowtie2 v2.3.4.1 was used to map quality-filtered, trimmed reads to the transcriptome and transcript abundance was estimated with RSEM v1.3.1 as implemented in the align_and_estimate_abundance.pl script in Trinity v2.8.4., To alleviate potential biases in mapping accuracy derived from reads with different lengths, reads generated from 100 bp paired-end sequencing were trimmed to 75 bp using FASTX_Trimmer v0.0.14 prior to mapping (http://hannonlab.cshl.edu/fastx_toolkit/). We used removeBatchEffect in limma v3.48.1 to test for potential batch effects caused by sequencing runs, sequencing lanes, and tissue preparation by different experimenters. Inspection of data with and without controlling for potential batch effects did not reveal any batch effects. To avoid overconfidence in gene expression differences between e.g., brooding and non-brooding O. eversi, we continued with the original data. Differential expression analyses were carried out on the gene-level (Trinity genes) and differentially expressed genes (False discovery rate (FDR) ≤ 0.05; Fold change (FC) ≥ 2) were determined for all possible comparisons of organs and brooding states with edgeR v3.28.1 using TMM-normalized transcript abundance estimations. If a gene was highly significantly overexpressed (FDR ≤ 10−10; FC ≥ 2) in one organ compared to all other organs, it was classified as organ-associated. We chose a strict FDR to focus on the most significantly overexpressed genes, which are characteristic of one tissue. The fold change cutoff was not set as stringent, because that would disproportionately remove highly expressed genes. Since the plug is only present during pelvic brooding, organ-associated gene sets were determined based on samples from brooding individuals only.

Weighted gene correlation network analysis

Modules of co-expressed genes were identified using weighted gene correlation network analysis with the R package WGCNA v. 1.68. A gene wise TMM normalized FPKM table was generated within the Trinity v2.8.4 pipeline and log2 transformed to serve as input for WGCNA. We used signed correlation and a minimum cluster size of 20 genes. An ANOVA was used to evaluate correlation of modules with organs and brooding status.

Immune cell marker gene expression

Immune cell marker genes for T cells, myeloid cells and NK cells were chosen based on the top 100 differentially expressed genes of one cell type over the others in the zebrafish Danio rerio provided by Carmona et al. (2017). One-to-one medaka (HdrR, Oryzias latipes) orthologs (Orthology confidence = 1; ranges from 0 (low) to 1(high)) of Danio rerio (GRCz11) immune cell marker genes were retrieved from Ensembl biomart (https://www.ensembl.org/biomart/martview/9cd1b625298fc050fe23948c248b03fa, release 100 - April 2020 © EMBL-EBI) based on their gene stable IDs.

All transcripts in the final transcriptome of O. eversi were searched in the medaka transcriptome (O. latipes (HdrR), ASM223467v1) using BLASTN v2.7.1. An O. eversi gene was chosen as a putative ortholog to a zebrafish (D. rerio) gene, if a transcript of that gene had a significant BLAST hit (E-value ≤ 10−10; nucleotide sequence identity ≥ 85%) that covered at least 50% of a medaka transcript from a medaka/zebrafish one-to-one ortholog.

Oryzias eversi one-to-one orthologs were present for 39 NK, 28 myeloid cell and 40 T cell marker genes in our reference transcriptome. Gene expression data of these genes were extracted from the TMM normalized gene expression matrix and their gene expression across all samples was visualized using the ptR script provided in the Trinity v2.8.4 pipeline. Fisher’s exact test was run in R v3.6.3 to estimate whether marker genes were enriched among organ-associated genes and organ-correlated WGCNA modules, compared to the rest of the transcriptome. In addition to immune cell marker genes, we randomly sampled 10,000 datasets of zebrafish genes with one-to-one orthologs/homologs to medaka and significant blast hits to genes in our dataset. We used these data to estimate how likely it is to observe an outcome that is equal to or more extreme than what we observed for myeloid cell marker genes.

Mammalian core placenta transcriptome gene expression

To explore the expression of genes with an essential role in the mammalian placenta, we retrieved one-to-one medaka orthologs and all homologs (O. latipes, HdrR, Orthology confidence = 1) of 115 human (GRCh38.p13) genes that are part of the mammalian core placenta transcriptome, but no housekeeping genes from Ensembl biomart (release 100 - April 2020 © EMBL-EBI, https://www.ensembl.org/biomart/martview/9cd1b625298fc050fe23948c248b03fa). Oryzias eversi genes were chosen as putative orthologous core placenta genes, if a transcript of these genes had a significant BLAST hit (E-value ≤ 10−10; nucleotide sequence identity ≥ 85%) that covered at least 50% of a medaka transcript from a medaka/human one-to-one ortholog that is part of the mammalian core placenta transcriptome. Gene expression data of putative orthologs of mammalian placenta genes were extracted from the TMM normalized gene expression matrix and visualized using the ptR script of the Trinity v2.8.4 pipeline. To estimate whether mammalian core placenta genes were enriched among organ-associated genes and organ-correlated WGCNA modules compared to the rest of the transcriptome Fisher’s exact test was run in R v3.6.3. In addition to testing mammalian placenta genes, we randomly sampled 10,000 datasets of human genes with one-to-one orthologs/homologs to medaka and significant blast hits to genes in our dataset. This allowed us to estimate a background distribution of similarity and to estimated how likely it is to observe an outcome that is equal to or more extreme than what we observed for mammalian placenta genes.

Gene ontology enrichment

To identify over-represented gene functions in characteristic genes of each organ, gene ontology (GO) enrichment analyses were carried out with organ-associated gene sets and with modules that were highly significantly correlated (p ≥ 0.001) with organs or brooding stages within a certain organ. GO assignments were based on Pfams and BLASTP results of open reading frames against UniProtKB/Swiss-Prot database and assigned gene ontologies were extracted from the Trinotate-SQLite database together with their parental terms. Enriched GO terms were identified in organ-associated gene sets compared to the rest of the transcriptome using GOseq v3.8. For visualization of enriched gene ontologies (FDR ≤ 0.05), GO terms were summarized and redundant terms were removed (allowed similarity: 0.5) with REVIGO.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Biological samples
Oryzias eversi tissue samples	This paper	Table S2; NCBI BioSample SAMN16604494 - SAMN16604534 and SAMN16058878
Chemicals, peptides, and recombinant proteins
0.5% Tricaine mesylate	Sigma-Aldrich, A5040	https://www.sigmaaldrich.com/DE/de/search/tricaine?focus=products&page=1&perPage=30&sort=relevance&term=tricaine&type=product_name
BOUIN-HOLLANDE’s Fixative for IHC	Morphisto	https://www.morphisto.de/en/shop/detail/d/BOUIN-HOLLANDE_Fixierl%C3%B6sung_f%C3%BCr_die_IHC/5/
Histosec Pastilles (without DMSO)	Sigma-Aldrich	https://www.merckmillipore.com/DE/de/product/Histosec-pastilles-without-DMSO,MDA_CHEM-115161
MASSON GOLDNER Trichrom Kit	Morphisto	https://www.morphisto.de/shop/detail/d/F%C3%A4rbekit%3A_MASSON_GOLDNER_Trichrom//9119/
Deposited data
Oryzias eversi RNaseq data	This paper	NCBI BioProject: PRJNA661525, NCBI Sequence Read Archive: SRR13078056 - SRR13078097
Software and algorithms
sickle v1.33	83	https://github.com/najoshi/sickle
Trinity v2.8.4	84,85	https://github.com/trinityrnaseq/trinityrnaseq
BLAST v2.7.1	86	https://blast.ncbi.nlm.nih.gov/Blast.cgi
MCSC	87	https://github.com/Lafond-LapalmeJ/MCSC_Decontamination
BUSCO v3.0.2	88	https://busco.ezlab.org/
TransDecoder v5.5.0	http://transdecoder.github.io	N/A
Trinotate v3.0.1	89	https://github.com/Trinotate/Trinotate.github.io/
HMMER v3.1b2	90	http://hmmer.org/
RNAmmer	91	https://services.healthtech.dtu.dk/service.php?RNAmmer-1.2
SignalP-5.0	92	https://services.healthtech.dtu.dk/service.php?SignalP-5.0
TmHmm 2.0c	93	https://services.healthtech.dtu.dk/service.php?TMHMM-2.0
Bowtie2 v2.3.4.1	94	http://bowtie-bio.sourceforge.net/bowtie2/index.shtml
RSEM v1.3.1	95	https://github.com/deweylab/RSEM
FASTX_Trimmer v0.0.14	http://hannonlab.cshl.edu/fastx_toolkit/	N/A
limma v3.48.1	96	https://bioconductor.org/packages/release/bioc/html/limma.html
edgeR v3.28.1	97	https://bioconductor.org/packages/release/bioc/html/edgeR.html
WGCNA v. 1.68	98	https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/
GOseq v3.8	99	https://bioconductor.org/packages/release/bioc/html/goseq.html
REVIGO	14	http://revigo.irb.hr/
Affinity Designer v1.5.3.69	N/A	https://affinity.serif.com/de/designer/
Zerene Stacker v1.04	N/A	https://zerenesystems.com/cms/stacker