Differential gene expression analysis and cytological evidence reveal a sexual stage of an amoeba with multiparental cellular and nuclear fusion

Sex is a hallmark of eukaryotes but its evolution in microbial eukaryotes is poorly elucidated. Recent genomic studies revealed microbial eukaryotes possess a genetic toolkit necessary for sexual reproduction. However, the mechanism of sexual development in a majority of microbial eukaryotes including amoebozoans is poorly characterized. The major hurdle in studying sex in microbial eukaryotes is a lack of observational evidence, primarily due to its cryptic nature. In this study, we used a tractable fusing amoeba, Cochliopodium, to investigate sexual development using stage-specific Differential Gene Expression (DGE) and cytological analyses. Both DGE and cytological results showed that most of the meiosis and sex-related genes are upregulated in Cochliopodium undergoing fusion in laboratory culture. Relative gene ontology (GO) category representations in unfused and fused cells revealed a functional skew of the fused transcriptome toward DNA metabolism, nucleus and ligases that are suggestive of a commitment to sexual development. However, the GO categories of unfused cells were dominated by metabolic pathways and other processes indicative of a vegetative phase. Our study provides strong evidence that the fused cells represent a sexual stage in Cochliopodium. Our findings have further implications in understanding the evolution and mechanism of inheritance involving multiparents in other eukaryotes with a similar reproductive strategy.


Introduction
Sexual reproduction is common and considered to have originated in the last common ancestors of all eukaryotes; however, its evolution still remains a mystery particularly among microbial eukaryotes [1][2][3][4][5][6]. Sexual reproduction can be defined as a stage in the life cycle involving meiosis-a biological process that reduces the genome complement by one-half (haploid); this is followed by fusion of these haploid cells (gametes), in a process commonly known as fertilization, to form a diploid zygote. This definition is primarily based on observations in macroscopic eukaryotes (e.g. animals and plants) that are usually dimorphic with two distinct sexes. Sexual reproduction in macroscopic eukaryotes is well defined with recognizable cellular and molecular signatures [5,7]. While some variations of sexual reproduction are generally known [1,8], the nature or existence of sex in most microbial eukaryotes is poorly documented [1]. Microbial eukaryotes including amoebozoans display diverse quality of life cycles that involve various types of asexual reproduction and sexual stages (when present) that are usually cryptic [1,9,10]. Recent comparative genomic studies of microbial eukaryotes including amoebozoans demonstrate that these microbes contain genes for sexual reproduction that are actively expressed [11][12][13][14][15]. These discoveries questioned the long-held view that most microbial eukaryotes were strictly asexual. However, sexual development and mechanisms of sex in most of these putative sexual lineages is still not well understood [1,16]. The main hurdle in understanding sexual development in microbial eukaryotes has been a lack of observed sexual activity, probably as a result of their complex and diverse life cycles [9,17]. Sources of difficulty in observing sexual phases include unculturability under normal laboratory conditions, poorly characterized life cycles, and in some cases the occurrence of sex during a dormant cyst stage, when the cell becomes condensed and opaque, precluding observation or manipulation [18][19][20]. We recently described the life cycle of a lens-shaped amoeba (Cochliopodium minus) that undergoes extensive multiple cellular and nuclear fusion during active growth (i.e. the vegetative stage) [16,21], hereafter referred to as Cochliopodium for brevity and since most of the members of the genus display similar behavior [22][23][24]. Cochliopodium possesses genes exclusively involved in meiosis [11,13]. This amoeba was suggested to be a likely useful model organism to circumvent some of the hurdles in studying sexual development in microbial eukaryotes with a similar life cycle [16]. The Cochliopodium has a well-defined and relatively short life cycle (~30 days) and it is easy to grow and maintain in the lab. Unlike most other amoebae (e.g., Entamoeba and Dictyostelium) that are assumed to engage in sex during the dormant, cyst stage, Cochliopodium exhibits sexual-like developmental stages (cellular and nuclear fusion) during the vegetative stage, which makes it an ideal candidate for investigative work involving live experimentation.
The supergroup Amoebozoa includes lineages with diverse life cycles that may or may not involve observable sexual stages. Few lineages of amoebozoans have been confirmed to engage in sexual reproduction [25][26][27][28][29][30][31]. Recent genetic studies demonstrate that all amoebae studied possess or express most of the genes that are recognized as exclusive for sex, albeit scanty physical evidence [11]. The four most common life cycle strategies that are assumed to involve sexual stages in microbial eukaryotes are found in the Amoebozoa. These include those that form sexual cysts [18,19,31], amoeboid or ciliated reproductive cells [32], trophozoite (vegetative) fusing cells [16], and those that alternate between sexual and asexual morphs [33]. In this study, we will elucidate the sexual development of one of the reproductive strategies with a tractable amoeba, Cochliopodium, using cytology and stage-specific differential gene expression (DGE) analysis.
Cochliopodium species grow as single cells with a single nucleus through most of their life cycle. However, in sufficiently dense cultures, they fuse forming larger plasmodial stages, their nuclei migrate within the plasmodium, come into contact and fuse, forming polyploid nuclei [16,22,23]. Subsequently, the merged nuclei undergo division within the plasmodium presumably resulting in a new mix of chromosomes. The plasmodium eventual divides back to uninucleate amoeba cells through plasmotomy (cell fission). The molecular aspects of this behavior have not been studied. Similar cellular fusion behavior is also known among other lineages belonging to the three major subclades of Amoebozoa [1,16]. Cellular fusion and nuclear depolyplodization is also reported in other eukaryotic lineages including some mammalian and cancer cells [34,35]. This cellular behavior is likely ancestral in eukaryotes. Elucidating the molecular aspects of this interesting process in Cochliopodium as a model eukaryotic microorganism will give insights into how microbes with dramatically different life cycles use unorthodox ways to adapt to, and evolve in, changing environments. In this study, we present some of the first evidence, from both cytology and differential gene expression data, that the fused stage in Cochliopodium is sexual. Fused amoebae were observed to express meiosis genes and other cellular processes that are consistent with sexual development. Our findings also have implications in understanding the evolution and mechanism of inheritance involving multiple parents in Cochliopodium and other eukaryotes that use a similar reproductive strategy.

Sample collection, single-cell transcriptome sequencing and assembly
Cultures of Cochliopodium minus (syn. C. pentatrifurcatum [35]) were grown in plastic Petri dishes with bottled natural spring water (Deer Park 1 ; Nestlé Corp. Glendale, CA) with added autoclaved grains of rice as an organic nutrient source to support bacterial growth as prey for the amoebozoans. Cultures were maintained until they reached maximum-growth density that we have consistently achieved in our laboratory cultures. Samples of unfused amoeba cells were collected during the first week of culturing before cellular fusion was at its peak (see [16]). Fused cell samples were collected during the steady state of cellular fusion and fission (8-12 days). Fused and unfused samples were primarily distinguished by their cell size, since fusion status based on ploidy or number of nuclei is difficult to achieve when working with live cells of Cochliopodium. Although this is not an ideal approach, based on our previous experience this method gives a rough starting point to classify cells as fused/unfused [16]. Individual cells (unfused or fused) were picked using a platinum wire loop (tip) or mouth pipetting techniques. Cells collected (1-10) using platinum a wire loop (tip) were directly transferred into 0.2-mL PCR tubes and processed for sequencing as described below. Cells collected using mouth pipetting were transferred to a clean glass slide or sterile agar medium to wash off freely floating bacteria, carried over from culture, with sterile water. This washing step reduces sequence contamination from bacteria during transcriptome sequencing [36]. Cleaned individual amoeba cells (1-10) were transferred into 0.2-mL PCR tubes and processed for sequencing using Seq 1 v4 Ultra 1 Low Input RNA Kit (Takara Bio USA, Mountain View, CA) as described in Wood et al. [13]. We also isolated total RNA from 100 cells, from both unfused and fused cell conditions, using NucleoSpin 1 RNA kit (Macherey-Nagel, Düren, Germany) according to the manufacturer's protocol. The total RNA was processed for sequencing using the SMART-Seq 1 kit as above. For sequencing, libraries were prepared from 1 ng of cDNA using the Nextera 1 XT DNA Library Preparation Kit (Illumina Inc., San Diego, CA USA) according to the manufacturer's instructions. Libraries were quantified on the Qubit 1 Fluorometer using the DNA HS Assay. All libraries were sequenced at Yale Center for Genomic Analysis on a HiSeq 2500 in paired-end, high-output mode with 75 bp reads.
FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to inspect reads of each sample for quality and length. Illumina adaptor sequences and low-quality reads with a score below 28 were removed using BBDuk (Joint Genome Institute, U.S. Department of Energy, Walnut Creek, CA USA). The trimming of low quality reads from both ends ("rl" trim mode) is based on the Phred algorithm implemented in BBDuk. Using the same program, we also removed reads shorter than 35 bp after trimming. The remaining reads were assembled de novo using rnaSPAdes-version 0.1.1 [37] with default parameters for gene inventory and DGE purposes.

Differential gene expression (DGE) analysis
To quantify the abundances of transcripts in our samples, a rapid and accurate program Kallisto [38] was used. Kallisto uses a pseudoalignment procedure that is robust to error detection. Since there is currently no reference genome data available for species of Cochliopodium, a comprehensive de-novo transcriptome was assembled from multiple samples of Cochliopodium minus [13,36,39,40]. This master transcriptome was used for mapping the sequencing reads to generate the counts for each sample. The quantification data was then imported using tximport package in R (programming language) for downstream DGE analysis. Two DGE tools, DESeq2 [41] and EdgeR [42] that were shown to perform well for smaller number of biological replicates, were initially applied [43]. Based on the results of this initial analysis, we designed two additional experimental conditions by selecting samples that showed good clustering and removing bad quality samples (see Results). All analyses were performed in DESeq2 packages in R. We used DESeq2 because it has a better moderation in log fold change (FC) for low count transcripts [43]. Fragments per kilobase of transcripts per million base pairs sequenced (FPKM) was used to estimate the level of gene expression. The Benjamini-Hochberg algorithm was used to adjust the resulting p-values for multiple testing with false discovery rate (FDR<0.05) [44]. Genes with an adjusted p-value (padj) of < 0.05 was determined as statistically significant. Up-and down-regulations of the significantly differentially expressed genes (DEGs) were determined using a threshold of 1.5 for log2 FC. Gene clustering of all the DEGs was plotted using function pheatmap in R.

Gene inventory and functional annotation with Blast2GO
To check the expression of potential sexual related genes in different cell conditions, the BLAST program was used to search our transcripts with a set of sexual related ortholog groups (OG) obtained from OrthoMCL database (http://www.orthomcl.org/orthomcl/) with E-value cut-off of 10 −15 . The best hits were chosen and checked manually for paralogs using phylogenetic tree. We also used HMMER [45] based homology search as a further confirmation of the genes discovered in our BLAST results. A total of 94 genes were investigated, which included 15 meiosis specific genes, 29 fusion and karyogamy related genes and 50 general sex-related genes based on the previous studies of amoebae [11][12][13]. The expression of the genes detected in our transcriptome were then further investigated in terms of their fold change and adjusted p-value among different cell conditions using functions in DESeq2 R packages.
To understand the biological activities among different cell conditions, functional annotation of DEGs identified between fused and unfused samples were performed using Blast2GO [46]. Blast2GO is a bioinformatics tool for high-quality protein functional prediction and provides Gene Ontology (GO) annotation visualization and statistical networks for genetics research. The DEGs that were upregulated in fused and unfused samples were analyzed using Blast2GO for functional annotation. We first used CloudBlast and InterProScan to search homologs to the DEGs against non-redundant databases from NCBI. GO terms were then retrieved by mapping the hits to the functional information stored in the GO databases. Finally, the obtained GO terms were assigned to the query genes for annotation. InterProScan GOs obtained through motifs/domains were also merged with the results of Blast to the annotations. Enzyme Codes (EC) were also provided when EC number or equivalent GOs are available. Statistical information of the function analysis was checked using the 'charts' and 'graphs' functions in Blast2GO. Enrichment analysis was performed to check GO terms that were abundant with DEGs upregulated in unfused or fused cells as a test set, and all of the DEGs as a reference set.

RNA in situ hybridization (RNA-ISH) of selected sex genes
We used a highly streamlined RNA in situ hybridization, RNAscope (Advanced Cell Diagnostics), method to detect six representative sexual related genes in fixed amoeba cells in all stages. These target genes include meiosis specific (Spo11, Mer3, Dmc1, Pch2), karyogamy (Kem1) and plasmogamy (Bni1). RNAscope involves a series of steps including pretreatments, hybridization, signal amplifications, and detection of target genes. The RNAscope 1 Multiplex Fluorescent Reagent Kit v2 (Cat. No. 323100; Advanced Cell Diagnostics) was performed using target probes specifically designed for Cochliopodium minus except one gene (Spo11), which was designed using a target sequence from Acanthamoeba castellanii. The probes were designed according to the manufacturer's instructions and include: 20ZZ probes (Cpe-BNI1cust targeting 221-1315, Cpe-DMC1-cust targeting 2-987 and Cpe-KEM1-cust targeting 2-1386), 17ZZ probe Cpe-MER3-cust targeting 1801-2725, 18ZZ probe Cpe-PCH2-cust targeting 2-965 and 13ZZ probe Cpe-SPO11-cust 2-676 targeting 2-676 (A. castellanii). Actin (15ZZ probe Cpe-Actin-cust targeting 104-1076) was used as a positive control. Adherent cells representing all stages (fused, fusing, and unfused) of the amoeba were transferred into a twowell glass chamber slide (Thermo ScientificTM, Nunc Lab-Tek, Rochester, NY) and fixed using either −80˚C cooled methanol or 4% formaldehyde (Ladd Research, Williston, VT) for 3-10 min. Fixed cells were dehydrated in a series of ethanol concentrations and then reverse rehydrated before treatment with hydrogen peroxide and Protease for permeabilization as stated in the manufacturer's manual. Hybridization of the probes to the RNA targets was performed by incubation in a HybEZ Oven for 2 hours at 40˚C. After two washes, the slides were processed for standard signal amplification steps for single or multiple channels. Chromogenic detection of probes was performed using OpalTM 520 (FP1487001KT) and OpalTM 570 (FP1488001KT) dyes (Akoya Biosciences) and Hoechst 33358 (1:1500) for DNA. The preparations were examined with a Zeiss LSM 700 Inverted Confocal Microscope (Carl Zeiss MicroImaging, Gottingen, Germany).

Single cells transcriptome data
A total of 20 samples designated as fused (9 replicates) and unfused (11 replicates) were collected (S1 Table). The yield and quality of sequence data differed among our single cells transcriptomes (S1 Table). Generally fused or large cells yielded more sequence data than smaller cells. Since transcriptome sequencing from an individual (single) cell did not generate enough data, likely due to the small amount of starting RNA material, we increased the number of cells to 5-10 and collected transcriptome data from these cells using the same approach. We also collected transcriptome data obtained using a regular RNA extraction from 100 cells for both unfused and fused samples (S1 Table). The designations of unfused and fused state of amoeba cells were primarily made on the basis of amoeba size. However, this was used as a guide rather than a predictor of an amoeba stage. Cochliopodium cells show a range of sizes during their life cycle (see [16]). The nucleus is not always easy to visualize due to opaqueness of the cytoplasm and cell inclusions that obscure nuclear observation in live cells. These problems make it difficult in identifying unfused or fused cells based on size and number of nuclei. Given these challenges, we made an effort to select cells based on our observation in previous publications [21-23] as a general guide to discriminate between unfused (~20-30 μm) and fused (> 50 μm) amoeba cells.

Differential gene expression (DGE) analysis
PCA (principal component analysis) plot functions in DESeq2 were used to visualize the overall effect of experimental and batch effects and check for possible sample outliers. Two samples that failed the quality check were removed and a plot of the remaining 18 samples showed good clustering, on the whole. However, some samples either formed a separate cluster of their own or were mixed with another group (S1 Fig). Taking into account the difficulties in distinguishing cell conditions during our sampling, we applied a strict clustering pattern to select four of the unfused samples (YT14-YT17) and four of the fused samples (YT26-YT29) that showed good clustering for two experimental conditions (see S1 Table). This two-sample condition analysis, each comprised of 4 replicates, showed good clustering with PC1 representing 76% of variance (Fig 1). Further inspection of the clustering of the 18 samples showed a possibility of an intermediate group consisting of 3 samples (see below).

MA-plot (log-intensity ratios [M-values] versus log-intensity averages [A-values])
showed that the majority of the transcripts between two-conditions have similar expression values (Fig 2). In total, 881 differentially expressed genes (DEGs) were identified between unfused and fused conditions in our DESeq2 analysis. Among them, 521 DEGs were tested as upregulated in the unfused condition samples and 360 DEGs tested as upregulated in the fused condition samples. The clustered heatmap of all DEGs that were identified between the unfused and fused conditions was generated using transformed normalized counts from DESeq2 R package and showed a correlative pattern of upregulation (downregulation) between the two conditions across all samples analyzed (Fig 3). A closer examination of the expression data structure from our initial analyses revealed a potentially third data cluster (YT24, YT25, and YT30) between unfused and fused sample conditions (S1 Table,

Expression of meiosis and sexual related genes
We examined the expression of the genes previously determined to be involved in sexual development of Cochliopodium and other microbial eukaryotes among different sampling conditions [11,[13][14][15]. The gene inventory conducted in all of our samples showed that 68% (64) of the 94 genes inventoried were present (S2 Table). Gene clustering of the 64 detected DEGs was performed using a transformed, normalized counts from variance stabilizing transformation generated in the DESeq2 package (Fig 4). The genes were categorized based on the role they play in sexual development. These include meiosis, plasmogamy, karyogamy and other sexual and developmental related processes (Fig 4, S2 Table). Results of the DGE analysis for all these genes were shown in supplementary S2 Table with the log2 FC and adjusted p-value for each detected gene between unfused and fused conditions. We also included the results from comparisons that included an intermediate stage to explore our data in depth.
Among all the sexual related genes available for expression examination, 33 of the 64 genes analyzed were upregulated in fused cells, while unfused cells had 30 of these genes upregulated (see S2 Table for adjusted p-value for each gene). One of the meiosis specific genes, Pch2, had no result due to its low expression (S2 Table). Among the 33 upregulated genes in fused samples, six were significant with threshold of adjusted p-value of 0.05 (S2 Table). These include genes that are involved in meiosis (Mer3 and Zip1-like), nuclear congression (Kem1), DNA damage sensing (Rad17) and sister chromatid cohesion (Smc1/3) (Table 1, S2 Table, Fig 4). In unfused samples, five of the 30 upregulated genes were significant (S2 Table), which includes cellular congression (Bni1), plasmogamy (Myo2), double-strand break repair (Lig4/Dnl1) and recombinational repair (Smc5, Msh2, Msh6) (Fig 4, S2 Table). When looking into the expression for each functional category, all four genes involved in plasmogamy and some genes involved in double-strand break repair had upregulation in the unfused condition samples; albeit, most of them didn't pass the significance threshold (Fig 4, S2 Table). Upregulation of genes in recombinational repair and sister chromatid cohesion categories were more prominent in fused samples, most of them with padj >0.05 (S2 Table). The expression pattern of some detected meiosis genes (e.g. Pch2) and other genes in S2 Table could not be analyzed due to low expression among samples or their lack in the reference transcriptome (e.g. Spo11). In general, we saw an overall trend that differentially expressed (upregulated) genes in the fused condition play an important role in sexual development. However, some meiosis specific genes such as Msh4/5 and Dmc1 were observed to be upregulated, without statistical significance, in some of the unfused cell samples (S2 Table).

Functional annotation of DEGs
Blast2GO was applied to perform functional annotation of the two sets of DEGs identified among the unfused (521 genes) and fused (360 genes) samples. Three GO categories (Molecular function, Biological process, Cellular components) were inferred by mapping the sequence information to existing annotation sources. Based on these analyses, we were able to identify key functional categories that can be used to inform about sexual development in Cochliopodium (S7-S9 Figs). GO category representations in fused samples were dominated by nucleic acid (DNA) metabolism, nucleus and associated nuclear components (chromosome, nuclear pore, nucleolus), RNA processing, ATP binding and helicase activities. These are based on the most frequent GO terms and sequence distribution for each GO category (S7-S9 Figs). In addition, the fused samples consisted of some notable functional categories directly or indirectly related to sexual development (cell division) including ligase, histone acetylation, peroxisome fission (PHD-finger protein) and MCM complex (S7-S9 Figs). Enzyme Class (EC) distribution comparison of the two sets of DEGs showed that all the enzyme classes had more sequences in unfused samples than in fused samples, except for the EC class Ligases (Fig 5). Various ligase enzymes were identified in fused samples including E3 ubiquitin-protein ligase previously shown to play a role in sexual development [47]. GO category representations in unfused samples were mostly dominated with metabolic activities (carbohydrate, lipid and protein), oxidation-reduction processes, vesicle transport, cytoskeleton and cellular and intra-

Cytological detection of selected sexual DEGs
We performed RNA in situ hybridization (RNA-ISH) assay of selected genes representative of sexual development in the life cycle of Cochliopodium. RNA-ISH analysis gives both qualitative expression levels and spatial distribution of DEGs within intact cells. Our RNA-ISH results are mostly congruent to the DGE analysis of the selected genes. The meiosis genes, Mer3, Pch2 and Dmc1, and karyogamy gene (Kem1) were expressed mostly, but not exclusively, in fused cells (Figs 6-8A and 8B). The expression of Mer3 was the most prominent in fused cells at different stages; that is, pre-karyogamy (Fig 6A), during karyogamy ( Fig 6B) and in polyploidy nuclear, post-karyogamy (Fig 6C). Mer3 was also expressed in unfused cells but at a much lower level (S10 Fig). The detection of Pch2 and Dmc1 was not as prominent as Mer3 (Fig 7). Pch2 was mostly detected in fused cells with multiple and polyploid nuclei (Fig 7C and 7D). Similarly, Dmc1 was detected mostly in fused samples but some unfused condition samples also expressed it (Fig 7A and 7B). The expression of Kem1 was noticeably visible around the nucleus periphery (Fig 8A and 8B), while the detection of Bni1 was not that prominent; it was expressed both in fused and unfused cells, but mostly in unfused cells (Fig 8C and 8D). We also attempted a colocalization (co-expression) experiment, but this experiment had many technical challenges. We were able to show that Mer3 and Dmc1 were co-expressed in fused cells (S11 Fig).

Fused cells are a sexual stage in Cochliopodium
The present study demonstrates that the fused cells in Cochliopodium both express and upregulate most of the meiosis and sexual related genes. Corroborative evidence from the cytological (RNA-ISH) study also provides spatial and temporal expression for some of these genes during the life cycle of Cochliopodium (Figs 6-8). Evidence for sex in a majority of amoebozoans was based on gene inventory of genomic or transcriptomic data, without supporting physical or stage specific DGE analysis [11]. Limited ultrastructural evidence for sex is available for a few non-model amoebae [25,26,31]. Most of the detailed genetic, DGE and cytological studies are limited to the model organism, Dictyostelium discoideum (e.g. [48][49][50]) and the pathogenic lineage, Entamoeba (e.g. [17,20,[51][52][53]). However, both of these amoebae have dramatically different life cycles and could not be representative of the diverse life cycle observed in the supergroup: Amoebozoa. Our study is the first to report comprehensive DGE and cytological studies of sexual development in amoebae characterized by multiple cellular and nuclear fusions. The findings of this study have important implications in understanding molecular mechanisms of cryptic sexual processes in microbial eukaryotes that share similar behavior. The study of sexual processes in microbial eukaryotes has been curtailed mainly due to lack of observational evidence. This study demonstrates that Cochliopodium can serve as a model to understand some aspect of sexual processes shared among microbial eukaryotes due to its well-defined life cycle, tractability and suitability for genetic manipulation studies [16].

Importance of cytological data in understanding DGE analysis in Cochliopodium
Our DGE analysis is largely corroborated by our RNA-ISH data for the selected sex genes. Even though most of our results from DGE and RNA-ISH were consistent with the expected sexual development of the amoeba, we observed some variations (albeit non-significant) in the DGE analysis among our samples (S2 Table). Among the variations noteworthy of mention include the detection and upregulation of meiosis specific genes in some unfused samples. PLOS ONE some unfused cells expressing sex genes. These unfused cells are likely preparing to enter the sexual stage and can be considered intermediates. This observation is supported by the detection of meiosis specific genes in our RNA-ISH staining in some unfused uninucleate cells, although these genes were consistently detected in higher numbers in fused samples (Fig 6B,  S10 Fig). All of our data sources are based on transcriptome (RNAseq), the detection of these genes in unfused cells, possibly entering sexual stage, might also be attributed to the time lag between transcription and translation. The distinction between unfused and fused cells would have been more apparent if we used proteomic data. Immunostaining of meiosis protein could give substantial informative results pertaining to the actual stage of the amoeba. Our attempt to immunostain proteins encoding SPO11 and DMC1 using commercially available antibodies was not successful. Sex genes are among the fastest evolving genes in eukaryotes [11,14,54]. Some sex genes that serve a similar function have evolved far beyond recognition making conventional sequence homology comparisons impractical. These genes can only be identified through a structural domain comparison or their localization in cytological studies (see [54]). Sex genes show great divergence even among closely related species. For example, the RNA-ISH probes designed for Cochliopodium did not yield good results in closely related species of amoebae. Proteomic research of sex genes is critical for determining sexual stages by providing physical (structural) localization of sexual stages as it occurs (e.g., synaptosomal complexes). However, this will require development of species-specific antibodies, which is challenging due to lack of genome data and the associated cost and time required to develop species-specific antibodies. With genome sequencing of many amoebae underway, a more comprehensive program of proteomic and gene manipulation experiments will be feasible in the near future, which will play an integral role in unraveling the various molecular mechanism of sexual strategies employed in the supergroup.
Another variation in DGE analysis relates to upregulation of some meiosis genes in unfused samples (Table 1, S1 Table, Fig 4). Although the upregulations of these genes (Msh5, Msh4 and Dmc1) were not significant (S2 Table), their detection in unfused samples requires an explanation. As described above, our data for unfused samples come from 5-10 cells; and hence it is likely that our unfused samples might have included a mix of cells among which some cells were preparing to enter a sexual stage. Our RNA-ISH data showed that unfused samples to some extent were observed to express meiosis genes (Fig 6B, S10 Fig). These unfused samples cannot be distinguished by cell size alone in live cultures and might partly explain the discrepancies between our DGE and RNA-ISH results. The detection of lowly expressed genes (e.g., Pch2) and a meiosis gene (Dmc1), shown to be upregulated in DGE analysis in unfused samples, were seen in our RNA-ISH to be more abundant in fused compared to unfused samples (Figs 6 and 7, S10 Fig). Our RNA-ISH results are consistent with the expected life cycle of the amoeba. Another possible explanation for the observed discrepancies between DGE and RNA-ISH could be explained by the high adjusted p-value (18.9%) of Dmc1, which can be interpreted as a false positive upregulation [44]. The cytological data greatly helped us to interpret DGE results and overcome some of the practical limitations of cell stage identification in Cochliopodium. Below we highlight the role of key genes, differentially expressed (detected) in our samples, in the sexual development of Cochliopodium.

DGE and cytological evidence point to the main crossover pathway employed in Cochliopodium
One of the most significantly upregulated DEGs in fused samples is Mer3 (Table 1, Fig 4). This gene is a meiosis specific gene and a member of the ZMM group described in budding yeast [55,56]. Mer3, together with other members ZMM genes (Zip1-4, Msh4 and Msh5), plays a crucial role in sexual organisms by providing a link between recombination and the synaptonemal complex (SC) assembly during meiosis [57,58]. Mer3 is among the highly expressed genes and was easily detected in cytological RNA-ISH staining in higher quantities (Fig 6). It was expressed at various phases of the fused cells including prior to, during and post-nuclear fusion (karyogamy, Fig 6), which is indicative of its critical role in sexual development of Cochliopodium. Mer3, a highly evolutionary conserved DNA helicase, is involved in ZMM genes dependent interference cross-over (class I cross-over [CO] pathway), where double-Holliday junctions are preferentially resolved toward CO. Interference in this pathway prevents nearby crossovers from occurring, thereby ensuring widely spaced crossovers along chromosomes [59,60]. Most of the genes reported to play role in the class I CO pathway were either differentially expressed or detected in fused samples. These include meiosis specific ZMM genes (Zip4, Msh4 and Msh5) and other genes not specific to meiosis, including Mlh1-3, Exo1 and Sgs1 (Fig 4, S1 Table). We also found a significant DEG, Zip1-like, gene; which is one of the most important ZMM group members that forms the central component of SC [61]. However, the reference sequences of the Zip1 gene retrieved from OrthoMCL database were few and quite divergent, which created a challenge for accurate homology assessment of this gene in our samples. Our homology search result for this gene in HMMER passed the inclusion threshold (HMMER c-evalue e-6) with a matching length of 414 bp, but failed the Blast (evalue e-7) threshold. HMMER is designed to detect remote homologs and takes into account shared protein domains in its search algorithm [45]. Based on this result, we decided to report the expression of Zip1-like genes from our transcriptome, which showed significant upregulation in fused cells (Table 1). However, this requires further investigation by searching similar genes in related species and acquisition of the full length of the gene in the genome of Cochliopodium, currently unavailable. Most of the ZMM group genes were described from budding yeast; similar genes playing the same role are known in plants and animals [62,63]. It should be noted that the Zip1 gene of the budding yeast and other eukaryotes shares little sequence homology; however, their proteins share biochemical/structural similarities and localize in SC [54]. Given these observed divergences and the significant upregulation of the Zip1-like gene in Cochliopodium (Table 1), further analysis is needed to determine the identity and role of this gene in sexual development of the amoeba.
The expression and detection of Mer3 and most of ZMM group genes suggests that Cochliopodium predominantly employs class I CO and synaptonemal complexes for meiotic recombination. We also detected Pch2, a protein reported to play a regulatory role associated with other meiosis specific genes such as Zip1 and Hop1 (not detected here) in the maintenance of synaptonemal complex organization [64,65]. Pch2 is one of the lowly expressed genes in our DGE analysis but was shown in our RNA-ISH data to localize more in fused samples (Fig 7C  and 7D). Class I CO is a common pathway of eukaryotes [54,58]; however, the exact mechanism in which class I CO pathway and SC operate in an amoeba characterized by multiparental nuclear fusion awaits further investigation. We also detected Mus81, one of the genes involved in the non-interference class II CO pathway (S2 Table). Even though Mus81 is detected at very low levels in Cochliopodium samples, its presence indicates that class II CO might be used as an alternative crossover pathway in this amoeba.

Significance of cohesin complex in multiparental nuclear fusion
Most of the genes comprising the cohesin complex genes (Smc1-4 and Rad21 as well as the meiotic specific cohesin subunit, Rec8) that are essential in sister chromatid cohesion and their faithful separation, were expressed in fused samples [66][67][68]. Understanding the functionality and molecular mechanism of cohesin proteins is particularly of interest for further investigation in Cochliopodium. Cellular fusion in Cochliopodium is followed by mutiparental nuclear fusion and likely involves mixing of chromosomes from multiple individuals [16]. In conventional mitosis and meiosis, cohesins hold sister chromatids together before anaphase during cell division to ensure equal distribution of chromatids to dividing daughter/gamete cells. Any failure in cohesins results in aneuploidy, which has drastic health and viability consequences [69]. Cochliopodium presents an atypical parasexual process where fusion can result in over 30 nuclei in one large fused plasmodium. Investigation of how Cochliopodium sorts such a large number of multiparental chromosomes precociously during subsequent nuclear and cellular fissions of fused cells will likely unravel a highly evolved or novel mechanism for the roles of cohesins.

Other meiosis and sex related genes
We detected several meiosis and sex-related genes that play a role in mismatch repair, initiation of double-strand break (DSB) and their repair as well as homologous recombination (Fig  4, S2 Table). Although these genes are detected in our DGE analysis, their patterns of expression were neither significant nor consistent in our samples (S2 Table). We report them to provide potentially interesting trends that need to be explored further. For example, genes involved in homologous recombination (HR), Rad51 and Dmc1 (meiosis specific), were expressed and upregulated in our unfused samples. Although Dmc1 was shown to localize with Mer3 in our fused samples (S11 Fig), the temporal expression of these genes needs further investigation. Both of these genes code for recombinases that play a critical role in DNA lesion repair, including double-strand breaks (DSBs), single-strand DNA gaps and interstrand crosslinks during meiosis [70]. Interestingly, genes involved in the alternative DSBs repair mechanism, nonhomologous DNA end joining (NHEJ), were all upregulated in unfused samples. Proteins involved in the NHEJ pathway, including KU70 and KU80, are operational during the vegetative stage (G0/G1 phases of the cell cycle) when sister chromatids are not available [71]. These genes are known to be downregulated during meiosis [72], which is consistent with the DGE results of our fused samples (Fig 4, S2 Table). Lastly, one of the key meiosis genes, Spo11, that initiates meiosis by programed DSB was not detected in our DGE analysis. Spo11 is detected in several species of amoebae, but some amoebae (e.g., Dictyostelium) lack it in their genome [11,12,73]. Spo11 is likely one of the lowly expressed genes, since we only detected it one time in one sample; and our attempt to detect it using RNA-ISH with a probe designed from a closely related species rendered no results.

Other DEGs and gene ontology (GO) provide information on the developmental and sexual stage of Cochliopodium
Additional DEGs and GO categories, mirroring cellular physiology, shed light on the life cycle of Cochliopodium. Genes engaged in plasmogamy and karyogamy were upregulated in unfused and fused samples, respectively (Fig 4, S2 Table). This result is consistent with the observed life cycle of Cochliopodium, where unfused cells entering a sexual stage would be expected to produce proteins enabling them to fuse, and expression of karyogamy genes in fused cells supports the observed nuclear fusion [16]. Among the most dominant DEGs in the GO analyses in unfused samples were cellular processes related to metabolic activity and vesicle-mediated transportation (S7-S9 Figs). These are indicative of an active vegetative phase of our unfused samples. Particularly, detection of a large number of proteolytic enzymes indicates that the cells might be engaging in high protein metabolism (S7 Fig). The most dominate GO terms in fused samples included DNA metabolism, nucleus, ATP binding, PHD-Finger protein and ligases (S7-S9 Figs). These results can be interpreted as fused cells primarily engaged in ATPmediated, DNA-related processes, likely reflecting their sexual nature. Among these, E3 ubiquitin ligase [47] and PHD-Finger [74] have been implicated to play a role in sexual development of plants. E3 ubiquitin ligase has a similar functional domain (i.e., RING-finger to Zip3), a meiosis specific gene and a member of the ZMM group known to play roles in crossover and SC assembly [54,75]. While this finding will require in depth analysis of these genes, the overall physiological and cellular processes of the GO analysis lend support to the sexual nature of the fused samples.

Understanding sexuality in Amoebozoa and future directions
We are just starting to scratch the surface of the enormous life cycle diversity in the supergroup Amoebozoa. The current study and previous studies indicated that members of Amoebozoa might use various mechanisms of sexual pathways that reflect their varied life cycles and behaviors [11,76]. For example, Dicytostelium has lost Spo11, but most amoebae seem to possess this gene [11,12,73]; this observation indicates the existence of variation in the mechanism of meiosis initiation in amoebozoans. A closer look of the known (yet to be described) diversity of the amoebozoans will likely uncover even additional, and novel, forms of life cycle and sexual strategies in this supergroup. Another layer of complexity in understanding sexuality in amoebae is the observation that some amoebae that display no signs of sexual activity constitutively, nonetheless express meiosis genes in actively growing cultures. Such an example is Acanthamoeba, a well-studied amoebozoan lineage that shows no evidence of fusion or other sexual like behavior during its life cycle [76]. Acanthamoeba and other amoebae are reported to change their ploidy during their life cycle [77][78][79][80][81]. This observation has led to a suggestion that amoebae might still be asexual lineages and have been using ancestral genes to perform recombination through other means such as polyploidization and gene conversion [82]. This study suggested that the genes comprising the meiosis toolkit had similar roles in recombination and mitosis in asexual lineages before the evolution of meiosis. Meiosis genes have been used as a blueprint for sex; however, it has been suggested that detection of meiosis genes should not be interpreted as evidence for sexual reproduction [82]. More investigation is required to unravel the origin and roles of meiosis genes in the life cycle of amoebzoans and other microbial eukaryotes.
Despite the progress we have made in the current study, many fundamental questions about the sexual nature and details of the life cycle of Cochliopodium remain unknown. Mating types, the nature of multiparental inheritance and polyploidization in Cochliopodium are still poorly elucidated. Cochliopodium has been shown to fuse even in monoclonal cultures (population derived from a single amoeba cell), though fusion frequency was lower compared to mixed cultures (S12 Fig). Mating types are known in Dictyostelium and other related amoebae [50,[83][84][85], further investigation is required to examine if fusion in Cochliopodium occurs randomly or among compatible mating types. The mechanism of inheritance involving multiple partners is a challenge to the classic mechanism of inheritance in dimorphic systems-well known in various eukaryotes. Triparental inheritance involving more than two mating types has been reported in Dictyostelium [50]. Further investigation using genome data, gene manipulation and cell biology studies are required to elucidate these fundamental questions in Cochliopodium and other amoebae showing similar life cycle behaviors.