physarum polycephalum (  slime mold  ) is a free - living amoebozoan protist . due to its rich cell biology , physarum is a model organism for cell motility1 and cell differentiation.24 more recently in rna editing,5 physarum is also a model organism for dna replication,6 optimization of cell morphology7 and behavioral intelligence.8 during its branched life cycle , physarum differentiates into several cell types that occur in temporal order ( fig . the multinucleate , macroscopic plasmodium , for example , may develop from microscopic , single nucleate amoebal cells or develop by germination of spherules.9 the number of nuclei contained in a plasmodium depends on its size . a plasmodium in a 9 cm petri dish , for example , contains approximately 10 nuclei that display natural synchrony in terms of cell cycle and differentiation status.10,11 during sporulation , the plasmodium differentiates into fruiting bodies containing spores which undergo a meiotic recombination to become the precursor cells of amoebae that may develop into gametes later on . sporulation occurs only in starved plasmodial cells and can be experimentally triggered , e.g. by exposure to far red light or to heat shock . well - established mendelian genetics and the occurrence of naturally synchronous macroscopic cells as source of homogeneous cell material for biochemical analyses have made physarum an attractive model organism for studying the regulatory control of cell differentiation . recent advances in sequencing technologies motivated detailed analyses of transcriptome dynamics during the life cycle of physarum.12,13 these studies were performed using rna isolated from cell populations . in general , differentiation follows spatial and temporal changes in transcript abundance in a cell type - specific manner . stochastic variations in gene expression presumably do impact cell - fate decisions,14 and therefore the time - resolved analysis of gene expression patterns in individual cells would provide valuable insight as compared to averaged data from measurements obtained from cell populations.14,15 expression patterns of single - cells have been analyzed using deep rna sequencing ( rna seq)16 to characterize the transcriptomes of individual embryonic mouse cells separated by technically complex procedures , and relying on the mouse genomic information for transcript assembly and mapping.17,18 the physarum genome has been sequenced but the data are still not assembled into a complete genome sequence ( the genome institute , washington university school of medicine ; http://genome.wustl.edu/genomes/view/physarum_polycephalum/ ) . therefore , here we evaluated the possibility of studying the global transcriptional changes during the differentiation of a single cell without relying on genomic information , and developed an approach to analyze the differential expression at several time points during the commitment of a plasmodial cell to sporulation . physarum macroplasmodia ( apogamic strain wt3119 ) were cultured as previously described.13,20 cells were grown and collected under two different conditions : ( i ) a plasmodium starved for 6 days ( competent d1 and d2 individual cell samples ) ; and ( ii ) a plasmodium starved for 6 days , exposed to far red light for 30 minutes , and returned to the dark for 6.5 hours ( l1 and l2 photoinduced cells ; table 1 ) . during this time period the cell becomes irreversibly committed to sporulation.21 samples were frozen with liquid nitrogen and poly(a)+ rna was isolated from the total rna samples ( by two rounds of oligo(dt ) affinity chromatography ) , and fragmented with ultra - sound ( 4 pulses of 30 sec at 4 c ) . subsequently , the rna fragments were poly(a)-tailed using poly(a ) polymerase , followed by treatment with tobacco acid pyrophosphatase ( tap ) . first - strand complementary dna ( cdna ) synthesis was performed using an oligo(dt)-adapter primer and the maloney murine leukemia virus ( m - mlv ) reverse transcriptase . the resulting cdnas were polymerase chain reaction ( pcr)-amplified for 1415 cycles to about 2030 ng/l , including distinctive 4-bp 5-barcodes for each sample ( table 1 ) , and using a high fidelity dna polymerase . the pcr products were purified with the agencourt ampure xp kit ( beckman coulter genomics ) , and pooled in equivalent amounts . cdnas in the range of 200 to 400 bp were fractionated from agarose gels and sequenced using the illumina hiseq 2000 system . poly(a)+ rna isolation , cdna synthesis and cdna library preparation as described here the 100-bp sequencing outputs were then trimmed for quality ( phred score > 33 ) , and later assembled de novo , using velvet22 and oases ( k - mers : 31 , 41 , 51 ) . the annotation of this assembly was carried out first through blast24 searches ( e - value 1e-3 ) against the swissprot25 protein database . a search for physarum noncoding rnas was not included due to the lack of complete gene models and a finished reference genome in this species . afterwards , domains and protein signature patterns were associated from matches to the interpro database , and gene ontology ( go ) based annotations were assigned using blast2go,26 from annotations pertaining to orthologs ( annotation e - value cutoff < 1e-6 ) . gene names and descriptions were filtered using the blast description annotation tool from blast2go.26 significant differences in go annotations between sets of up- and downregulated genes from each cdna library were evaluated using fisher exact tests , as implemented in blast2go . to assess the differential expression between the several single cells , the sequencing output was split using the barcode information for each sample . then the decoded outputs were mapped to the novel assembly with bowtie.27 samtools28 and tablet29 were later used to obtain mapped read counts . for expression comparisons , we obtained the following for each transcript : ( i ) the number of mapped reads , and ( ii ) the normalized expression value , as measured in reads per kilobase per million mapped reads ( rpkm).30 to identify differentially expressed transcripts between starved and photoinduced cells , the non - normalized mapped read count data was analyzed using the r - based package deseq.31 transcript abundances for each gene were estimated as a weighted mean of mapped read counts from each replicate , normalized to the library size . p - values ( adjusted for false discovery rate ) were generated for each gene in pairwise comparisons between different conditions ( competent and induced cells ) . four datasets consisting a total of 77.07 million 100-base reads from single cell physarum plasmodia were obtained from the illumina sequencing ( 77.02 million reads with phred score > 33 ; 7.12 gb ) . this rna - seq output was deposited in the european nucleotide archive32 ( study accession erp001220 ) . replicate data distributions were 1.85 and 1.82 gb corresponding to the starved plasmodium ( cdna library replicates d1 and d2 ) , and 1.67 and 1.78 gb for the cells collected 6 hours after photoinduction ( libraries l1 and l2 ; table 1 ) . therefore , assuming a 10% of protein - encoding genes , we obtained a 237.32x coverage for the 300 mb genome of physarum.13 later , from the de novo assembly of the trimmed sequencing output , we obtained 909,505 sequences , with a n50 contig size of 371 bp . large cdnas from this assembly ( > 500 bp ) were then clustered into 16,822 contigs ( n50 length : 778 bp ) with cap3 , to create a comprehensive set of representative transcripts . 10,278 of these contigs included transcript abundance data in the form of mapped reads in at least one cell sample , which was used as a measure of gene expression . transcriptionally active genes were then defined as contigs with at least one mapped read present in all four samples or differentiation stages ; in this regard , we observed 8,149 transcripts with mapped reads in all libraries . in order to assess the reproducibility of this single - cell approach , we normalized the mapped read counts for each transcriptionally active transcript as reads per kilobase per million mapped reads ( rpkm).30 we then compared the rpkm - normalized data for competent plasmodia ( d1 and d2 cells ) and photoinduced cells ( l1 and l2 ) that were separately processed ( fig . cells from related developmental stages exhibit very similar transcriptomic profiles ( competent cells : r = 0.99 ; light - induced plasmodia : r = 0.98 ) . on the other hand , lower correlations were found between competent and photoinduced cells ( r = 0.96 between d2 and l1 ; r = 0.97 in all other cases ; fig . 3 ) . however , further studies involving comparisons between several cell types and cell stages are required to establish whether these lower correlations account for the variations between individual cells . afterwards , we estimated the differentially expressed transcripts using deseq31 from the raw mapped reads . for these analyses , only contigs with a combined count of 300 mapped reads among all the samples were considered ; i.e. , 3,164 contigs were then selected that fitted these criteria . this mapped read count threshold was selected to reduce noise caused by spurious contigs and alignments . upon normalization , the distribution of mapped reads reflected the presence of differentially expressed transcripts , and genes with other kinds of regulation , with a slightly greater set of genes with higher expression in light - induced cells ( fig . 4 ) . specifically , we identified 556 upregulated transcripts ( p - value < 0.05 ; 504 of these with a false discovery rate ( fdr ) < 0.1 ) and 531 downregulated transcripts ( 475 with fdr < 0.1 ) between the photoinduced and competent cell libraries for transcriptionally active contigs with mapped reads ( fig . to assign functions to the novel sequences , annotations were associated to the transcriptome assembly . in this way we obtained 92,641 gene ontology ( go ) terms , corresponding to 5,722 swissprot orthologs , where 64,730 go terms belong to 4,222 sequences with mapped reads . cdnas were linked biological processes ( n = 1,135 ) , molecular functions ( n = 1,558 ) , and cellular components ( n = 576 ) . from the transcriptionally active genes with mapped reads , 231 annotated transcripts were upregulated , and 264 transcripts were downregulated ( supplementary data ) . a comparison of go annotations between sets of up- and downregulated transcripts revealed two terms exclusively featured in upregulated genes ( conversely , six go terms were identified only in downregulated transcripts (  synapse, go:0045202 ;  synapse part, go:0044456 ;  anti - oxidant activity, go:0016209 ;  translation regulator activity, go:0045182 ;  immune system process, go:0002376 ; and  viral reproduction, go:0016032 ; fig . next , we tested for the enrichment of go terms in up- and downregulated contigs against the full list of annotated transcripts , using the fisher s exact test as implemented in blast2go.26 in this manner , significant overrepresentation was found only in upregulated transcripts ( fdr = 0.037 ; p - value = 0.046 ) , with all go terms belonging to the molecular function category of ontologies : metal ion binding ( go:0046872 ) , calcium ion binding ( go:0005509 ) , ion binding ( go:0043167 ) , and cation binding ( go:0043169 ) . all these belong to the same hierarchy of ontologies , so all these can be summarized with the lower and more specific category , i.e. , the  calcium ion binding go term . both analyses of go annotations correlate well with previous studies that point to the upregulation of genes associated with the ion transport in the light - induced plasmodium.13,20 later , to evaluate the genes with similar regulation , we clustered the fully annotated transcripts for the highest statistically significant up- and down - regulation levels , as compared to the starved plasmodium cell libraries ( tables 2 and 3 ) . in spite of the apparent diversity of the annotations , we inferred potential functions based on ortholog identities and gene ontology assignments . in this way , we identified upregulated transcripts encoding endopeptidases ( physa ) , phospholipases ( pldg ) and stress response proteins ( bpm1 , nah1 ) , as well as genes related to biosynthetic processes ( coad , iod1 , pyr1 ) , development ( stx3 ) , chromatin remodeling ( ya27 ) , and signaling ( arf1 , cyh4 , sar1 , ltbp2 ) , that are highly expressed 6.5 hours after photoinduction ( table 2 ) . on the other hand , a different group of genes is downregulated upon light exposure . in this case , we found transcripts associated with actin - binding ( mys2 , coma ) , fmn - binding ( nos , ncpr ) , signaling ( vwka ) , sugar- ( tct1 ) and cation - binding proteins ( xanp , bot2 ) , a transporter ( pep3 ) and transferases ( set5 , hmnt ) , that were downregulated 6.5 hours after light induction ( table 3 ) . these measurements of transcriptional regulation at different time points correlate well with previous results in cell pools , where actin - binding and signaling proteins were identified as core members of the regulatory network during sporulation.13 we also noticed the expression of multiple transcript isoforms in the same cell at the same time point , both in upregulated transcripts ( pyr1 , bpm1 , physa ; table 2 ) , and in downregulated transcripts ( coma , vwka , nos , ncpr , xanp , mys2 , set5 ; table 3 ) . this phenomenon has also been observed in previous single - cell studies , and has been attributed to the complexity of transcript variants.33 whether these genes encode isoforms controlling stage - specific signaling pathways remains to be studied in detail . to date , two studies have reported the rna - seq analysis of transcriptomes in eukaryotic organisms , using single embryonic cells as models.17,18 in these works , both the assembly and mapping procedures were achieved using the mouse genome as a reference . previously , we succeeded in measuring the expression on a set of 35 differentially regulated genes in populations of single - cells , using the gexp rt - pcr platform.20 here , using the power of rna - seq to obtain whole transcriptomes without relying on previous genomic information , we were able to characterize a large set of expressed genes in different samples during the sporulation of physarum , an organism without a known genomic sequence . furthermore , in order to obtain single cells , all former studies on single - cell multiplex gene expression analysis required complex separation methods , such as pipetting cells manually , or using laser micro - dissection or fluorescence - activated cell sorting.15 in this study , we used the plasmodium , a natural macroscopic multinucleate single - cell stage from physarum , whose culture and handling is straightforward , and for which there are several well established methods for genetic manipulation.3 by combining the power of next - generation sequencing technologies with the simplicity for obtaining single cells from physarum , we developed an approach to characterize the whole transcriptome through the differentiation of this lower eukaryote at the single - cell level . the observed regulation patterns correlate well with previous studies on differential gene expression during the commitment to sporulation in the slime mold , particularly with respect to proteins involved in signaling and actin binding . we expect that improvements in single - cell transcriptomics , such as discrimination in sense and antisense transcripts , the ability to sequence a more diverse range of nucleic acid species , and other future developments , will help us to display a more precise picture of the regulatory network controlling differentiation in this organism . transcript abundance data and annotation (  supplementary.xls  ) . annotations and mapped reads for all the assembled cdnas ( excel spreadsheet ) .
physarum polycephalum is a unicellular eukaryote belonging to the amoebozoa group of organisms . 
 the complex life cycle involves various cell types that differ in morphology , function , and biochemical composition . 
 sporulation , one step in the life cycle , is a stimulus - controlled differentiation response of macroscopic plasmodial cells that develop into fruiting bodies . 
 well - established mendelian genetics and the occurrence of macroscopic cells with a naturally synchronous population of nuclei as source of homogeneous cell material for biochemical analyses make physarum an attractive model organism for studying the regulatory control of cell differentiation . here 
 , we develop an approach using rna - sequencing ( rna - seq ) , without needing to rely on a genome sequence as a reference , for studying the transcriptomic changes during stimulus - triggered commitment to sporulation in individual plasmodial cells . 
 the approach is validated through the obtained expression patterns and annotations , and particularly the results from up- and downregulated genes , which correlate well with previous studies .