BACKGROUND
spermatogenesis is a complex developmental process with a unique cell division, meiosis, as a major defining event.
the entire process includes maintenance of a small population of spermatogonial stem-cells, mitotic divisions of differentiating spermatogonia, meiotic prophase and ensuing divisions of spermatocytes, and post-meiotic differentiation of spermatids, by a process known as spermiogenesis.
in mammalian testes, spermatogenesis occurs within seminiferous tubules, where all germ cells associate with one kind of somatic cell, the sertoli cell, which provides the appropriate niche and microenvironment for the spermatogenic process.
the adult testis is characterized by presence of all of the cells types in the spermatogenic lineage, with waves of differentiation throughout the testis propelled by retinoic acid signaling  <cit> .
in mice, the first wave of spermatogenesis is initiated by spermatogonia shortly after birth, producing a sequential and orderly appearance of each of the more differentiated stages at regular intervals though the first four weeks of life.
although semi-synchronous with respect to the advancing wave of the most differentiated cells, the juvenile onset of spermatogenesis also includes regular and asynchronous initiation of subsequent waves of spermatogenic differentiation.
this asynchronous and continuous process of spermatogenesis has made it difficult to achieve molecular characterization of specific cell types in the lineage.
this has been particularly the case with respect to analysis of the defining process of gametogenesis, meiosis, which occurs in spermatocytes.
the complex events of meiosis i prophase include recombination, homologous chromosome pairing, and synapsis, taking place as spermatocytes progress through the leptotene, zygotene, pachytene and diplotene substages.
these events culminate in the first meiotic division, a reductive division in which homologous chromosomes are separated, producing secondary spermatocytes that rapidly undergo the second, equational, meiotic division to produce haploid round spermatids.
because of the genetic importance of meiotic recombination for the production of chromosomally normal gametes and offspring, there has been great interest in elucidating the molecular hallmarks and their underlying transcriptional signatures that define the meiotic spermatocyte substages.
toward the goal of achieving a molecular understanding of spermatogenesis, considerable effort has been devoted to separation of specific germ-cell differentiation substages from the histologically complex seminiferous epithelium.
a widely used approach, commonly known as the “sta-put” sedimentation process , involves enzymatic dissociation of germ and somatic cells and enrichment for specific stages by cell-size-based sedimentation at unit gravity on a bovine albumin gradient.
reasonably good enrichment of the most uniquely sized cells  can be obtained from testes of adult mice.
however, many cells of interest  are not retrieved from adults because of their anatomical position bounded by sertoli cell tight junctions.
while these early cell types can be retrieved from juvenile testes, the amount of biological material required is daunting and enrichment is not robust.
isolation of specific cell types by fluorescence-activated cell sorting, facs  is promising and becoming more widely applied; however sample sizes are small.
finally, rather than cell separation, total testis germ-cell populations can be collected in a developmental continuum during the first two to three weeks of juvenile development of mouse testes, to take advantage of the first wave of spermatogenesis.
in this way, molecular features have been defined with respect to the sequential appearance of more advanced cell stages of spermatogenesis.
this approach is useful primarily for correlating the appearance of a molecular entity to the developmental appearance of a specific cell type.
however, the degree of resolution has been suboptimal, because gene or protein expression has not been related to absolute frequencies of cell stages, a challenge we tackle in this study.
together, these methods for enriching or inferring spermatogenic  substages have contributed to studies over the past decade on the developmental transcriptome of mammalian spermatogenesis, as recently reviewed  <cit> .
most such studies have taken advantage of microarrays querying known coding sequences , and thus sequence-biased, but also provided interesting views of alternative splicing and other features of the spermatogenic transcriptome  <cit> , and identified previously unknown potential targets for contraception and fertility  <cit> .
more recently, methods for unbiased high-throughput deep sequencing of the transcriptome by rna-seq have been employed .
these studies have revealed unappreciated regulation of pirnas  <cit> , and global whole-genome views of the male germ-cell transcriptome  <cit> .
however, the challenge for all transcriptomic analyses, particularly rna-seq approaches, has been to computationally deconstruct the entire testis or germ-cell transcriptome into substage-specific transcriptomes.
this is an important goal, given both the abundance and complexity of rna species in the testis, particularly with respect to the imputed contribution of both coding and non-coding rna from spermatocytes and spermatids  <cit> .
in one example of such a computational approach, frequencies of specific germ-cell stages throughout the first wave of spermatogenesis previously published in another analysis  <cit>  were used to estimate cell type-specific expression patterns in a separate data set  <cit> .
while coming closer to the goal of substage-specific transcriptomes, this study relied on a low sample size and on integrating non-contemporaneous data.
other approaches  <cit> , employed cell sorting by facs and subsequent validation of purity by meiotic markers.
while these studies have yielded important insights into global gene expression switches, they rely on facs, not always available or practical for small samples, such as from infertile mutant models.
here we have tackled the problem of deconvolving transcriptomes of complex germ-cell populations into stage-specific transcriptomes by computationally integrating highly detailed and combinatorial cytological staging of the same cell samples as subjected to rna-seq analysis.
deep and accurate cell stage phenotyping was conducted using antibodies to stra <dig>  sycp <dig>  phosphorylated h2afx, and histh1t, all exhibiting meiotic substage specificity of expression and/or localization.
by collecting highly enriched germ-cell populations from testes at two-day intervals through the first wave of spermatogenesis, and by multiple sampling of individual mice  for both cytological composition and rna composition by high-throughput sequencing, we developed an unusually deep data set that portrays variation in cell substage composition and transcript abundance.
to decompose these data into substage-specific transcriptome patterns in an unbiased way, we chose a nonparametric solution with minimal assumptions about data structure, thus developing a novel permutation-based maximum covariance analysis.
this method enabled us to assign co-expressed transcripts to one or more meiotic substages, thereby linking distinct molecular programs to cytologically defined cell states.
moreover, to better understand the regulation of each germ-cell substage transcriptome, we integrated transcription factor information to identify some of the key molecular regulators driving these cell stages.
this approach provides a model for deconvolving transcriptomes of complex cell populations with well-defined cytological attributes.
together, these data provide an unprecedented view of the complexity of meiotic transcription programs and their coordinate regulation.
RESULTS
experimental design
germ cells were obtained from individual mice at two day intervals from  <dig> dpp to  <dig> dpp  and divided into two aliquots, one for cytology and one for rna-seq .
for the c57bl/6j  strain used in this study, this developmental window captures the initial differentiation of spermatogonia into early leptonema, through the stages of pachynema and diplonema.
for each sample, we determined the relative numbers of cells of each substage by cytological criteria ; by scoring germ cells from each individual mouse, we captured variability between mice.
the purity of germ cells in each sample, computed as germ-cell count divided by total dapi-stained cell count, was greater than 90 % for all samples.
upon inspection of both cytological and rna-seq data, two germ-cell samples collected at  <dig> dpp revealed poor concordance with age-matched samples, likely due to insufficient cell numbers for successful rna-seq library preparation ; thus these samples were omitted from all subsequent analyses .fig.
<dig> decomposition of meiotic cells and gene expression.
a experimental design.
germ cells were isolated from whole testes from each of five juvenile male mice at  <dig>   <dig>   <dig>   <dig>   <dig> and  <dig> dpp.
each sample was analyzed for both cytology and gene expression by rna-seq.
pmca was developed to identify the meiotic substage-specific transcriptomes.
b cytological classification and cellular populations.
isolated germ cells were immunolabeled for stage-specific maker proteins: stra <dig>  expressed in differentiating spermatogonia; sycp <dig>  a component of synaptonemal complex shows meiotic substage-specific labeling morphology; phosphorylated histone h2afx , marking dsbs throughout early prophase nuclei and restricted to the xy chromosomes during pachynema and diplonema; and histone histh1t, expressed in post-mid-pachytene spermatocytes.
nuclei were counter-stained with dapi.
c meiotic cell composition during the first wave of spermatogenesis.
the contribution by specific germ-cell stages for each developmental time point and mouse is shown, with colors representing specific meiotic substages.
abbreviations: sp’gonia = spermatogonia, prelep = pre-leptonema, early lep = early leptonema, lep = leptonema, late lep = late leptonema, zyg = zygonema, early pach = early pachynema, late pach = late pachynema, dip = diplonema.
d concordance of expression and cytological data.
principal component  <dig>  versus principal component  <dig>  from independent pca of cytological and rna-seq data.
colors correspond to time points and icon shape corresponds to data type as indicated
cytological deconvolution of the meiotic substage composition of germ cells during the first wave of spermatogenesis
germ-cell substage frequencies in each cytological preparation were determined by combinatorial immunolabeling with antibodies to meiotic substage marker proteins, scoring  <dig> germ cells per mouse.
marker proteins assessed were stra <dig>  a meiosis-initiating factor  <cit>  present in differentiated spermatogonia and some leptotene spermatocytes; sycp <dig>  a synaptonemal complex  protein present in the chromosomal axes of leptotene and zygotene spermatocytes, in the lateral elements of the mature sc in pachytene spermatocytes, and disassembling in diplotene spermatocytes  <cit> ; the phosphorylated form of histone h2afx , which localizes to chromatin surrounding dna double-strand breaks  in characteristic patterns that discriminate early meiotic prophase from the pachytene and diplotene stages  <cit> ; and the testis germ cell-specific histone h <dig> variant hist1h1t , which is a marker of the mid-to-late pachytene stage  <cit> .
the combinatorial labeling patterns for each marker protein allowed categorization of germ cells of each sample into eight substages, achieving a higher degree of meiotic substage discrimination than previous transcriptome analyses .table  <dig> immunolabeling criteria for cell types and meiotic substages
spermatogoniaa
aspermatogonia are a separate mitotic stage, not a meiotic substage
together, these data provide a comprehensive picture of postnatal spermatogenic progress through meiotic prophase of the first wave of differentiating cells for comparison to data previously obtained by histological analyses  <cit> .
at every time point in the current data, greater than 50 % of the retrieved germ cells were spermatogonia , reflecting continually initiating separate waves of spermatogenesis.
representation of these cells decreased over time , likely due to initiation of additional waves of meiotic cells, increasing numbers of post-spermatogonial spermatocytes, and establishment of cell junctions that impede cell retrieval.
throughout this period, due in part to initiation of subsequent waves of spermatogenesis, the average contribution to the total germ-cell population by preleptotene and leptotene spermatocytes remained relatively constant , thus reducing the power to identify transcripts specific to just these substages.
the first appreciable numbers of h1t-negative early pachytene spermatocytes were observed at  <dig> dpp, and h1t-positive late-pachytene cells were abundant by 16- <dig> dpp, consistent with previously published data  <cit> .
by  <dig> dpp, less than 10 % of the cells had progressed to the diplotene stage, and only about  <dig>  % to metaphase.
overall substage frequencies were similar at  <dig> and  <dig> dpp, and also at  <dig> and  <dig> dpp .
for computational analyses , we combined substages that exhibited similar frequency patterns across the developmental time span, specifically late leptotene and zygotene substages, and late pachytene and diplotene substages.
high concordance of gene expression and cytological data
to assess the feasibility of associating gene expression with cytologically-defined substages, we performed independent principal component analyses  on each data set.
the high concordance between the cytological frequency of cell substages and rna expression data  suggested that changes in gene expression in pooled germ cells can be explained by variation in cytological proportions of the differentiation states during spermatogenesis.
although cells in different substages may contribute different amounts of rna to each sample, this analysis demonstrates an overall quantitative agreement between the two data types.
this supported the validity of combining these two datasets to identify meiotic substage-specific transcriptional programs.
meiotic substage-specific gene expression derived by covariance analysis
to identify signatures for meiotic substage-specific gene expression, we developed a novel permutation-based maximum covariance analysis , which maps groups of co-expressed genes to combinatorial cytological marker-based staging based antibodies to stra <dig>  sycp <dig>  phosphorylated h2afx, and histh1t.
this statistical approach demonstrates concordance of distinct cellular programs to each meiotic substage based on the similarity of “cyto-pattern” and “gene-expression pattern” across samples.
in brief, the concordant patterns are derived from the most preeminent features of covariance between cytological and transcriptome data across all samples .
in this rna-seq dataset, a total of  <dig>  ensembl genes  and  <dig> noncode lncrna genes  <cit>  were detected  ≥ 3) in at least one sample of the isolated germ cells.
we also detected pirna precursors transcripts that had previously been identified  <cit>  with appropriate substages, although they are not the focus in the current study.
through pmca, we identified  <dig> spermatogonia genes and  <dig> meiosis substage-specific genes.
expression of many transcripts could not be assigned to distinct substages and instead, were shared across several consecutive substages;  <dig> genes were shared among late leptotene, zygotene, and early pachytene stages while  <dig> genes were shared among early pachytene, late pachytene, and diplotene stages .
notably, increased numbers of expressed genes were observed at  <dig> dpp when cells first reach the late pachytene stage .fig.
<dig> expression of substage- concordant and negatively concordant genes.
a heat maps show gene expression of substage-concordant genes at each time point.
b heat maps show gene expression of substage-negatively concordant genes at each time point.
color bars on the left side of the heat maps represent meiotic substages .
gene expression is shown as average log <dig> over time point replicates.
c overlap among substage-concordant gene lists, illustrating the number of transcripts in common between each pair of lists.
d overlap among substage-negatively concordant gene lists, illustrating the number of transcripts in common between each pair of lists
spermatogoniaa
aspermatogonia are a separate mitotic stage, not a meiotic substage
in addition to identifying genes concordant with specific substages, pmca also identified genes that were negatively concordant with specific substages.
the negatively concordant genes have patterns of expression that are opposite to the cytological patterns; thus negatively concordant defines gene expression that is lower when the cytological proportion of a specific substage is higher and vice versa .
genes detected as negatively concordant with one stage are often detected as concordant with another substage.
for example, genes negatively concordant with spermatogonia are genes expressed across the meiotic cell substages, but not enriched in spermatogonia.
not surprisingly, the  <dig> genes negatively concordant with spermatogonia were enriched for known meiotic genes.
indeed, among all genes expressed, many of those with known meiotic function and/or gene ontology annotation for meiotic function were not concordant with any single meiotic substage, but were expressed across meiotic substages.
validation of substage-specific gene expression patterns
we bolstered and validated these computationally derived results in three distinctly different ways.
first, gene expression analyses by qrt-pcr supported the rna-seq expression data.
second, the validity of the pmca procedure was queried by analysis of genes expressed in a highly enriched cell population of mid- to late-pachytene spermatocytes retrieved from adult testes by unit gravity sedimentation.
finally, we used our data to confirm a known pattern of male germ-cell specific gene expression, x-chromosome silencing during meiotic prophase.
we compared transcript expression patterns detected by rna-seq to independent qrt-pcr assays.
we tested representative genes among those increasing, decreasing, or showing little change in expression across the sequential time points, as well as genes known to be highly expressed in pachynema.
germ cells of  <dig> and  <dig> dpp males were subjected to qrt-pcr to determine expression of each tested gene relative to a reference housekeeping gene ; data are shown in additional file 2: figure s <dig>  the qrt-pcr assays reflected rna-seq findings.
for example, in  <dig> dpp samples we found elevated expression of pachytene-enriched genes and of genes shown by rna-seq to increase in expression throughout the juvenile period tested .
overall, patterns of gene expression are mostly concordant between the two different quantification methods and between sample sets, although, as expected, rna-seq provides finer resolution and higher information content.
we further validated the pmca-derived meiotic substage transcriptomes by comparison to highly enriched adult pachytene spermatocytes obtained by sedimentation at unit gravity , which sorts cells based on size, not cytology or dna content .
the purity of late pachytene/diplotene spermatocytes in this preparation was 90 % based on cytological immunostaining with antibodies to stage-specific proteins phosphorylated h2afx, sycp <dig>  and hist1h1t.
rna isolated from the size-enriched pachytene spermatocytes was subjected to rna-seq, allowing us to compare genes expressed in the pachytene spermatocytes to the meiotic substage-specific gene lists derived by pmca.
cross-referencing from this dataset, there was a highly significant enrichment of pachynema-expressed genes among the gene lists from later meiotic substages , but there was no enrichment for pachytene genes in the gene lists for spermatogonia or the early substages of meiosis .
in fact, 99 % of the genes in the late pachytene/diplotene list were also found in the enriched adult pachytene spermatocyte samples.table  <dig> representation of enriched pachytene spermatocyte genes in meiotic substages
spermatogoniaa
aspermatogonia are a separate mitotic stage, not a meiotic substage
finally, we queried the meiotic substage-specific patterns of gene expression determined by pmca by confirming meiotic sex-chromosome inactivation , a well-known feature of meiotic gene expression.
in spermatocytes, most of the axial length of x and y chromosomes, which are non-homologous, is not paired or synapsed.
by the onset of the pachytene stage of meiotic prophase, the unpaired regions of the sex chromosomes become transcriptionally inactivated  <cit> .
we assessed expression of x-linked genes from our rna-seq data; overall,  <dig> x-linked genes  were detected in the substage-associated gene lists.
of these, only  <dig> genes were detected in the highly enriched adult pachytene spermatocytes, providing evidence of robust msci.
these  <dig> transcripts, which do not show any specific regional localization on the x chromosome, might reflect ongoing transcription and escape from msci, but, perhaps more likely, may represent stable transcripts persisting after inactivation.
in the developmental transcriptome, the strongest x-linked gene signals are those represented in the negatively concordant gene lists , that is, genes with expression pattern opposite to the cytologically determined frequencies of each substage.
this strong signal is particularly apparent for genes in the late-pachytene negatively concordant list, which exhibit relatively constant expression from  <dig> to  <dig> dpp but then sharply decrease in expression at  <dig> and  <dig> dpp; this is evidence for pachytene substage msci.
in addition to msci being revealed by the negatively concordant gene lists, there is also diminished representation of x-linked genes in gene lists concordant with later stages of meiotic prophase .
overall,  <dig> x-linked genes appear to be down-regulated by pachynema to late pachynema, or even earlier  and only  <dig> x-linked genes with early expression are detected at later stages .table  <dig> substage-specific x-linked gene analysis
spermatogoniaa
aspermatogonia are a separate mitotic stage, not a meiotic substage
fig.
<dig> x-linked gene expression during the first wave of spermatogenesis.
a boxplots of the difference from mean gene expression at  <dig> dpp are shown for each time point, with genes on the autosomes in purple and those on the xy chromosomes in green.
gene expression is log <dig>  b heatmap of x-linked gene expression at each time point.
genes are clustered based on the pattern of expression: an increase of expression at  <dig> dpp and then stable throughout , steadily decreasing over time , or low expression from 8- <dig> dpp and high expression from 12- <dig> dpp followed by low expression at 16- <dig> dpp .
gene expression is shown as average log <dig> over time-point replicates
in conclusion, evidence from both msci and the transcriptome of highly enriched adult pachytene spermatocytes validates the pmca-derived substage- transcriptome signatures derived by pmca.
moreover, as discussed below, we find good concordance among our stage-specific gene lists and functions with those derived by other studies.
enrichment analysis of meiotic substage gene expression
the pmca-derived meiosis substage-specific gene lists, coupled with gene ontology  analysis  <cit> , can provide insight into meiotic and spermatogenic function.
genes assigned to each substage were ranked by their similarity score, a measure of how closely the gene’s expression pattern follows the substage’s cytological pattern.
using ranked lists, we performed a ranked go analysis in gorilla  <cit> ; these ranked lists provide approximations of meiotic substage gene ontology.
significantly, many meiosis-specific genes are not prominently represented among the meiosis substage-specific gene lists derived by pmca.
however, go terms for genes negatively concordant with spermatogonia were predominantly meiotic terms , confirming that meiotic genes tend to be expressed across multiple meiotic substages and not unique to a particular meiotic substage .
thus meiotic substage transitions are probably not acutely regulated at the level of transcription of many of the known meiotic genes, although among late leptonema/zygonema-associated go terms,  <dig> out of  <dig> go “biological process” terms were associated with meiosis regulation , and similar go terms were also identified in the early pachynema gene list  .
we extracted  <dig> genes for m. musculus that are associated with any go terms containing “meiosis” or “meiotic” , and of these genes we considered the  <dig> genes which are expressed in our developmental time series .
we found that  <dig>  genes were not concordant with any specific substage, while  <dig>  were concordant with late leptonema and zygonema,  <dig>  were concordant with early pachynema, and  <dig>  were concordant with late pachynema and diplonema.
of meiotic terms that were negatively concordant with substages,  <dig>  were negatively concordant with spermatogonia and  <dig>  were negatively concordant with late pachynema and diplonema .
go terms for genes expressed in late meiotic prophase were significantly enriched for spermatogenesis, spermiogenesis and fertilization, reflecting transcription of mrnas to be stored for later translation during the haploid phase of spermatogenesis  <cit> .
for example, early pachynema gene lists are enriched for go terms associated with spermiogenesis , and the majority of late meiotic prophase- or late-pachynema/diplonema-associated genes were associated with spermiogenesis or fertilization-related go terms  .
although genes concordant with early leptonema had transcription-related go terms in  <dig> out of  <dig> go terms for process , we found that many genes negatively concordant with leptotene through zygotene substages have go terms for transcription and related processes , possibly reflecting the cytologically diminished incorporation of rna precursors during the earliest meiotic prophase substages  <cit> .
also among the genes lists for these early stages we found extracellular membrane- or molecular transport-related go terms prominently represented , which may related to the fact that these cells transit through the sertoli cell junctions that create the blood-testis barrier  <cit> .
additionally we performed feature enrichment analysis with the hypergeometric test on the substage-specific gene lists.
interestingly, protein-coding genes are enriched in gene lists concordant with the late-leptotene/zygotene stages  and the early pachytene stage .
moreover, they are enriched in the set of genes that is negatively-concordant with late-pachytene and diplotene substages .
however, the genes concordant with late-pachytene and diplotene substages are not depleted in protein-coding gene, which comprise the expected majority of this gene set.
this suggests that while certain protein-coding genes are abundant in the late pach/dip substage, large numbers of other protein-coding genes are downregulated in this stage, perhaps reflecting cessation of mrna transcription in preparation for the meiotic division stage, or, alternatively, post-transcriptional regulation by piwi-interacting rnas .
recent studies show that pirnas play important roles in genome stability by suppressing harmful transposons as well as by regulating mrnas  <cit> , and future analyses could integrate pirna expression with these data.
transcription factor analysis
because our results point to a large and diverse meiotic germ-cell transcriptome, also noted by others  <cit> , we inferred the underlying regulatory networks accounting for these patterns, using the iregulon bioinformatic approach to identify transcription factors  potentially regulating substage-specific genes .
tfs were identified for each substage-specific gene list with high normalized enrichment scores , corresponding to an estimated false discovery rate of less than  <dig>   <cit> .
we then determined which tf genes were unique to, or shared among, the substages.
as can be seen in additional file 1: table s <dig> and additional file 2: figure s <dig>  mrna transcripts for some tfs are substage-specific.
for example, pou5f <dig> transcript is specific to the preleptotene stage while zbtb <dig> transcript is specific to the pachytene/diplotene substages.
other tf transcripts, similar to meiosis-specific genes in general, are shared across substages but may be specific to early or late meiosis ; zfp <dig> transcript is shared across late meiosis substages while tpbl <dig> transcript is common to the preleptotene  and late leptotene/zygotene/early pachytene  substages.
we also found that some tfs, for example mybl <dig>  also have target genes that are negatively concordant with substages .
we used the mgi bioinformatics batch query tool to determine that  <dig> of the  <dig> tfs in this analysis were previously found to be highly expressed specifically and/or significantly in the testis or male germ cells and  <dig> are annotated to male-reproduction-related phenotypes .
to gain a deeper insight into the regulatory patterns of meiosis, we selected the highest scoring tf genes for further analysis, based on iregulon’s nes score which corresponds to a low false discovery rate.
of these highest scoring tfs, half of them are annotated to meiosis-related functions: etv <dig>  e2f <dig>  gata <dig>  rfx <dig>  and zfp <dig>  we included four other well-known meiosis tfs: rarg, mybl <dig>  etv <dig>  and tbpl <dig>  we compared the mrna expression pattern of these tf genes with the expression patterns of their target genes in order to identify putative regulatory relationships, an analysis that is based on the assumption that a relevant tf protein appears more or less contemporaneously with its transcript .
in cases where the expression pattern of the tf gene had the opposite expression pattern of its target genes, we infer that the tf acts as a repressor on these targets; but in cases where the expression patterns of the tf gene and its target genes were concordant, we postulate that the tf enhances expression of its targets.
for tfs with genes that did not change over our time course, such as the known candidate nrf <dig>  we could not infer a relationship to target genes and excluded them from further analysis.
following these assumptions, as illustrated in figs.
<dig> and  <dig>  we suggest that ovol <dig> and yy <dig> act as repressors of target genes in early leptonema and etv <dig>  etv <dig>  and zfp <dig> act as repressors on target genes after initiation of meiosis, beginning with the late leptotene and zygotene substages.
evidence suggests that zbtb <dig>  gata <dig>  and etv <dig> also act as repressors on gene targets in pachynema ; while rfx <dig> appears to be an enhancer of target genes in late pachytene and diplotene substages , but may repress its target genes in preleptonema, although this relationship is less clear.
interestingly, these underlying assumptions based on relative expression levels suggest that mybl <dig>  tbpl <dig>  and e2f <dig> act as both enhancers and repressors of their target genes.
mybl <dig> and tbpl <dig> appear to repress target genes in early meiosis while enhancing target genes in late meiosis, while e2f <dig> has the opposite pattern and may repress target genes in late meiosis while enhancing target genes in early leptonema .
similar regulatory switching has been shown to result from changes in protein co-factors  <cit>  and post-translational modifications  <cit> .
moreover, genes for some of the tfs described above are also themselves target genes associated with specific substages.
by considering these relationships, we inferred candidate regulatory interactions among tfs; for example, zfp <dig> and etv <dig> suppress mybl <dig> gene expression, while mybl <dig> enhances expression of rfx <dig>  ovol <dig>  yy <dig>  and tbpl <dig> genes .fig.
<dig> expression of tfs and target genes for each substage.
the difference from mean gene expression at  <dig> dpp of the tf and its target genes at each time point.
to illustrate the overall pattern, a smoothed line was fit to the substage-specific gene expression.
color of lines represents the substage of target genes.
gene expression is average log <dig>  averaged over time point replicates.
expression is plotted for  zpf <dig>   etv <dig>   mybl <dig>  and  zbtb33
fig.
<dig> regulation of substage-concordant gene expression during meiotic prophase.
analyses suggest that substage-concordant gene expression is regulated both positively and negatively by a network of tfs.
tfs are shown by their gene names, and the meiotic substages of predicted tf target genes are shown in the center.
black arrows represent a tf acting as an enhancer on its target genes and red bars represent a tf acting as a suppressor on its target genes.
, lp/d = late pachynema and diplonema, anti-el = genes negatively concordant with early leptonema, anti-ll/z/ep = genes negatively concordant with late leptonema, zygonema, and early pachynema)
overall, 41 % of the genes associated with specific meiotic substages  are predicted targets of one or more of these eight meiosis-related tfs.
this noteworthy observation suggests that this concise regulatory network can account for a substantial portion of the meiotic program of transcription.
discussion
here we identify a male meiotic germ-cell transcriptome using a novel analysis based on a dense dataset of cytological substage-specific signatures.
the high concordance between rna-seq expression data and cytological proportions of isolated germ cells across all samples we analyzed allowed us to develop a novel pmca to identify the substage-specific transcriptomes for meiotic prophase.
this computational method does not require the use of facs or sedimentation sorting of cells and can be applied to other complex cell populations for which there are well-defined cytological criteria.
cytological deconvolution of rna-seq data from the first wave of spermatogenesis
in mammalian males, spermatogenesis is continuous and asynchronous, ensuring daily capacity to deliver sperm backed up by a testis comprised of abundant numbers of germ cells at all stages of spermatogenesis.
while advantageous to reproductively active males, this biological imperative has frustrated attempts to isolate stage-purified spermatogenic cells.
methods for stratification of cell samples by either sedimentation  <cit>  or flow cytometry  <cit>  enrich specific spermatogenesis cell stages, while the juvenile first wave of spermatogenesis  provides a leading edge of differentiating cells against a background of less differentiated cells.
in the laboratory mouse, meiosis is initiated at about  <dig> dpp in a subset of spermatogonia by retinoic acid stimulation of stra <dig>  <cit> ; however, initiation of subsequent waves of meiosis make the cellular population increasingly heterogeneous , resulting in a complex histology  <cit> .
the cell-stage frequencies presented here  are based on a powerful combinatorial application of antibodies recognizing well-characterized and highly stage-specific marker proteins.
we find a continuous and relatively constant background of the earliest spermatogenic cells, namely, spermatogonia and the preleptonema and leptonema , and, not surprisingly, these “background” spermatogenic cells are the ones most refractory for computational assignment of a substage-specific transcriptome.
on the other hand, pachytene spermatocytes are not present at the earliest juvenile stages, and appear in a discrete time period; consequently, it was possible to assign a robust gene list to this stage.
gene assignments to specific substages were guided by the covariance of transcriptomic data and the cytological data.
assignment of expressed transcripts to cellular subpopulations in heterogeneous samples has long been a computational challenge .
in this study, we exploited the advantage of paired rna-seq and cytological data to develop the novel pmca approach that isolates changes of gene expression specific to each meiotic substage.
maximum covariance analysis  was first developed in the meteorological sciences  <cit>  and was popularized in the climatological sciences in the 1990s  <cit> .
more recently, an mca approach has been used in a bioinformatics context to clarify relationships between gene and protein expression  <cit> .
while mca is often an effective tool for detecting common signals in two sets of variables, it can be limited by a tendency to fit spurious patterns when faced with increased sampling variation  <cit> .
current methods  <cit>  employ a parametric smoothing model using principal component regression, which requires a normality assumption or rely on go analysis  <cit> .
we developed a novel pmca method that not only overcomes the spurious pattern identification liability but that does so without the need for any parametric assumptions or reliance on go analysis.
our pmca approach is broadly applicable to multi-dimensional data derived from a common set of samples.
substage specificity of meiotic gene expression and regulation
historically, meiotic prophase substages are characterized by the morphology of chromatin, and correlated genetic mechanisms have been revealed in the past decades  <cit> .
it is known that there is widespread transcription of protein-coding genes in the testis  <cit> , and indeed, we found that among germ cells alone, more than  <dig>  genes are transcribed, suggesting that a significant portion of testis transcriptional complexity is due to germ cells.
using pmca, we successfully identified genes with expression patterns that matched the “cyto-patterns” of stage-specific cell frequencies determined by antibody labeling .
some genes were shared among substages; for example, late leptotene and zygotene shared  <dig> genes with early pachytene and early pachytene shared  <dig> genes with late pachytene and diplotene.
due to the strong temporal signal of contribution of the late pachytene germ cells, which do not appear until  <dig> dpp and greatly increase in representation by  <dig> dpp, we identified over  <dig> late pachytene genes.
however, since the proportion of preleptotene and early leptotene cells did not vary greatly over time, we were unable to deconvolve stage-specific transcriptomes for the earliest meiotic substages as confidently as for the later substages.
we compared these substage-specific gene lists to those developed in other recent analyses of the germ-cell testis transcriptome.
the striking increase in number of genes unique to the pachytene stage was also observed by soumillon, et al.
<cit> , where, in fact, greater numbers of genes were considered expressed than in this study.
however, we have used a more stringent cutoff  to eliminate from the data transcripts expressed at low levels.
when we relaxed this requirement and compared the number of genes expressed with a cutoff of tpm >  <dig>  we found the number of expressed transcripts to be similar to that reported by soumillion, et al.
<cit> .
we also compared the substage-specific gene lists to the relevant clusters identified by margolin, et al.
<cit> .
in their study, gene lists were temporally clustered from single samples taken at  <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig>  and  <dig> dpp.
while our study was based on multiple replicates at each time point as well as fine-grained cytological analysis, our pmca gene lists substantially overlapped with their derived gene list clusters .
in all,  <dig> of the  <dig> genes in clusters identified by margolin, et al.
<cit>  were also expressed in our dataset.
we also compared our results with those derived from flow-cytometry based methods of sorting by dna content  <cit> .
for both of these studies, we found overall concordance in genes assigned to meiotic substages .
these similarities were especially pronounced for late pachytene and diplotene substages due to the large numbers of genes in these sets.
however, we also observed a broader trend of alignment between early and late substages.
furthermore, our use of multiple markers to cytologically characterize cells paired with our pmca analysis allowed for more precise substage assignment than possible by cell sorting.
for instance, our spermatogonia and pre-leptotene sets were the top two strongest overlaps with the “secondary spermatocyte”  cell fraction in da cruz, et al.
<cit>  , but split these assignments to provide improved substage resolution.
similarly, our joint late-leptotene and zygotene genes and early pachytene genes had the two greatest overlaps with the combined “leptotene and zygotene” and “pachytene spermatocytes”  genes , and further partition these transcripts into more precise substages.
analysis of fallahi et al.
<cit>  data was limited by different experimental platforms , different assay timing , and generally low numbers of uniquely assigned genes, but also revealed significant overlaps in late substages .
by assessing functional concordance by similarity in go annotations, we determined that our early leptotene and fallahi et al.
<cit>  pre-leptotene sets share a common enrichment in rna transcription genes.
finally, the most prominent divergences between our data and these previous results were in the spermatogonia fraction.
we note that cell sorting techniques are susceptible to inclusion of somatic-cells in this fraction, whereas our spermatogonia genes were significantly negatively-concordant with many standard meiosis genes .
in sum, these comparisons suggest that our transcript sets encompass and substantially expand these previous findings.
in this study, we characterized each meiotic substage using the list of substage-specific genes from pmca analysis.
some meiosis genes were associated to particular substages.
for example, spata <dig>  <cit>  and msh4/ <dig>  <cit>  are both highly associated with late leptotene/zygotene substages, and both are required for recombination.
however, many canonical meiosis genes were found throughout all meiotic prophase substages; these genes include rad <dig>  rec <dig>  and syce <dig>  this may well reflect a lack of acute transcriptional regulation for these important transcripts.
rather than transcription at the precise meiotic substage of use, quality and efficiency of meiosis may be ensured by having transcripts present and available for translation throughout meiotic prophase, with substage transitions regulated post-transcriptionally and/or post-translationally.
further, although we confirmed that spermatogonia negatively concordant genes are enriched with meiotic genes, we found several meiotic genes in the gene list concordant with spermatogonia, including fign, which was reported to be required for meiotic recombination  <cit> , suggesting that some of the meiotic program is set up prior to meiotic entry.
not surprisingly, this developmental transcriptome analysis also revealed that many genes expressed during meiosis do not have known functions directly contributing to meiosis.
these may instead be part of a parallel program of spermatogenic gene expression.
for example, many of the genes expressed in early leptonema are associated with transcription or rna metabolism, as well as with cellular processes such as cell-cell interactions, which are of considerable importance in the biology of the seminiferous tubule.
pmca also revealed that late pachytene/diplotene-expressed genes are significantly associated with go terms associated with spermiogenic processes, correlating well with findings of other studies  <cit> .
during post-meiotic stages of nuclear condensation, transcription is globally repressed  <cit> , underlying the biological rationale for prior transcription to support the active protein synthesis during spermiogenesis.
in addition, the late pachytene/diplotene-associated gene lists are significantly enriched with non-protein-coding rnas, including lincrnas  and pirna precursors.
meiotic regulatory networks
analysis of predicted tf targets  <cit>  generated a regulatory network that potentially governs the meiotic and spermatogenic programs of gene expression identified by this study .
the protein zfp <dig> was central in our network, and is required for embryo development in zebrafish  <cit>  and human  <cit> .
the human ortholog, znf <dig>  is ubiquitously expressed  <cit> , and binds to the promoter region of target genes where it is required for formation of chromatin loops by interacting with polii and rad <dig>  <cit> .
because rad <dig> is differentially expressed after the pachytene stage of meiosis  <cit> , it is possible that zfp <dig> is crucial for spermatogenic processes by regulating the transcription of genes during pachytene/diplotene substage, an idea that should be tested experimentally since our analysis indicates that zfp <dig> may target genes strongly associated with genes in these substages.
reinforcing the validity of our approach, our regulatory network also involved mybl <dig> .
our results support previous findings that a subset of late pachynema/diplonema-expressed genes, many involved with the spermatogenic and spermiogenic processes, are associated with mybl <dig>  <cit> .
mybl <dig> regulates transcription of cell cycle-related genes.
germ cells with a homozygous mutation in the mybl <dig> gene exhibit defects in meiotic chromosome synapsis  <cit> .
moreover, our results suggest that mybl <dig> also is in a network associated with the transcription of genes encoding proteins required for pirna biogenesis, including tdrd <dig>  tdrd <dig> and tdrd <dig>  consistent with previous reports  <cit> .
going beyond our core regulatory network, we also identified several other tfs, whose target genes include the tdrd family and piwi genes, potentially involved in pirna processing: atf <dig>  elf <dig>  elk <dig>  elk <dig>  fli <dig>  nfya, nfyb, nrf <dig>  pbx <dig>  rfx <dig>  rfx <dig>  sp <dig>  yy <dig>  zbtb <dig>  zfp <dig> and zfp <dig>  this analysis not only revealed a core network for transcriptional regulation of meiotic progression  but also suggested that a significant proportion of the genes expressed in the meiotic transcriptome may be controlled by a concise entourage of transcription factors.
CONCLUSIONS
this study has untangled in part the complexity and parallel process of spermatogenesis and meiosis by focusing on associating gene expression with highly specific cytological signatures defining meiotic prophase substages.
this unique and powerful approach to deconvolving transcriptomes of complex cell populations is applicable for discovery of transcriptional signals in other such complex cell populations or heterogeneous tissues.

 spermatogenesis is a complex developmental process with a unique cell division, meiosis, as a major defining event. the entire process includes maintenance of a small population of spermatogonial stem-cells, mitotic divisions of differentiating spermatogonia, meiotic prophase and ensuing divisions of spermatocytes, and post-meiotic differentiation of spermatids, by a process known as spermiogenesis. in mammalian testes, spermatogenesis occurs within seminiferous tubules, where all germ cells associate with one kind of somatic cell, the sertoli cell, which provides the appropriate niche and microenvironment for the spermatogenic process. the adult testis is characterized by presence of all of the cells types in the spermatogenic lineage, with waves of differentiation throughout the testis propelled by retinoic acid signaling. in mice, the first wave of spermatogenesis is initiated by spermatogonia shortly after birth, producing a sequential and orderly appearance of each of the more differentiated stages at regular intervals though the first four weeks of life. although semi-synchronous with respect to the advancing wave of the most differentiated cells, the juvenile onset of spermatogenesis also includes regular and asynchronous initiation of subsequent waves of spermatogenic differentiation. this asynchronous and continuous process of spermatogenesis has made it difficult to achieve molecular characterization of specific cell types in the lineage. this has been particularly the case with respect to analysis of the defining process of gametogenesis, meiosis, which occurs in spermatocytes. the complex events of meiosis i prophase include recombination, homologous chromosome pairing, and synapsis, taking place as spermatocytes progress through the leptotene, zygotene, pachytene and diplotene substages. these events culminate in the first meiotic division, a reductive division in which homologous chromosomes are separated, producing secondary spermatocytes that rapidly undergo the second, equational, meiotic division to produce haploid round spermatids. because of the genetic importance of meiotic recombination for the production of chromosomally normal gametes and offspring, there has been great interest in elucidating the molecular hallmarks and their underlying transcriptional signatures that define the meiotic spermatocyte substages. toward the goal of achieving a molecular understanding of spermatogenesis, considerable effort has been devoted to separation of specific germ-cell differentiation substages from the histologically complex seminiferous epithelium. a widely used approach, commonly known as the “sta-put” sedimentation process, involves enzymatic dissociation of germ and somatic cells and enrichment for specific stages by cell-size-based sedimentation at unit gravity on a bovine albumin gradient. reasonably good enrichment of the most uniquely sized cells can be obtained from testes of adult mice. however, many cells of interest are not retrieved from adults because of their anatomical position bounded by sertoli cell tight junctions. while these early cell types can be retrieved from juvenile testes, the amount of biological material required is daunting and enrichment is not robust. isolation of specific cell types by fluorescence-activated cell sorting, facs is promising and becoming more widely applied; however sample sizes are small. finally, rather than cell separation, total testis germ-cell populations can be collected in a developmental continuum during the first two to three weeks of juvenile development of mouse testes, to take advantage of the first wave of spermatogenesis. in this way, molecular features have been defined with respect to the sequential appearance of more advanced cell stages of spermatogenesis. this approach is useful primarily for correlating the appearance of a molecular entity to the developmental appearance of a specific cell type. however, the degree of resolution has been suboptimal, because gene or protein expression has not been related to absolute frequencies of cell stages, a challenge we tackle in this study. together, these methods for enriching or inferring spermatogenic substages have contributed to studies over the past decade on the developmental transcriptome of mammalian spermatogenesis, as recently reviewed. most such studies have taken advantage of microarrays querying known coding sequences, and thus sequence-biased, but also provided interesting views of alternative splicing and other features of the spermatogenic transcriptome, and identified previously unknown potential targets for contraception and fertility. more recently, methods for unbiased high-throughput deep sequencing of the transcriptome by rna-seq have been employed. these studies have revealed unappreciated regulation of pirnas, and global whole-genome views of the male germ-cell transcriptome. however, the challenge for all transcriptomic analyses, particularly rna-seq approaches, has been to computationally deconstruct the entire testis or germ-cell transcriptome into substage-specific transcriptomes. this is an important goal, given both the abundance and complexity of rna species in the testis, particularly with respect to the imputed contribution of both coding and non-coding rna from spermatocytes and spermatids. in one example of such a computational approach, frequencies of specific germ-cell stages throughout the first wave of spermatogenesis previously published in another analysis were used to estimate cell type-specific expression patterns in a separate data set. while coming closer to the goal of substage-specific transcriptomes, this study relied on a low sample size and on integrating non-contemporaneous data. other approaches, employed cell sorting by facs and subsequent validation of purity by meiotic markers. while these studies have yielded important insights into global gene expression switches, they rely on facs, not always available or practical for small samples, such as from infertile mutant models. here we have tackled the problem of deconvolving transcriptomes of complex germ-cell populations into stage-specific transcriptomes by computationally integrating highly detailed and combinatorial cytological staging of the same cell samples as subjected to rna-seq analysis. deep and accurate cell stage phenotyping was conducted using antibodies to stra sycp phosphorylated h2afx, and histh1t, all exhibiting meiotic substage specificity of expression and/or localization. by collecting highly enriched germ-cell populations from testes at two-day intervals through the first wave of spermatogenesis, and by multiple sampling of individual mice for both cytological composition and rna composition by high-throughput sequencing, we developed an unusually deep data set that portrays variation in cell substage composition and transcript abundance. to decompose these data into substage-specific transcriptome patterns in an unbiased way, we chose a nonparametric solution with minimal assumptions about data structure, thus developing a novel permutation-based maximum covariance analysis. this method enabled us to assign co-expressed transcripts to one or more meiotic substages, thereby linking distinct molecular programs to cytologically defined cell states. moreover, to better understand the regulation of each germ-cell substage transcriptome, we integrated transcription factor information to identify some of the key molecular regulators driving these cell stages. this approach provides a model for deconvolving transcriptomes of complex cell populations with well-defined cytological attributes. together, these data provide an unprecedented view of the complexity of meiotic transcription programs and their coordinate regulation.  experimental design germ cells were obtained from individual mice at two day intervals from dpp to dpp and divided into two aliquots, one for cytology and one for rna-seq. for the c57bl/6j strain used in this study, this developmental window captures the initial differentiation of spermatogonia into early leptonema, through the stages of pachynema and diplonema. for each sample, we determined the relative numbers of cells of each substage by cytological criteria; by scoring germ cells from each individual mouse, we captured variability between mice. the purity of germ cells in each sample, computed as germ-cell count divided by total dapi-stained cell count, was greater than 90 % for all samples. upon inspection of both cytological and rna-seq data, two germ-cell samples collected at dpp revealed poor concordance with age-matched samples, likely due to insufficient cell numbers for successful rna-seq library preparation; thus these samples were omitted from all subsequent analyses.fig. decomposition of meiotic cells and gene expression. a experimental design. germ cells were isolated from whole testes from each of five juvenile male mice at and dpp. each sample was analyzed for both cytology and gene expression by rna-seq. pmca was developed to identify the meiotic substage-specific transcriptomes. b cytological classification and cellular populations. isolated germ cells were immunolabeled for stage-specific maker proteins: stra expressed in differentiating spermatogonia; sycp a component of synaptonemal complex shows meiotic substage-specific labeling morphology; phosphorylated histone h2afx, marking dsbs throughout early prophase nuclei and restricted to the xy chromosomes during pachynema and diplonema; and histone histh1t, expressed in post-mid-pachytene spermatocytes. nuclei were counter-stained with dapi. c meiotic cell composition during the first wave of spermatogenesis. the contribution by specific germ-cell stages for each developmental time point and mouse is shown, with colors representing specific meiotic substages. abbreviations: sp’gonia = spermatogonia, prelep = pre-leptonema, early lep = early leptonema, lep = leptonema, late lep = late leptonema, zyg = zygonema, early pach = early pachynema, late pach = late pachynema, dip = diplonema. d concordance of expression and cytological data. principal component versus principal component from independent pca of cytological and rna-seq data. colors correspond to time points and icon shape corresponds to data type as indicated cytological deconvolution of the meiotic substage composition of germ cells during the first wave of spermatogenesis germ-cell substage frequencies in each cytological preparation were determined by combinatorial immunolabeling with antibodies to meiotic substage marker proteins, scoring germ cells per mouse. marker proteins assessed were stra a meiosis-initiating factor present in differentiated spermatogonia and some leptotene spermatocytes; sycp a synaptonemal complex protein present in the chromosomal axes of leptotene and zygotene spermatocytes, in the lateral elements of the mature sc in pachytene spermatocytes, and disassembling in diplotene spermatocytes; the phosphorylated form of histone h2afx, which localizes to chromatin surrounding dna double-strand breaks in characteristic patterns that discriminate early meiotic prophase from the pachytene and diplotene stages; and the testis germ cell-specific histone h variant hist1h1t, which is a marker of the mid-to-late pachytene stage. the combinatorial labeling patterns for each marker protein allowed categorization of germ cells of each sample into eight substages, achieving a higher degree of meiotic substage discrimination than previous transcriptome analyses.table immunolabeling criteria for cell types and meiotic substages spermatogoniaa aspermatogonia are a separate mitotic stage, not a meiotic substage together, these data provide a comprehensive picture of postnatal spermatogenic progress through meiotic prophase of the first wave of differentiating cells for comparison to data previously obtained by histological analyses. at every time point in the current data, greater than 50 % of the retrieved germ cells were spermatogonia, reflecting continually initiating separate waves of spermatogenesis. representation of these cells decreased over time, likely due to initiation of additional waves of meiotic cells, increasing numbers of post-spermatogonial spermatocytes, and establishment of cell junctions that impede cell retrieval. throughout this period, due in part to initiation of subsequent waves of spermatogenesis, the average contribution to the total germ-cell population by preleptotene and leptotene spermatocytes remained relatively constant, thus reducing the power to identify transcripts specific to just these substages. the first appreciable numbers of h1t-negative early pachytene spermatocytes were observed at dpp, and h1t-positive late-pachytene cells were abundant by 16- dpp, consistent with previously published data. by dpp, less than 10 % of the cells had progressed to the diplotene stage, and only about  % to metaphase. overall substage frequencies were similar at and dpp, and also at and dpp. for computational analyses, we combined substages that exhibited similar frequency patterns across the developmental time span, specifically late leptotene and zygotene substages, and late pachytene and diplotene substages. high concordance of gene expression and cytological data to assess the feasibility of associating gene expression with cytologically-defined substages, we performed independent principal component analyses on each data set. the high concordance between the cytological frequency of cell substages and rna expression data suggested that changes in gene expression in pooled germ cells can be explained by variation in cytological proportions of the differentiation states during spermatogenesis. although cells in different substages may contribute different amounts of rna to each sample, this analysis demonstrates an overall quantitative agreement between the two data types. this supported the validity of combining these two datasets to identify meiotic substage-specific transcriptional programs. meiotic substage-specific gene expression derived by covariance analysis to identify signatures for meiotic substage-specific gene expression, we developed a novel permutation-based maximum covariance analysis, which maps groups of co-expressed genes to combinatorial cytological marker-based staging based antibodies to stra sycp phosphorylated h2afx, and histh1t. this statistical approach demonstrates concordance of distinct cellular programs to each meiotic substage based on the similarity of “cyto-pattern” and “gene-expression pattern” across samples. in brief, the concordant patterns are derived from the most preeminent features of covariance between cytological and transcriptome data across all samples. in this rna-seq dataset, a total of ensembl genes and noncode lncrna genes were detected  ≥ 3) in at least one sample of the isolated germ cells. we also detected pirna precursors transcripts that had previously been identified with appropriate substages, although they are not the focus in the current study. through pmca, we identified spermatogonia genes and meiosis substage-specific genes. expression of many transcripts could not be assigned to distinct substages and instead, were shared across several consecutive substages; genes were shared among late leptotene, zygotene, and early pachytene stages while genes were shared among early pachytene, late pachytene, and diplotene stages. notably, increased numbers of expressed genes were observed at dpp when cells first reach the late pachytene stage.fig. expression of substage- concordant and negatively concordant genes. a heat maps show gene expression of substage-concordant genes at each time point. b heat maps show gene expression of substage-negatively concordant genes at each time point. color bars on the left side of the heat maps represent meiotic substages. gene expression is shown as average log over time point replicates. c overlap among substage-concordant gene lists, illustrating the number of transcripts in common between each pair of lists. d overlap among substage-negatively concordant gene lists, illustrating the number of transcripts in common between each pair of lists spermatogoniaa aspermatogonia are a separate mitotic stage, not a meiotic substage in addition to identifying genes concordant with specific substages, pmca also identified genes that were negatively concordant with specific substages. the negatively concordant genes have patterns of expression that are opposite to the cytological patterns; thus negatively concordant defines gene expression that is lower when the cytological proportion of a specific substage is higher and vice versa. genes detected as negatively concordant with one stage are often detected as concordant with another substage. for example, genes negatively concordant with spermatogonia are genes expressed across the meiotic cell substages, but not enriched in spermatogonia. not surprisingly, the genes negatively concordant with spermatogonia were enriched for known meiotic genes. indeed, among all genes expressed, many of those with known meiotic function and/or gene ontology annotation for meiotic function were not concordant with any single meiotic substage, but were expressed across meiotic substages. validation of substage-specific gene expression patterns we bolstered and validated these computationally derived results in three distinctly different ways. first, gene expression analyses by qrt-pcr supported the rna-seq expression data. second, the validity of the pmca procedure was queried by analysis of genes expressed in a highly enriched cell population of mid- to late-pachytene spermatocytes retrieved from adult testes by unit gravity sedimentation. finally, we used our data to confirm a known pattern of male germ-cell specific gene expression, x-chromosome silencing during meiotic prophase. we compared transcript expression patterns detected by rna-seq to independent qrt-pcr assays. we tested representative genes among those increasing, decreasing, or showing little change in expression across the sequential time points, as well as genes known to be highly expressed in pachynema. germ cells of and dpp males were subjected to qrt-pcr to determine expression of each tested gene relative to a reference housekeeping gene; data are shown in additional file 2: figure s the qrt-pcr assays reflected rna-seq findings. for example, in dpp samples we found elevated expression of pachytene-enriched genes and of genes shown by rna-seq to increase in expression throughout the juvenile period tested. overall, patterns of gene expression are mostly concordant between the two different quantification methods and between sample sets, although, as expected, rna-seq provides finer resolution and higher information content. we further validated the pmca-derived meiotic substage transcriptomes by comparison to highly enriched adult pachytene spermatocytes obtained by sedimentation at unit gravity, which sorts cells based on size, not cytology or dna content. the purity of late pachytene/diplotene spermatocytes in this preparation was 90 % based on cytological immunostaining with antibodies to stage-specific proteins phosphorylated h2afx, sycp and hist1h1t. rna isolated from the size-enriched pachytene spermatocytes was subjected to rna-seq, allowing us to compare genes expressed in the pachytene spermatocytes to the meiotic substage-specific gene lists derived by pmca. cross-referencing from this dataset, there was a highly significant enrichment of pachynema-expressed genes among the gene lists from later meiotic substages, but there was no enrichment for pachytene genes in the gene lists for spermatogonia or the early substages of meiosis. in fact, 99 % of the genes in the late pachytene/diplotene list were also found in the enriched adult pachytene spermatocyte samples.table representation of enriched pachytene spermatocyte genes in meiotic substages spermatogoniaa aspermatogonia are a separate mitotic stage, not a meiotic substage finally, we queried the meiotic substage-specific patterns of gene expression determined by pmca by confirming meiotic sex-chromosome inactivation, a well-known feature of meiotic gene expression. in spermatocytes, most of the axial length of x and y chromosomes, which are non-homologous, is not paired or synapsed. by the onset of the pachytene stage of meiotic prophase, the unpaired regions of the sex chromosomes become transcriptionally inactivated. we assessed expression of x-linked genes from our rna-seq data; overall, x-linked genes were detected in the substage-associated gene lists. of these, only genes were detected in the highly enriched adult pachytene spermatocytes, providing evidence of robust msci. these transcripts, which do not show any specific regional localization on the x chromosome, might reflect ongoing transcription and escape from msci, but, perhaps more likely, may represent stable transcripts persisting after inactivation. in the developmental transcriptome, the strongest x-linked gene signals are those represented in the negatively concordant gene lists, that is, genes with expression pattern opposite to the cytologically determined frequencies of each substage. this strong signal is particularly apparent for genes in the late-pachytene negatively concordant list, which exhibit relatively constant expression from to dpp but then sharply decrease in expression at and dpp; this is evidence for pachytene substage msci. in addition to msci being revealed by the negatively concordant gene lists, there is also diminished representation of x-linked genes in gene lists concordant with later stages of meiotic prophase. overall, x-linked genes appear to be down-regulated by pachynema to late pachynema, or even earlier and only x-linked genes with early expression are detected at later stages.table substage-specific x-linked gene analysis spermatogoniaa aspermatogonia are a separate mitotic stage, not a meiotic substage fig. x-linked gene expression during the first wave of spermatogenesis. a boxplots of the difference from mean gene expression at dpp are shown for each time point, with genes on the autosomes in purple and those on the xy chromosomes in green. gene expression is log b heatmap of x-linked gene expression at each time point. genes are clustered based on the pattern of expression: an increase of expression at dpp and then stable throughout, steadily decreasing over time, or low expression from 8- dpp and high expression from 12- dpp followed by low expression at 16- dpp. gene expression is shown as average log over time-point replicates in conclusion, evidence from both msci and the transcriptome of highly enriched adult pachytene spermatocytes validates the pmca-derived substage- transcriptome signatures derived by pmca. moreover, as discussed below, we find good concordance among our stage-specific gene lists and functions with those derived by other studies. enrichment analysis of meiotic substage gene expression the pmca-derived meiosis substage-specific gene lists, coupled with gene ontology analysis, can provide insight into meiotic and spermatogenic function. genes assigned to each substage were ranked by their similarity score, a measure of how closely the gene’s expression pattern follows the substage’s cytological pattern. using ranked lists, we performed a ranked go analysis in gorilla; these ranked lists provide approximations of meiotic substage gene ontology. significantly, many meiosis-specific genes are not prominently represented among the meiosis substage-specific gene lists derived by pmca. however, go terms for genes negatively concordant with spermatogonia were predominantly meiotic terms, confirming that meiotic genes tend to be expressed across multiple meiotic substages and not unique to a particular meiotic substage. thus meiotic substage transitions are probably not acutely regulated at the level of transcription of many of the known meiotic genes, although among late leptonema/zygonema-associated go terms, out of go “biological process” terms were associated with meiosis regulation, and similar go terms were also identified in the early pachynema gene list. we extracted genes for m. musculus that are associated with any go terms containing “meiosis” or “meiotic”, and of these genes we considered the genes which are expressed in our developmental time series. we found that genes were not concordant with any specific substage, while were concordant with late leptonema and zygonema, were concordant with early pachynema, and were concordant with late pachynema and diplonema. of meiotic terms that were negatively concordant with substages, were negatively concordant with spermatogonia and were negatively concordant with late pachynema and diplonema. go terms for genes expressed in late meiotic prophase were significantly enriched for spermatogenesis, spermiogenesis and fertilization, reflecting transcription of mrnas to be stored for later translation during the haploid phase of spermatogenesis. for example, early pachynema gene lists are enriched for go terms associated with spermiogenesis, and the majority of late meiotic prophase- or late-pachynema/diplonema-associated genes were associated with spermiogenesis or fertilization-related go terms. although genes concordant with early leptonema had transcription-related go terms in out of go terms for process, we found that many genes negatively concordant with leptotene through zygotene substages have go terms for transcription and related processes, possibly reflecting the cytologically diminished incorporation of rna precursors during the earliest meiotic prophase substages. also among the genes lists for these early stages we found extracellular membrane- or molecular transport-related go terms prominently represented, which may related to the fact that these cells transit through the sertoli cell junctions that create the blood-testis barrier. additionally we performed feature enrichment analysis with the hypergeometric test on the substage-specific gene lists. interestingly, protein-coding genes are enriched in gene lists concordant with the late-leptotene/zygotene stages and the early pachytene stage. moreover, they are enriched in the set of genes that is negatively-concordant with late-pachytene and diplotene substages. however, the genes concordant with late-pachytene and diplotene substages are not depleted in protein-coding gene, which comprise the expected majority of this gene set. this suggests that while certain protein-coding genes are abundant in the late pach/dip substage, large numbers of other protein-coding genes are downregulated in this stage, perhaps reflecting cessation of mrna transcription in preparation for the meiotic division stage, or, alternatively, post-transcriptional regulation by piwi-interacting rnas. recent studies show that pirnas play important roles in genome stability by suppressing harmful transposons as well as by regulating mrnas, and future analyses could integrate pirna expression with these data. transcription factor analysis because our results point to a large and diverse meiotic germ-cell transcriptome, also noted by others, we inferred the underlying regulatory networks accounting for these patterns, using the iregulon bioinformatic approach to identify transcription factors potentially regulating substage-specific genes. tfs were identified for each substage-specific gene list with high normalized enrichment scores, corresponding to an estimated false discovery rate of less than. we then determined which tf genes were unique to, or shared among, the substages. as can be seen in additional file 1: table s and additional file 2: figure s mrna transcripts for some tfs are substage-specific. for example, pou5f transcript is specific to the preleptotene stage while zbtb transcript is specific to the pachytene/diplotene substages. other tf transcripts, similar to meiosis-specific genes in general, are shared across substages but may be specific to early or late meiosis; zfp transcript is shared across late meiosis substages while tpbl transcript is common to the preleptotene and late leptotene/zygotene/early pachytene substages. we also found that some tfs, for example mybl also have target genes that are negatively concordant with substages. we used the mgi bioinformatics batch query tool to determine that of the tfs in this analysis were previously found to be highly expressed specifically and/or significantly in the testis or male germ cells and are annotated to male-reproduction-related phenotypes. to gain a deeper insight into the regulatory patterns of meiosis, we selected the highest scoring tf genes for further analysis, based on iregulon’s nes score which corresponds to a low false discovery rate. of these highest scoring tfs, half of them are annotated to meiosis-related functions: etv e2f gata rfx and zfp we included four other well-known meiosis tfs: rarg, mybl etv and tbpl we compared the mrna expression pattern of these tf genes with the expression patterns of their target genes in order to identify putative regulatory relationships, an analysis that is based on the assumption that a relevant tf protein appears more or less contemporaneously with its transcript. in cases where the expression pattern of the tf gene had the opposite expression pattern of its target genes, we infer that the tf acts as a repressor on these targets; but in cases where the expression patterns of the tf gene and its target genes were concordant, we postulate that the tf enhances expression of its targets. for tfs with genes that did not change over our time course, such as the known candidate nrf we could not infer a relationship to target genes and excluded them from further analysis. following these assumptions, as illustrated in figs.  and we suggest that ovol and yy act as repressors of target genes in early leptonema and etv etv and zfp act as repressors on target genes after initiation of meiosis, beginning with the late leptotene and zygotene substages. evidence suggests that zbtb gata and etv also act as repressors on gene targets in pachynema; while rfx appears to be an enhancer of target genes in late pachytene and diplotene substages, but may repress its target genes in preleptonema, although this relationship is less clear. interestingly, these underlying assumptions based on relative expression levels suggest that mybl tbpl and e2f act as both enhancers and repressors of their target genes. mybl and tbpl appear to repress target genes in early meiosis while enhancing target genes in late meiosis, while e2f has the opposite pattern and may repress target genes in late meiosis while enhancing target genes in early leptonema. similar regulatory switching has been shown to result from changes in protein co-factors and post-translational modifications. moreover, genes for some of the tfs described above are also themselves target genes associated with specific substages. by considering these relationships, we inferred candidate regulatory interactions among tfs; for example, zfp and etv suppress mybl gene expression, while mybl enhances expression of rfx ovol yy and tbpl genes.fig. expression of tfs and target genes for each substage. the difference from mean gene expression at dpp of the tf and its target genes at each time point. to illustrate the overall pattern, a smoothed line was fit to the substage-specific gene expression. color of lines represents the substage of target genes. gene expression is average log averaged over time point replicates. expression is plotted for zpf etv mybl and zbtb33 fig. regulation of substage-concordant gene expression during meiotic prophase. analyses suggest that substage-concordant gene expression is regulated both positively and negatively by a network of tfs. tfs are shown by their gene names, and the meiotic substages of predicted tf target genes are shown in the center. black arrows represent a tf acting as an enhancer on its target genes and red bars represent a tf acting as a suppressor on its target genes., lp/d = late pachynema and diplonema, anti-el = genes negatively concordant with early leptonema, anti-ll/z/ep = genes negatively concordant with late leptonema, zygonema, and early pachynema) overall, 41 % of the genes associated with specific meiotic substages are predicted targets of one or more of these eight meiosis-related tfs. this noteworthy observation suggests that this concise regulatory network can account for a substantial portion of the meiotic program of transcription. discussion here we identify a male meiotic germ-cell transcriptome using a novel analysis based on a dense dataset of cytological substage-specific signatures. the high concordance between rna-seq expression data and cytological proportions of isolated germ cells across all samples we analyzed allowed us to develop a novel pmca to identify the substage-specific transcriptomes for meiotic prophase. this computational method does not require the use of facs or sedimentation sorting of cells and can be applied to other complex cell populations for which there are well-defined cytological criteria. cytological deconvolution of rna-seq data from the first wave of spermatogenesis in mammalian males, spermatogenesis is continuous and asynchronous, ensuring daily capacity to deliver sperm backed up by a testis comprised of abundant numbers of germ cells at all stages of spermatogenesis. while advantageous to reproductively active males, this biological imperative has frustrated attempts to isolate stage-purified spermatogenic cells. methods for stratification of cell samples by either sedimentation or flow cytometry enrich specific spermatogenesis cell stages, while the juvenile first wave of spermatogenesis provides a leading edge of differentiating cells against a background of less differentiated cells. in the laboratory mouse, meiosis is initiated at about dpp in a subset of spermatogonia by retinoic acid stimulation of stra; however, initiation of subsequent waves of meiosis make the cellular population increasingly heterogeneous, resulting in a complex histology. the cell-stage frequencies presented here are based on a powerful combinatorial application of antibodies recognizing well-characterized and highly stage-specific marker proteins. we find a continuous and relatively constant background of the earliest spermatogenic cells, namely, spermatogonia and the preleptonema and leptonema, and, not surprisingly, these “background” spermatogenic cells are the ones most refractory for computational assignment of a substage-specific transcriptome. on the other hand, pachytene spermatocytes are not present at the earliest juvenile stages, and appear in a discrete time period; consequently, it was possible to assign a robust gene list to this stage. gene assignments to specific substages were guided by the covariance of transcriptomic data and the cytological data. assignment of expressed transcripts to cellular subpopulations in heterogeneous samples has long been a computational challenge. in this study, we exploited the advantage of paired rna-seq and cytological data to develop the novel pmca approach that isolates changes of gene expression specific to each meiotic substage. maximum covariance analysis was first developed in the meteorological sciences and was popularized in the climatological sciences in the 1990s. more recently, an mca approach has been used in a bioinformatics context to clarify relationships between gene and protein expression. while mca is often an effective tool for detecting common signals in two sets of variables, it can be limited by a tendency to fit spurious patterns when faced with increased sampling variation. current methods employ a parametric smoothing model using principal component regression, which requires a normality assumption or rely on go analysis. we developed a novel pmca method that not only overcomes the spurious pattern identification liability but that does so without the need for any parametric assumptions or reliance on go analysis. our pmca approach is broadly applicable to multi-dimensional data derived from a common set of samples. substage specificity of meiotic gene expression and regulation historically, meiotic prophase substages are characterized by the morphology of chromatin, and correlated genetic mechanisms have been revealed in the past decades. it is known that there is widespread transcription of protein-coding genes in the testis, and indeed, we found that among germ cells alone, more than genes are transcribed, suggesting that a significant portion of testis transcriptional complexity is due to germ cells. using pmca, we successfully identified genes with expression patterns that matched the “cyto-patterns” of stage-specific cell frequencies determined by antibody labeling. some genes were shared among substages; for example, late leptotene and zygotene shared genes with early pachytene and early pachytene shared genes with late pachytene and diplotene. due to the strong temporal signal of contribution of the late pachytene germ cells, which do not appear until dpp and greatly increase in representation by dpp, we identified over late pachytene genes. however, since the proportion of preleptotene and early leptotene cells did not vary greatly over time, we were unable to deconvolve stage-specific transcriptomes for the earliest meiotic substages as confidently as for the later substages. we compared these substage-specific gene lists to those developed in other recent analyses of the germ-cell testis transcriptome. the striking increase in number of genes unique to the pachytene stage was also observed by soumillon, et al., where, in fact, greater numbers of genes were considered expressed than in this study. however, we have used a more stringent cutoff to eliminate from the data transcripts expressed at low levels. when we relaxed this requirement and compared the number of genes expressed with a cutoff of tpm >  we found the number of expressed transcripts to be similar to that reported by soumillion, et al.. we also compared the substage-specific gene lists to the relevant clusters identified by margolin, et al.. in their study, gene lists were temporally clustered from single samples taken at and dpp. while our study was based on multiple replicates at each time point as well as fine-grained cytological analysis, our pmca gene lists substantially overlapped with their derived gene list clusters. in all, of the genes in clusters identified by margolin, et al. were also expressed in our dataset. we also compared our results with those derived from flow-cytometry based methods of sorting by dna content. for both of these studies, we found overall concordance in genes assigned to meiotic substages. these similarities were especially pronounced for late pachytene and diplotene substages due to the large numbers of genes in these sets. however, we also observed a broader trend of alignment between early and late substages. furthermore, our use of multiple markers to cytologically characterize cells paired with our pmca analysis allowed for more precise substage assignment than possible by cell sorting. for instance, our spermatogonia and pre-leptotene sets were the top two strongest overlaps with the “secondary spermatocyte” cell fraction in da cruz, et al., but split these assignments to provide improved substage resolution. similarly, our joint late-leptotene and zygotene genes and early pachytene genes had the two greatest overlaps with the combined “leptotene and zygotene” and “pachytene spermatocytes” genes, and further partition these transcripts into more precise substages. analysis of fallahi et al. data was limited by different experimental platforms, different assay timing, and generally low numbers of uniquely assigned genes, but also revealed significant overlaps in late substages. by assessing functional concordance by similarity in go annotations, we determined that our early leptotene and fallahi et al. pre-leptotene sets share a common enrichment in rna transcription genes. finally, the most prominent divergences between our data and these previous results were in the spermatogonia fraction. we note that cell sorting techniques are susceptible to inclusion of somatic-cells in this fraction, whereas our spermatogonia genes were significantly negatively-concordant with many standard meiosis genes. in sum, these comparisons suggest that our transcript sets encompass and substantially expand these previous findings. in this study, we characterized each meiotic substage using the list of substage-specific genes from pmca analysis. some meiosis genes were associated to particular substages. for example, spata and msh4/ are both highly associated with late leptotene/zygotene substages, and both are required for recombination. however, many canonical meiosis genes were found throughout all meiotic prophase substages; these genes include rad rec and syce this may well reflect a lack of acute transcriptional regulation for these important transcripts. rather than transcription at the precise meiotic substage of use, quality and efficiency of meiosis may be ensured by having transcripts present and available for translation throughout meiotic prophase, with substage transitions regulated post-transcriptionally and/or post-translationally. further, although we confirmed that spermatogonia negatively concordant genes are enriched with meiotic genes, we found several meiotic genes in the gene list concordant with spermatogonia, including fign, which was reported to be required for meiotic recombination, suggesting that some of the meiotic program is set up prior to meiotic entry. not surprisingly, this developmental transcriptome analysis also revealed that many genes expressed during meiosis do not have known functions directly contributing to meiosis. these may instead be part of a parallel program of spermatogenic gene expression. for example, many of the genes expressed in early leptonema are associated with transcription or rna metabolism, as well as with cellular processes such as cell-cell interactions, which are of considerable importance in the biology of the seminiferous tubule. pmca also revealed that late pachytene/diplotene-expressed genes are significantly associated with go terms associated with spermiogenic processes, correlating well with findings of other studies. during post-meiotic stages of nuclear condensation, transcription is globally repressed, underlying the biological rationale for prior transcription to support the active protein synthesis during spermiogenesis. in addition, the late pachytene/diplotene-associated gene lists are significantly enriched with non-protein-coding rnas, including lincrnas and pirna precursors. meiotic regulatory networks analysis of predicted tf targets generated a regulatory network that potentially governs the meiotic and spermatogenic programs of gene expression identified by this study. the protein zfp was central in our network, and is required for embryo development in zebrafish and human. the human ortholog, znf is ubiquitously expressed, and binds to the promoter region of target genes where it is required for formation of chromatin loops by interacting with polii and rad. because rad is differentially expressed after the pachytene stage of meiosis, it is possible that zfp is crucial for spermatogenic processes by regulating the transcription of genes during pachytene/diplotene substage, an idea that should be tested experimentally since our analysis indicates that zfp may target genes strongly associated with genes in these substages. reinforcing the validity of our approach, our regulatory network also involved mybl. our results support previous findings that a subset of late pachynema/diplonema-expressed genes, many involved with the spermatogenic and spermiogenic processes, are associated with mybl. mybl regulates transcription of cell cycle-related genes. germ cells with a homozygous mutation in the mybl gene exhibit defects in meiotic chromosome synapsis. moreover, our results suggest that mybl also is in a network associated with the transcription of genes encoding proteins required for pirna biogenesis, including tdrd tdrd and tdrd consistent with previous reports. going beyond our core regulatory network, we also identified several other tfs, whose target genes include the tdrd family and piwi genes, potentially involved in pirna processing: atf elf elk elk fli nfya, nfyb, nrf pbx rfx rfx sp yy zbtb zfp and zfp this analysis not only revealed a core network for transcriptional regulation of meiotic progression but also suggested that a significant proportion of the genes expressed in the meiotic transcriptome may be controlled by a concise entourage of transcription factors.  this study has untangled in part the complexity and parallel process of spermatogenesis and meiosis by focusing on associating gene expression with highly specific cytological signatures defining meiotic prophase substages. this unique and powerful approach to deconvolving transcriptomes of complex cell populations is applicable for discovery of transcriptional signals in other such complex cell populations or heterogeneous tissues. 
here the juvenile onset of spermatogenesis in mice is analyzed by combining cytological and transcriptomic data in a novel computational analysis that allows decomposition of the transcriptional programs of spermatogonia and meiotic prophase substages. independent principal component analyses of cytological and transcriptomic data yielded similar patterns for both data types, providing strong evidence for substage-specific gene expression signatures. germ cells from testes of individual mice were obtained at two-day intervals from to 18 days post-partum, prepared as surface-spread chromatin and immunolabeled for meiotic stage-specific protein markers. expression of meiosis-specific genes is not substage-limited, suggesting regulation of substage transitions at other levels. a novel permutation-based maximum covariance analysis was developed to map co-expressed transcripts to one or more of the eight meiotic prophase substages, thereby linking distinct molecular programs to cytologically defined cell states. here it revealed not only features of meiotic substage-specific gene expression, but also a network of substage-specific transcription factors and relationships to potential target genes. meiosisspermatogenesismaximum covariance analysismousetranscriptomerna-seqhttp://dx.doi.org/ /100000057national institute of general medical sciences gm99640http://dx.doi.org/ /501100001691japan society for the promotion of science http://dx.doi.org/ /100000054national cancer institute ca034196issue-copyright-statement© the author 2016 this integrated analysis provides a general method for resolving complex cell populations.
 the continuous and non-synchronous nature of postnatal male germ-cell development has impeded stage-specific resolution of molecular events of mammalian meiotic prophase in the testis. here the juvenile onset of spermatogenesis in mice is analyzed by combining cytological and transcriptomic data in a novel computational analysis that allows decomposition of the transcriptional programs of spermatogonia and meiotic prophase substages.  germ cells from testes of individual mice were obtained at two-day intervals from to 18 days post-partum, prepared as surface-spread chromatin and immunolabeled for meiotic stage-specific protein markers. eight stages were discriminated cytologically by combinatorial antibody labeling, and rna-seq was performed on the same samples. independent principal component analyses of cytological and transcriptomic data yielded similar patterns for both data types, providing strong evidence for substage-specific gene expression signatures. a novel permutation-based maximum covariance analysis was developed to map co-expressed transcripts to one or more of the eight meiotic prophase substages, thereby linking distinct molecular programs to cytologically defined cell states. expression of meiosis-specific genes is not substage-limited, suggesting regulation of substage transitions at other levels.  this integrated analysis provides a general method for resolving complex cell populations. here it revealed not only features of meiotic substage-specific gene expression, but also a network of substage-specific transcription factors and relationships to potential target genes. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords meiosisspermatogenesismaximum covariance analysismousetranscriptomerna-seqhttp://dx.doi.org/ /100000057national institute of general medical sciences gm99640http://dx.doi.org/ /501100001691japan society for the promotion of science http://dx.doi.org/ /100000054national cancer institute ca034196issue-copyright-statement© the author 2016 