Transcriptome analysis provides a blueprint of coral egg and sperm functions

Sperm sample collection and sequencing library preparation

M. capitata sperm were collected from the fringing reefs on the west side of Moku o Lo’e and RNA extracted as described in Putnam et al. (2017). Three individual RNA-seq libraries were generated using 200 ng of the total RNA from each sample using the Illumina TruSeq RNA Library Preparation Kit v2. The libraries were combined in equimolar concentrations and run on a single Illumina MiSeq flowcell using the Illumina MiSeq Reagent Kit v3 (150 cycles, single-end). The genomic and transcriptomic data created for this project are available under NCBI BioProject  PRJNA339779 and the three libraries were pooled together and can be found in BioSample: SAMN05607941; Sample name: Mcap_nonBleach_RNA; SRA: SRS1632867.

Transcriptomic and genomic data analysis of sperm

After trimming for quality (parameters shown in Fig. S1) using the CLC Genomics Workbench 8.5.1 (Qiagen, Hilden, Germany), the combined sperm RNA-seq data yielded 29,815,942 high quality reads for assembly. The M. capitata genome (Shumaker et al., 2019) and structural annotations used for the mappings are available at http://cyanophora.rutgers.edu/montipora/. The libraries were individually mapped to the reference genome using CLC Genomics Workbench (count data can be found in Table S1). Only “unique exon read” counts (i.e., the number of reads that match uniquely to exons, including across exon-exon junctions) were used for downstream analyses.

Egg sample collection

In the June 2018 (6/13/18) spawning period, egg-sperm bundles were collected from ambient conditions in the field on the fringing reefs at the HIMB under Special Activity Permit 2018-50 from the Hawai’i Department of Aquatic Resources. Bundles were brought back to the lab and buoyant eggs were separated from the dense sperm after the bundles broke up. The eggs were subsequently rinsed 3 times with 0.2 µm filtered seawater and then following the removal of the water, replicate tubes of eggs were snap frozen in liquid nitrogen and stored at −80 °C.

Egg RNA extraction and sequencing library preparation

RNA was extracted using a dual DNA/RNA extraction method (Quick-DNA/RNA Miniprep Plus Kit, Cat# D7003). A total of 300µl of Zymo DNA/RNA shield, 30 µl of PK digestion buffer was added to each sample tube, followed by the addition of 5µl Proteinase K. Next, samples were vortexed, spun down, and placed in a Thermomixer for 1 h at 55 °C, shaking at 1,100 rpm. After heating, samples were and centrifuged at 2,200 rcf for one minute to separate any remaining solids. The supernatant was transferred to a new microfuge tube and equal volume of DNA/RNA Lysis Buffer was added and mixed. All liquid was transferred to a new Spin-Away filter column with collection tube and spun at 16,000 rcf for 30 s. The flow through containing the RNA fraction was added to a 1.5 mL tube with an equal volume of 100% EtOH was and tubes were vortexed and spun down. A total of 700 µl of the resulting mixture was added to the RNA spin columns where it was centrifuged at 16,000 rcf (g) for 30 s and the flow through (Zymo kit waste) was subsequently discarded. This step was repeated and then 400 µl DNA/ RNA Wash Buffer was gently added to each RNA column. The samples were centrifuged at 16,000 rcf (g) for 30 s and the flow through (Zymo kit waste) was discarded. Next, 80 µl DNase I treatment master mix (75 µl DNA Digestion buffer × # of samples, 5 µl DNase I × # of samples) was added directly to the filter of the RNA columns and incubated at room temp for 15 min. Then, 400 µl of the DNA/ RNA Prep Buffer was gently added to each column and the mixture was centrifuged at 16,000 rcf (g) for 30 s and the flow through (Zymo kit waste) was discarded. The columns containing the bound RNA were then transferred to new 1.5 mL microcentrifuge tubes and 50 µl of warmed DNase/RNase free water was added to each RNA column by dripping slowly directly on the filter. The samples were incubated at room temperature for 5 min and centrifuged at 16,000 rcf (g) for 30 s. This step was repeated to obtain a final elution volume of 100 µl, and tubes were stored at −80 °C. The cDNA libraries were prepared and sequenced by Genewiz using standard Illumina strand-specific RNA-seq preparation with poly-A selection and then sequenced with the HiSeq instrument using 2 × 150 bp reagents with ∼15M raw paired-end reads per sample.

Transcriptomic and genomic data analysis of eggs

After trimming for quality (parameters shown in Fig. S1) using the CLC Genomics Workbench 8.5.1 (Qiagen, Hilden, Germany), the combined egg RNA-seq data yielded 125,054,213 high quality paired reads for assembly. The M. capitata genome (Shumaker et al., 2019) and structural annotations used for the mappings are available at http://cyanophora.rutgers.edu/montipora/. The libraries were individually mapped to the reference genome using CLC Genomics Workbench 8.5.1 (QIAGEN, Aarhus Denmark) (count data can be found in Table S2). Only “unique exon read” counts were used for downstream analyses.

Functional analysis

Using GO terms to organize KEGG pathway data

These sets of “expressed genes” were assigned KEGG Orthology (KO) terms using the KofamScan software (Aramaki et al., 2019), yielding 4,078 genes in the “egg” set with at least one KO term and 3,146 genes in the “sperm” set with at least one KO term. Each set of KOs (egg, sperm) was mapped onto well-studied biological pathways using the “Reconstruct Pathway” tool (Kanehisa, 2017). Next, all KEGG pathways associated with each KO (excluding human-specific pathways) were retrieved and cross-referenced with GO terms of the category “cellular component”. This allowed for assignment of each KO term and its associated pathway(s) to a putative location within the cell. The exhaustive list of all KEGG pathway-GO pairs is included in Tables S4 and S6 and the most common terms relevant to each cellular component were selected as the functional pathways to display in our analysis. This analysis was done to make an overall assessment of which metabolic processes are most frequently expressed in egg and sperm cells based on the normalized read counts for RNA-seq transcript data, and to identify if there are any major differences in overall functionality between each gamete type that can be further investigated.

Differential gene expression analysis

DEG analysis was done in RStudio 3.5.3 (RStudio Team, 2018) using the DESeq2 package (Love, Huber & Anders, 2014) using the raw counts as input (parameters shown in Fig. S2). The count matrices and column data used as input can be found in Tables S8, S9, and S12. Differentially expressed genes were identified in contrasts between egg and adult samples and sperm and adult samples, as well as between egg and sperm samples. Egg replicates corresponding to individual RNA-seq runs are labeled E1, E2, E3, sperm samples are labeled Ub2, Ub3, Ub4, and ambient control treatment adult samples from Shumaker et al. (2019) are labeled W1, W5, and W7 in all related tables and figures.

Functional analysis

Fisher’s exact test revealed similar functional enrichments among the “expressed” gene sets (TPM > 100) in both gamete types. Among the Biological Process and Molecular Function GO terms enriched in the “expressed” genes in the egg samples, the top ten are ‘ATP binding’, ‘GTP binding’, ‘GTPase activity’, ‘structural constituent of ribosome’, ‘chromatin binding’, ‘microtubule motor activity’, ‘negative regulation of transcription by RNA polymerase II’, ‘protein polyubiquitination’, ‘translation initiation factor activity’, and ‘ubiquitin ligase activity’. The top ten enriched GO terms for sperm “expressed” vs. all genes are ‘ATP binding’, ‘GTP binding’, ‘GTPase activity’, ‘structural constituent of ribosome’, ‘chromatin binding’, ‘microtubule binding’, ‘endonuclease activity’, ‘microtubule motor activity’, ‘negative regulation of transcription by RNA polymerase II’, and ‘protein polyubiquitylation’. The complete list of the GO terms for each of these enrichment profiles and their dataset of origin are presented in Table 1. These GO lists, their prevalence in the M. capitata genome, and degree of enrichment are shown in Fig. S3. All GO term enrichments beyond the top 30 shown here were minor and not reported.

Using GO terms to organize KEGG pathway data

After indexing KO terms by their corresponding Cellular Component GO terms, we determined how much of each transcriptome (TPM of all genes assigned to each GO term divided by total TPM of all genes in each “expressed” dataset) is associated with each cellular component and found this metric to be very similar between the egg and sperm datasets (see Fig. 1). To then determine the underlying functional pathways associated with these locations within the cell, the top KEGG pathways for selected parts of the cell were identified, excluding those that are human/ mammal-specific or disease-related. Each KO term was counted for each part of the cell and the most frequent terms were included in Fig. 2. The complete list of all of the raw data used in this analysis can be found in Tables S5 and S7. These results further support the Biological Process and Molecular Function GO data shown above by confirming the similarity between the egg and sperm with respect to which pathways are most prominently expressed and where they are active within the cell. The functional implications of these differences will be discussed below.

Differential gene expression analysis

Gamete vs. adult

From a broader perspective, and perhaps surprisingly, M. capitata egg and sperm appear to have similar patterns of gene expression and share functional enrichment patterns. However, in order to gain a more in-depth understanding of how each of these cells function, how those functions compare with gene expression in adult M. capitata tissue, and how egg and sperm gene expression is related to each other, it is crucial to statistically test for differential gene expression. As expected, given the divergent different tissue types analyzed, principle components analysis (PCA) shows strong differentiation in gene expression patterns between gametic and adult libraries (Fig. S4).

More specifically, with respect to the DEG analysis of egg vs. adult RNA-seq data, DESeq2 identified 13,890 transcripts that were significantly differentially expressed (defined in this study as FDR-adjusted p-value < 0.05), of which 4,487 were up-regulated and 9,403 down-regulated. With respect to the differential expression analysis of sperm vs. adult RNA-seq data, DESeq2 identified 9,717 transcripts that were significantly differentially expressed (defined in this study as FDR-adjusted p-value < 0.05), of which 2,985 were up-regulated and 6,732 down-regulated. Log2-fold change estimates (L2FC), all FDR-adjusted p-values, putative annotations, blastx hits and percent identity (PID), KO terms, and GO terms are provided for all of these genes the Tables S10 and S11. The top ten DEGs that had BLAST hits for egg up- and down-regulated and sperm up- and down-regulated are shown in Tables S3 and S4 with all of their differential expression and BLAST statistics and the full lists of DEGS and accompanying annotations are in Table S3. It is difficult to use these data to gain specific insights because many of the DEGs lack BLAST hits, and those that do, most are to predicted, hypothetical, or uncharacterized proteins in recently sequenced corals. This trend is illustrated in Fig. S5 that provides a snapshot of the degree of similarity between the egg and sperm DEG lists as well as the limitations associated with DEGs that have poor annotation.

To understand some of the basic functions of the egg and sperm transcriptomes, four sets of DEGs were used as the test sets for separate GO-enrichment analyses (Fisher’s Exact Test) against the reference set of all genes in the M. capitata genome (Blast2GO OmicsBox 1.1.164; Götz et al., 2008): (1) egg up-regulated genes (L2FC >1.5; 3,645 DEGs), (2) egg down-regulated genes (L2FC <−1.5; 8,760 DEGs), (3) sperm up-regulated genes (L2FC >1.5; 2,795 genes), and (4) sperm down-regulated genes (L2FC <−1.5; 6,606). The top enriched GO terms for egg up-regulated genes are “ATP binding”, “phosphatase activity”, “positive regulation of biosynthetic process”, “protein serine/threonine kinase activity”, and “GTP binding”. The top enriched GO terms for egg down-regulated genes are “signaling receptor activity”, “regulation of transcription by RNA polymerase II”, “lipid metabolic process”, “organonitrogen compound catabolic process”, and “regulation of signaling”. The top enriched GO terms for sperm up-regulated genes are “purine nucleotide binding”, “purine ribonucleoside triphosphate binding”, “enzyme binding”, “positive regulation of RNA metabolic process”, and “regulation of localization”. The top enriched GO terms for sperm down-regulated genes are “positive regulation of cellular metabolic process”, “regulation of transcription by RNA polymerase II”, “carboxylic acid metabolic process”, “regulation of cellular protein metabolic process”, and “animal organ development”. The full datasets are represented by the bar charts of combined Biological Process and Molecular Function “most specific” enriched GO terms in Fig. S6.

To further resolve the enrichment and DEG data, separate lists of up-regulated (L2FC >1.5) and down-regulated (L2FC <−1.5) genes (same sets as used above for Fisher’s Exact Test) were combined for each of the egg and sperm, indexed with their KO terms retrieved from Kofamscan (Aramaki et al., 2019), and then uploaded into the KEGG “Reconstruct Pathway” tool (Kanehisa, 2017). For every pathway, the up- and down-regulated genes were noted and based on trends seen across egg and sperm data in the above enrichment analysis, connections between KEGG pathway activity and known physiological mechanisms of egg and sperm data were made. One major finding from pathway analysis is that nearly every ribosomal protein in both the egg and sperm datasets is down-regulated (Fig. S7). Another intriguing result, when considering DEGs from the egg vs. adult tissue is that “lipid metabolic process” is among the most down-regulated GO terms. Eggs are rich in lipids and lipid activity (Arai et al., 1993; Figueiredo et al., 2012), therefore our data appeared to be anomalous, or resulting from annotation bias. However, upon further inspection, we find that nearly every down-regulated gene associated with this GO term is annotated as a lipase. Therefore, activities associated with lipid breakdown are those down-regulated in the egg.

Sperm vs. egg comparison

PCA analysis of egg vs. sperm cDNA libraries does not show as marked a difference, as found between the egg and sperm libraries compared to the adult (PC1 of egg vs. adult shows 95% variance, PC1 of sperm vs. adult shows 89% variance, and PC1 of egg vs. sperm shows 47% variance; Fig. S4). Specifically, with respect to the DEG analysis of egg vs. sperm RNA-seq data, these results are more limited in scope. DESeq2 identified 247 transcripts that were significantly differentially expressed (defined in this study as FDR-adjusted p-value < 0.05), of which 108 were up-regulated in the egg (down-regulated in the sperm) and 137 were down-regulated in the egg (up-regulated in the sperm). L2FC was >1.5, <−1.5 for all transcripts, therefore no further filtering was needed. L2FC values, all FDR-adjusted p-values, putative annotations, blastx hits and percent identity (PID), KO terms, and GO terms are provided for all of these genes in Table S13. The sets of egg vs. sperm DEGs were too small to do Fisher’s Exact Test for GO enrichment and too sparse to fill in KEGG pathways, therefore the most useful annotations for these genes are from BLAST. The top ten DEGs that had BLAST hits for this analysis are shown in Table 2 and all annotations are listed in Table S13.

Translation

One of the starkest differences to emerge from the DEG analysis is the marked down-regulation of nearly every ribosomal protein-encoding gene in both the egg and sperm datasets when compared with the adult (Fig. S7). This trend has been documented in human and mammalian systems where sperm are “translationally silent”; i.e., cytoplasmic ribosomal assembly and thus activity is not fully functional when the sperm reaches maturity and nuclear-encoded transcripts may be primarily translated on mitochondrial ribosomes (Gur & Breitbart, 2006; Zhao et al., 2009; Amaral et al., 2014). This trend has also been documented in the non-coral egg literature as well with studies of Xenopus eggs (Smits et al., 2014) and in mice oocytes, where ribosomal protein expression is repressed during late stage oocyte development (Taylor & Pikó, 1992). As highly specialized cells with much lower overall gene expression compared to adult cells, it is not surprising for ribosomal protein genes and thus, translation to be down-regulated in both the sperm and egg. However, this topic needs to be further studied to determine whether this phenomenon in sperm is due solely to its role as a gamete or because the cytoplasmic ribosomes may be translationally inactive.

Fertilization

In broadcast-spawning corals, there are pre-zygotic barriers that gametes must overcome to achieve fertilization (Monteiro, Serrão & Pearson, 2012; Monteiro et al., 2016), including successful chemical signaling between gametes. Chemical signals are secreted by egg cells to attract sperm. Whether this process is utilized to guide sperm to an egg within the body of an individual organism or the broader environment depends on the system being studied. Regardless, there are two major components to this process: secretion of the chemical attractant by the egg, and taxis initiated by chemical receptors, and made physically possible by motile cilia/ the flagellum of the sperm. In mammalian sperm, for example, motile cilia and the sperm flagellum develop in a similar fashion, have the same axoneme structure, and are virtually identical (Clermont, Oko & Hermo, 1993; Avidor-Reiss & Leroux, 2015; Wachten, Jikeli & Kaupp, 2017). The evolution of compartmentalized sperm cilia (via primary cytosolic ciliogenesis) is necessary for the flagellar movement that propels the sperm cell through its environment to the egg. This is thought to be a universal process in metazoans because it is present in both protostomes (Drosophila) and deuterostomes (mammals, humans), as well as in basal metazoans (corals, sponges) (Avidor-Reiss & Leroux, 2015). In broadcast-spawning corals like M. digitata (a relative of M. capitata), eggs play a significant role in the regulation of sperm activity and fertilization success, including sperm signaling and stimulation of flagellar motility (Coll et al., 1994; Morita et al., 2006). The signals produced by the egg must be species-specific to prevent hybridization during mass spawning events where many species release gametes into the water column at the same time. In organisms like ascidians and echinoids, Ca²⁺ has been shown to induce sperm flagellar motility (Yoshida et al., 2003; Morita et al., 2006). A similar effect has been shown experimentally in corals in Acropora species although this alone is insufficient to explain the species-specific nature of sperm and egg union (Morita et al., 2006). It is difficult to find evidence of these complex processes in the summary data of M. capitata discussed above due to the large number of genes lacking annotations. However, upon taking a closer look at the GO data used for Fig. 2, for “motile cilium”, 79 genes are associated with this GO term, 38 are present in the sperm with TPM counts >100, and two of those (g63277 and g9762) are also linked to the GO term “calcium ion binding”. Neither of these genes were assigned a KEGG annotation which is why the data for flagellar and motile cilia cellular components are scarce for this analysis. Both of these genes have the same top blastx hit: XP_015777656.1 PREDICTED: uncharacterized protein LOC107355583 isoform X1 in Acropora digitifera. Both of these genes have transcripts that are significantly up-regulated in the sperm vs. adult data and egg vs. adult data. Based on these results, it is possible these genes play a key role in communication associated with the fertilization process. Differential expression and BLAST annotation statistics for these two genes are presented in Fig. S8. We speculate that these genes are interesting targets for future CRISPR/Cas9 based gene knockdowns (recently developed for Scleractinia; Cleves et al., 2018; Cleves et al., 2020) to explore coral reproductive biology.

Differential expression between egg and sperm

As mentioned throughout this study, it is difficult to gain conclusive insights into the functional capabilities of different M. capitata cell types due to the lack of annotations for many DEGs. Of the 108 genes significantly up-regulated in the egg when compared to sperm, 26 of these DEGs (∼24%) have top blastx hits to uncharacterized proteins and 12 (∼11%) lack an annotation. Of the 139 significantly up-regulated genes in sperm, 32 (∼23%) have top blastx hits to uncharacterized proteins, whereas 13 (∼9%) lack an annotation. These results demonstrate that ca. 1/3 of the genes in this dataset, including many of the most highly differentially expressed genes (based on L2FC), do not provide functional information. It is worth noting however, that many of the annotated DEGs in the egg set encode proteins that are localized to the nucleus and take part in mitosis-related or DNA-based processes, whereas the proteins encoded by sperm DEGs often have functions related to motility and ATP-binding. These limited findings are consistent with general gamete physiology, whereby eggs express cell division, transcription, and DNA repair functions and sperm are optimized for energy production to support motility.