Germline development in rat revealed by visualization and deletion of Prdm14

The germline is a unique cell lineage capable of transmitting genetic information to the next generation. In mammals, all germline cells, including sperm and eggs, originate from primordial germ cells (PGCs) specified from the pluripotent epiblast in pre-gastrulating embryos (Saitou and Yamaji, 2012; Tang et al., 2016). Studying PGC specification and development is important for understanding the fundamentals of cell fate determination as well as the mechanisms underlying disease and infertility.

Germline development in mammals has been extensively studied in the mouse model system. Loss-of-function studies using knockout mice has led to the identification of essential inductive signals for PGC fate: WNT and BMP (Aramaki et al., 2013; Lawson et al., 1999; Ohinata et al., 2009; Ying and Zhao, 2001), and the downstream core transcriptional regulators Prdm1 (also called Blimp1), Prdm14 and Tfap2c (also called Ap2γ) (Ohinata et al., 2005; Weber et al., 2010; Yamaji et al., 2008). Furthermore, remarkable progress in reconstituting germline development in vitro from mouse pluripotent stem cells (Hayashi et al., 2012, 2011; Hikabe et al., 2016) has enabled comprehensive analysis (Kurimoto et al., 2015; Shirane et al., 2016) and large-scale genetic screening (Hackett et al., 2018), which was difficult previously due to the limited number of nascent PGCs in vivo.

Here, we examine germline development in the rat model system. Rats are widely used as laboratory animals. Many genetic and physiological features of rats are much closer to humans than those of mice are (Aitman et al., 2008). The study of rats, as the closest laboratory animal to mice, is essential to determine conserved mechanisms of germline development among rodents. Nevertheless, only limited information on the rat germline is available, likely due to the lack of appropriate genetically modified lines. Histological analysis to identify rat PGCs by morphology (Beaumont and Mandl, 1963), alkaline phosphatase activity (Kemper and Peters, 1987) and immunostaining of some representative markers (Encinas et al., 2012; Leitch et al., 2010; Northrup et al., 2011) have provided a rough timeline for rat germline development. In a recent report, to visualize gonadal to adult germ cells in vivo, Ddx4-EGFP transgenic rats were generated (Gassei et al., 2017). However, owing to the absence of Ddx4 expression in earlier PGCs, the specification and migration stages of rat PGCs have not yet been documented. Thus, by using novel genetically modified rats, we aimed to collect fundamental information, such as developmental timing and gene expression changes during germline development in rats, to compare with the well-studied mouse model.

Of the unique and important regulators for mammalian germline, we focused on Prdm14, a PR domain-containing transcriptional regulator. In mice, Prdm14 is specifically expressed in pluripotent cells and in nascent PGCs during specification of the germline (Yamaji et al., 2008). Loss of Prdm14 in mice causes lack of germ cells by mid-gestation and subsequent infertility (Yamaji et al., 2008). Mechanistic analyses using mouse pluripotent stem cells have demonstrated the important role of Prdm14 in PGC fate induction, epigenetic reprogramming, and acquisition of naïve pluripotency (Grabole et al., 2013; Leitch et al., 2013; Mallol et al., 2019; Nakaki et al., 2013; Okashita et al., 2016; Yamaji et al., 2013).

Here, we genetically modified the rat Prdm14 locus. We introduce a fluorescent reporter, histone H2B fused to Venus (H2BVenus), in the Prdm14 locus to monitor and purify developing rat PGCs. We also generate Prdm14-deficient rats, which reveal a conserved function of Prdm14 among rodents.

To visualize the expression pattern of Prdm14 in rats, we replaced exons 1-4 of the endogenous Prdm14 with the H2BVenus fluorescent gene in rat embryonic stem cells (ESCs) via homologous recombination (Fig. 1A, Fig. S1A,B). All Prdm14-H2BVenus heterozygous knock-in rat ESC lines homogeneously expressed H2BVenus (Fig. S1C,D), akin to mouse ESCs cultured under defined conditions (Yamaji et al., 2013). We generated chimeric rats by injecting a knock-in rat ESC line into wild-type rat blastocysts. Crossing the chimeras with wild-type rats yielded heterozygous (Prdm14^+/HV) rats obtained through germline transmission of the ESC line (Fig. S1E). We then analyzed the expression pattern of Prdm14-H2B Venus in the embryos and fetuses that were obtained by crossing the Prdm14^+/HV rat with wild-type rats.

During preimplantation development, two-cell-stage embryos started to express H2BVenus in the nucleus and subsequently maintained expression until the blastocyst stage (Fig. 1B). In mice, Prdm14 is heterogeneously expressed in four-cell-stage embryos, and might be involved in promoting lineage allocation towards the inner cell mass (Burton et al., 2013). In agreement, some blastomeres in four-cell-stage rat embryos showed higher H2BVenus expression than others (Fig. 1B, arrowheads in 4cell). Moreover, the inner cell mass containing pluripotent cells in the blastocyst at embryonic day (E) 4.5 showed a brighter H2BVenus signal than the surrounding trophectoderm layer (Fig. 1B, blastocyst). Although further functional validation is necessary, the role of Prdm14 in preimplantation development might be well-conserved within rodents.

After implantation, although fluorescence was almost undetectable prior to gastrulation (∼E7.5), a cluster of cells at the posterior epiblast adjacent to the extra-embryonic ectoderm in gastrulating E8.5 embryos faintly, but specifically, showed expression of H2BVenus (Fig. 1B, arrowhead in E8.5). The number of H2BVenus⁺ cells increased at E9.5. The expression pattern of H2BVenus in whole-mount immunostained E9.5 embryos closely correlated with that of Tfap2c, an essential PGC specifier (Weber et al., 2010), suggesting that H2BVenus specifically marks PGCs (Fig. 1C). H2BVenus⁺ PGCs expressing Tfap2c migrated along the hindgut at E10.5-11.5 and reached the gonadal region at around E12.5 (Fig. 1B,C). In gonads, rat PGCs expressing Tfap2c as well as Ddx4 progressively expanded by E15.5, which is when the gonads show obvious morphological difference between sexes (Fig. 1B,C, Fig. S2). At later stages, the expression level of H2BVenus and the number of H2BVenus⁺ cells in the gonadal germ cells gradually declined until birth, and germ cells in testes or ovaries of postnatal day 8 rats did not show any expression (Fig. S2). Collectively, inserting the H2BVenus gene into the rat Prdm14 locus reveals its unique expression pattern, in pluripotent cells during preimplantation development and in PGCs during post-implantation development.

We next examined the timing of PGC specification in rats based on H2BVenus fluorescence in peri-gastrulating embryos (∼E7.75-8.75). We used morphological features to stage embryos (Fig. S3A) and whole-mount immunostaining for Oct3/4 (also known as Pou5f1) (epiblast) and brachyury (T) (primitive streak/mesoderm) to define developmental timing accurately.

At the pre-streak stage (E7.75), the epiblast uniformly expressed Oct3/4 but T was not detectable (Fig. 2A, E7.75). However, at the early streak stage (E8.0), the posterior epiblast started to express T, indicating onset of gastrulation (E8.0; Fig. 2A). Notably, a few T⁺ cells showed faint expression of Prdm14-H2BVenus (Fig. 2A, arrowheads in insets of E8.0 embryo), suggesting that nascent PGCs in rat might be specified from T⁺ cells at this stage. The number of H2BVenus⁺ nascent PGCs among T⁺ cells at the posterior epiblast gradually increased with elongation of the primitive streak at E8.25, defined as the mid-streak stage (Fig. 2A, E8.25). Then, at late primitive streak no bud (LSOB) stage (E8.5), the nascent PGCs formed an obvious cluster (Fig. 2A, E8.5), as shown in Fig. 1B. At E8.75, when the allantois bud appears and the amnion closes, the size of the cluster increased whereas T expression declined, indicating suppression of somatic genes (Fig. 2A,B, E8.75). The presence of H2BVenus⁺ cells distant from the cluster might reflect dynamic movement of rat PGCs, as has been observed in mice (McDole et al., 2018).

We next examined expression of the tripartite core PGC specifiers Prdm1, Tfap2c and Prdm14 (H2BVenus), which are sufficient to induce PGC fate in the mouse epiblast (Magnúsdóttir et al., 2013; Nakaki et al., 2013). At E8.5 (LSOB), cells at the posterior epiblast expressed either Prdm14-H2BVenus and/or Prdm1 but not Tfap2c, whereas the extra-embryonic ectoderm strongly expressed Tfap2c (Fig. 2C, Prdm14^+/HV). The presence of either Prdm1 or Prdm14 single-positive cells suggests that T might induce these transcription factors independently, as in mice (Aramaki et al., 2013). At early to late head fold stage (E-LHF: E9.5), PGCs expressed all three regulators concomitantly, suggesting that these factors have established the network for germ cell fate (Fig. 2D, Prdm14^+/HV). Tfap2c is a downstream target of Prdm1 and Prdm14 (Cheetham et al., 2018; Kurimoto et al., 2015). This chronological order of key transcription factor expression is reflected in our immunostaining with the reporter in rat.

We also checked for expression of the pluripotency genes Oct3/4 and Sox2, both of which are important for the maintenance of PGCs in mice (Campolo et al., 2013; Okamura et al., 2008). At E8.5 (LSOB), Oct3/4 expression was slightly reduced in H2BVenus⁺ nascent PGCs as well as surrounding T⁺ cells, whereas Sox2 expression is absent and restricted to the anterior epiblast (E8.5; Fig. 2A, Fig. S3B). At E9.5 (E-LHF), H2BVenus⁺ PGCs co-express both Oct3/4 and Sox2, as well as Nanog, another core pluripotency gene (Fig. S6C), suggesting that, like in mouse, rat PGCs also upregulate a pluripotency-associated gene program (Prdm14^+/HV; Fig. 2D, E9.5; Fig. S3B).

Because Prdm14 is the earliest marker for nascent PGCs, we next analyzed its function in rats. Prdm14-deficient rats were obtained by crossing Prdm14^+/HV with Prdm14^+/mut, another mutant line with a deletion in exon 1-4 like the Prdm14^+/HV (Fig. S1F,G). The use of two lines allows us to compare the expression level of H2BVenus precisely in both heterozygous and homozygous mutants, as both have a single copy of the H2BVenus gene.

In Prdm14^mut/HV embryos at E8.5 (LSOB), we detected similar numbers of Prdm14-H2BVenus⁺ nascent PGCs as were observed in Prdm14^+/HV embryos (Figs 2C and 3B). The nascent PGCs expressed Prdm1 even in the absence of Prdm14 in Prdm14^mut/HV embryos (Fig. 2C) concordant with Prdm1 expression being independent of Prdm14. In contrast, at E9.5 (EHF-LHF), H2BVenus expression in Prdm14^mut/HV embryos was not upregulated as it is in Prdm14^+/HV embryos (Fig. 2D). More importantly, the faint H2BVenus⁺ cells did not express Prdm1, Tfap2c or Sox2 at either the protein or mRNA level (Prdm14^mut/HV; Figs 2D and 3A). The number of H2BVenus⁺ cells was unchanged from E8.5, which is in contrast to heterozygous or wild-type embryos in which the PGC numbers increased by 3-fold within 1 day (Fig. 3B). Thus, loss of Prdm14 leads to failure in activation and/or maintenance of key regulators that are essential for the PGC fate as well as subsequent re-acquisition of pluripotency genes. Notably, even in the absence of Prdm14, mouse PGCs delay but retain the expression of several germ cell markers as well as pluripotency genes at the transcript level (Grabole et al., 2013; Shirane et al., 2016; Yamaji et al., 2008). Although further careful analysis matching the reporter and developmental stage between mouse and rat are necessary, Prdm14 might regulate PGC specification more stringently in rats than in mice.

As a result of the collapse in the pre-migratory PGCs, gonads of Prdm14-deficient embryos at E15.5 completely lost germ cells, whereas the supporting somatic cells were entirely normal (Fig. 3C, Fig. S4). Although the testes and ovaries at postnatal day 1 appeared similar between the knockout and wild type, they were smaller in young (4 weeks) and adult (13 weeks or 18 weeks) rats due to the lack of germ cells (Fig. 3D, Fig. S4). Therefore, loss of Prdm14 in rats causes infertility both in males and females like in Prdm14-deficient mice (Yamaji et al., 2008).

To obtain further information on rat PGCs at each developmental stage, we purified PGCs at various stages and used RNA sequencing (RNA-seq) to investigate changes to the transcriptome during development. We collected H2BVenus⁺ cells at E9.5 (early PGCs), E10.5 and E11.5 (migratory PGCs), and E12.5 and E15.5 (gonadal PGCs) (Fig. S5B). In addition to the PGC samples, we also isolated the post-implantation epiblast at E7.25 and E7.75 as reference, because the epiblast differentiates into all germ layers including the germline (Fig. S5A).

A heat map of the correlation coefficients shows high correlation between each replicate and cell type (Fig. S5C). Hierarchical clustering of 8966 detectable genes resulted in two distinct clusters consisting of the pluripotent epiblast and developing PGCs (Fig. S5D). In the principal component analysis (PCA), PC1 reflects dynamic changes in mitochondrial gene expression during PGC development (Fig. S5E,G) as shown for mice and human (Cree et al., 2008; Floros et al., 2018), and the PC2/PC3 scatterplot nicely reflects a temporal trajectory from epiblast to PGCs (Fig. 4A, Fig. S5E). From the PCA, we extracted positively and negatively expressed gene sets correlated with the PC loading (Fig. S5F) and performed K means clustering analysis using the gene sets (K=9) (Fig. 4B). Cluster 1/2/3 genes, upregulated in the epiblast samples, are enriched in gene ontogeny (GO) terms related to morphological and metabolic pathways of post-implantation epiblast (Fig. 4B, right). Importantly, de novo DNA methyl-transferase genes (Dnmt3a, Dnmt3b and Dnmt3l) were highly enriched in Cluster 1/2/3, and dramatically downregulated in the PGCs (Fig. 4B,C). In addition, expression of Uhrf1, which is essential for maintenance of DNA methylation, was also reduced in PGCs at E9.5 compared with the epiblast (Fig. 4C). Subsequently, methylation-sensitive germ cell-related genes, such as Ddx4, Dazl and Mael, were progressively upregulated (Fig. S6C), suggesting that rat PGCs at E9.5 initiate extensive epigenetic reprogramming. Cluster 4, transiently upregulated in early PGCs at E9.5, contained ‘embryo development’ genes including some mesoderm markers (Cdx2, Cdx4, Hoxd1 and Tbx3), which are suppressed by the germ cell program at later stages (Fig. 4B, Fig. S6C). Cluster 5/6 contained modestly (Cluster 5) or abundantly (Cluster 6) expressed key germ cell-related genes (Prdm1, Prdm14, Tfap2c, Nanos3, Dnd1, Nr5a2 and Zfp296), all transcription factors crucial for specification and/or maintenance of PGCs in mice and the other species (Fig. 4B,C) (Bleckwehl and Rada-Iglesias, 2019; Sybirna et al., 2019b). Cluster 7/8 contained highly expressed gametogenesis and meiosis genes associated with late PGCs (Fig. 4B). Taken together, transcriptome analysis of purified rat PGCs uncovers gene sets that mark key events and transition points during PGC specification and subsequent development.

Finally, to reveal conserved/species-specific gene expression patterns in PGCs, we performed a cross-species comparison among mouse, rat and human early/late gonadal PGCs. Hierarchical clustering of 980 homologous genes revealed separate clustering depending on the species (Fig. S6A). A heat map of correlation coefficients shows that each developmental stage of PGCs is highly correlated particularly within rodents (Fig. S6B). A 1- to 2-day lag in gestation period between mouse and rat reflects the difference in PGC development timeline.The expression pattern of key PGC genes exhibits some conservation across the mammals but there are significant differences between rodents and human (Fig. 4D). Sox17, a key transcriptional regulator of PGCs in non-rodent mammals (Irie et al., 2015; Kobayashi et al., 2017), is almost absent in mouse and rat gonadal PGCs, whereas it is temporarily upregulated in early PGCs in rats (Fig. S5F) and mice (Yabuta et al., 2006). Furthermore, expression patterns of pluripotency genes such as Sox2, Klf4 and Tfcp2l1 also show a similar divergence, confirming that the PGC program is highly conserved between rat and mice.

In summary, by genetically modifying a core germ cell specifier, Prdm14, in rats, we obtained a faithful reporter line to monitor and purify developing rat PGCs and a knockout line to reveal its crucial function in PGC fate. Overall, the dynamics and transcriptional regulatory network of PGCs are well conserved within the rodents, albeit delayed in rats due to a difference in gestation periods. However, some differences might impact the regulation of pluripotency and germline development. For instance, although lack of Prdm14 causes infertility in both mouse and rat, the influence on PGC specification might be more severe in rat than in mouse. As naïve pluripotent stem cells in rats can be maintained only under very stringent conditions (Buehr et al., 2008; Buehr et al., 2003; Li et al., 2008), rat PGCs might also be sensitive to changes in the transcriptional regulatory network. A recent work demonstrated that PRDM14 is crucial for human PGC specification and activates distinct transcription programs in mouse and human (Sybirna et al., 2019a preprint), playing context-dependent roles. Similar analysis of Prdm14 function in rats would extend our understanding of PGC networks across species.

Further investigation in rats using the models we have generated, and comparative studies across species will reveal the conserved principles of germline development, which will be informative for the establishment of a widely applicable in vitro system for reconstituting the rat germline cycle.

Crlj:WI rats (RGD ID: 2312504) were purchased from Charles River Laboratories Japan (Kanagawa, Japan). C57BL/6N mice were purchased from Clea Japan (Tokyo, Japan). All experiments were performed in accordance with the animal care and use committee guidelines of the National Institute for Physiological Sciences and Tokyo University of Agriculture.

For construction of targeting vector, 5′ and 3′ homology arms amplified from genomic DNA of rat ESCs (WDB/Nips-ES1/Nips, RGD ID: 10054010) (Hirabayashi et al., 2012), H2BVenus-IRES-Puro^r amplified from prRosa26-H2BVenus-IRES-Puro^r vector (Kobayashi et al., 2012), and MC1-DTA cassette amplified from a MC1-DTA vector were assembled and inserted into pBluescript KS(+) (Stratagene) with an in-fusion HD cloning kit (Takara Bio). For efficient gene targeting, we used the CRISPR/Cas9 system. Guide RNAs targeting around the start codon sequence of rat Prdm14 genes (A: 5′-caccgGGTCTCGCGAGAAGGCGATA-3′, 5′-aaacTATCGCCTTCTCGCGAGACCc-3′; B: 5′-caccgACATTCGTCAAGTAGACGAG-3′, 5′-aaacCTCGTCTACTTGACGAATGTc-3′) were cloned into eSpCas9(1.1) (Addgene plasmid #71814).

For gene targeting, undifferentiated rat ESCs (WDB/Nips-ES1/Nips, RGD ID: 10054010) were maintained on mitomycin-C-treated mouse embryonic fibroblasts (MEFs) in N2B27 medium containing 1 μM PD0325901 (Axon, Groeningen, The Netherlands), 3 μM CHIR99021 (Axon) and 1000 U/ml of rat leukemia inhibitory factor (Millipore). To introduce the vectors into rat ESCs, we carried out reverse transfection using Lipofectamine 2000 (Thermo Scientific) as described before (Kobayashi et al., 2017). In brief, 2×10e5 rat ESCs were suspended in 100-200 µl of Opti-MEM (Thermo Scientific) containing the targeting vector, CRISPR/Cas9 vector, and Lipofectamine 2000 (Thermo Scientific), and were incubated for 5 min at room temperature. Then, rat ESCs were seeded onto four-drug resistant (DR4) MEFs and 48 h later, 0.5 µg/ml puromycin (Sigma-Aldrich) was added to the culture medium for selection. After the selection, puromycin-resistant and H2BVenus⁺ rat ESC colonies were picked. Correct targeting was confirmed by PCR using primers shown in Fig. S1A and Table S1 and a representative result is shown in Fig. S1B.

For making chimeric rats and obtaining germline transmission, we performed blastocyst injection of targeted rat ESCs. Rat blastocysts were collected in mR1ECM medium from oviduct and uterus of rats 4.5 days post-coitum (dpc). For micro-manipulation, rat ESCs were trypsinized and suspended in rat ESC medium. A piezo-driven micro-manipulator (Prime Tech, Tokyo, Japan) was used to drill zona pellucida and trophectoderm under the microscope and eight rat ESCs were introduced into blastocyst cavities near the inner cell mass. After the injection, the blastocysts were transferred into the uterus of pseudopregnant recipient Wistar rats (3.5 dpc).

For efficient and rapid generation of Prdm14 mutant rats, we performed electroporation to introduce CRISPR/Cas9 system into fertilized rat embryos. Ribonucleic complexes composed of 100 ng/µl crRNA (the target is the same with A and B in Prdm14-H2BVenus rats), 200 ng/µl tracrRNA and 200 ng/µl Cas9 protein (all products were purchased from IDT) were prepared according to the manufacturer's protocol. Rat pronuclear stage embryos were collected in mKRB medium from the oviducts of rats 0.5 days dpc. Electroporation was carried out using NEPA21 (Nepagene) according to the manufacturer's protocol. After electroporation, the embryos underwent follow-up culture for 1 day in mKRB medium and were transferred into the oviduct of pseudopregnant recipient Wistar rats (0.5 dpc). Large deletion mutants were judged by PCR and subsequent Sanger sequencing. Of 13 mutant lines, we used two mutant lines (Mutant-3 and -10) with identical phenotypes for the following experiments (Fig. S1F,G).

For isolating rat epiblast (Prdm14^+/HV) at E7.25 and E7.75, the ectoplacental cone and extra-embryonic ectoderm was mechanically cut off. Then, the visceral endoderm surrounding the epiblast was removed by enzymatic treatment using 1.0 mg/ml collagenase type IV (Worthington Biochemical Corporation), 1.0 mg/ml dispase (FUJIFILM Wako, Osaka, Japan) and 0.3 mg/ml hyaluronidase (Sigma-Aldrich) (Rivera-Pérez et al., 2007) for 3-5 min at 37°C followed by repeated pipetting with a glass capillary tube. For purifying rat PGCs from Prdm14^+/HV rat embryos, the embryonic tissues indicated in Fig. S5A,B were dissected and trypsinized using 0.25% trypsin/EDTA for 5 min at 37°C. H2BVenus⁺ rat PGCs were collected with the SH800 cell sorter (Sony). The genotype of the embryos was determined by H2BVenus fluorescence or PCR using ectoplacental cone and/or extra-embryonic ectoderm.

Embryos were fixed in 4% paraformaldehyde for 10-30 min at room temperature (RT) or from 4 h to overnight at 4°C depending on the stages. For cryosections, the fixed embryos were embedded in OCT compound. Each sample was incubated with primary antibodies for 1-2 h at RT or overnight at RT or 4°C, and then with fluorescent-conjugated secondary antibodies with DAPI (FUJIFILM Wako) for 1 h at RT or overnight at 4°C. Antibodies used are listed on Table S2. For embryos, the samples were cleared with CUBIC1/2 or CUBIC2 alone according to a published protocol (Kubota et al., 2017). The embryos and cryosections were observed under a confocal laser-scanning microscope (Nikon A1Rsi).

cDNAs of single rat PGCs were amplified following a previously described method with a minor modification as follows (Kurimoto et al., 2006). For cell collection, Prdm14^+/HV and/or Prdm14^mut/HV embryos were obtained by crossing Prdm14^+/HV and Prdm14^+/mut rats. By using the SH800 cell sorter, H2BVenus⁺ cells in embryos at early head fold stage (E9.5) were collected in a 50 µl droplet of PBS containing 0.5% polyvinyl alcohol on the cover of a 96-well tissue culture plate (IWAKI). Then, the cells were picked up into 0.5 ml PCR tubes (Life Technologies) containing 4.5 µl of cell-lysis buffer with ERCC RNA spike-in Mix (Thermo Fisher) (Nakamura et al., 2015), using a mouth pipette as described (Kurimoto and Saitou, 2011), and were snap frozen in liquid nitrogen. The frozen cells were stored in dry ice or at −80°C until use. Genotypes of the embryos were judged by PCR using the ectoplacental cone as described above. We decided to compare gene expression levels of cells from littermate embryos. Just before use, the frozen cells were thawed on ice, and single cell cDNA amplification was performed as described (Kurimoto et al., 2006), with a minor modification. In brief, for the 20-cycle PCR amplification step after the second-strand synthesis, we added 1 µl of 0.375 µg/µl V1(dT)24 primer to each reaction tube (final concentration of 20 ng/µl), instead of 19 µl of a mixture containing ExTaq buffer, dNTPs, V1(dT)24 primer and ExTaq, which was used in the original protocol (Kurimoto et al., 2006). Expression levels of selected genes were analyzed by real-time PCR using the primers shown in Table S3. Samples with extremely low/high levels of spike RNAs (ERCC) and/or housekeeping genes [Arbp (Rplp0), Gapdh] were eliminated from the analysis as outliers.

Total RNA was extracted using the PicoPure RNA Isolation Kit (Thermo Fisher Scientific) according to the manufacturer's protocols. Two hundred picograms of total RNA were reverse transcribed using SMART-Seq v4 Ultra Low Input RNA Kit (Takara Bio). The synthesized cDNA was amplified with 15 cycles of PCR. Each of the cDNA libraries were indexed and prepared for Illumina sequencing using KAPA Hyper Prep Kit (Kapa Biosystems) with 13 cycles of PCR. The RNA-seq libraries were quantified by qPCR using KAPA Library Quantification Kit (Kapa Biosystems). All libraries were mixed and subjected to single-end 75 bp sequencing on NextSeq 500 system (Illumina) using High Output Kit v2. Basecalls were performed using NextSeq 500/550 RTA software (v. 2.4.11). FASTQ files were generated using bcl2fastq (v. 2.18.0.12).

Mouse male and female PGCs at E10.5 were isolated and subjected to RNA-sequencing preparation. Briefly, heterozygous Pou5f1-ΔPE-GFP mouse embryos (C57BL/6N background) were recovered at E10.5, and the sex was distinguished based on PCR amplification of Zfy gene (Kobayashi et al., 2013). GFP-positive PGCs were sorted by flow cytometry. Purified 500 pg to 1 ng of total PGC RNA were used to generate cDNA libraries with 12 cycles of PCR amplification using the SMARTer v4 Ultra Low Input RNA Kit. Pre-amplified cDNA was fragmented into 200-bp fragments using an S2 sonicator (Covaris) and then used to construct sequencing libraries with 15 cycles of PCR using a NEBNext Ultra DNA Library Prep Kit (New England BioLabs). After qPCR quantification using the KAPA Library Quantification Kit, all libraries were mixed and subjected to single-end 50 bp sequencing on HiSeq 2500 system (Illumina). Basecalls were performed using HiSeq RTA software (v. 1.17.20). FASTQ files were generated using bcl2fastq (v. 2.17.1.14). Raw sequence data of E13.5 male and female PGCs was downloaded from DDBJ (accession: DRA003597) (Sakashita et al., 2015).

For processing the rat and mouse RNA-seq data, the first three bases were trimmed out from raw reads using FASTX Toolkit 0.0.13.2. Trimmed reads were further quality-trimmed using the Cutadapt v. 1.14 program with the following configurations: trimming the primer ‘Clontech SMART CDS Primer II A: AAGCAGTGGTATCAACGCAGAGTAC’; minimal base quality: 20; minimal read length after trimming: 35 nt. RNA-seq reads were aligned to the rat (rn6) or mouse (mm10) genome assembly using TopHat2 (v. 2.1.0) with bowtie2 (v. 2.3.4.1). Cufflinks version 2.2.1 program together with Ensembl gene annotation (rn6, 32545 genes; mm10, 46073 genes) were used to assemble transcripts and estimate their abundances. The read densities [reads per kilobase of transcript per million reads (RPKM)] were converted to the log₂ (RPMK+1) scale (log₂ expression level), and averaged among the biological replicates.

For PCA and GO analysis, rat genes showing maximum log₂ (RPKM+1) values >3 in at least one replicate were selected (8966 genes). Hierarchical clustering was performed based on Ward's method using ‘hclust’ function of R package. The PCA was performed using the ‘prcomp’ function of the R package without scaling. Significant differentially expressed genes (DEGs) were defined as genes showing a more than 2 s.d. radius of scaled PC2 and PC3 values (1174 genes). Significant genes were grouped into K-means clusters (K=9) using the ‘kmeans’ function of R package. GO analysis was performed using DAVID (v. 6.8).

For cross-species comparison between human, raw sequence data of human PGCs was downloaded from NCBI: SRP057098 (PRJNA280935). Raw reads were quality- and adaptor-trimmed using the Trim_galore v. 0.6.4 program with the following configurations: minimal base quality: 25; minimal read length after trimming: 35 nt. RNA-seq reads were aligned to the hg38 genome assembly using TopHat2 (v. 2.1.0) with bowtie2 (v. 2.3.4.1). Cufflinks version 2.2.1 program together with Ensembl gene annotation (hg38, 60594 genes) were used to assemble transcripts and estimate their abundances. The RPKM values were converted to the log₂ (RPMK+1) scale, and averaged across biological replicates. For comparison, one-to-one orthologous genes in the rat, mouse and human were listed using BioMart (Ensembl) with Z-score values. Orthologs of significant DEGs were selected (980 genes), and hierarchical clustering was performed using R package. Correlation coefficients between two variables (mouse/rat/human samples) were calculated using Excel.

We thank members of the Hirabayashi lab, particularly Keiko Yamauchi, Megumi Hashimoto and Naoko Niizeki for help with animals, and Minako Ohnishi for secretarial support. We thank Prof. Yayoi Obata, Prof. Tomohiro Kono and Dr Akihiko Sakashita of Tokyo University of Agriculture for their help in preparing mouse PGCs and valuable discussion regarding data analysis. We also thank Dr Roopsha Sengupta for editing and providing critical input to the manuscript. We thank the Spectrography and Bioimaging Facility, NIBB Core Research Facilities for technical support.

The peer review history is available online at https://dev.biologists.org/lookup/doi/10.1242/dev.183798.reviewer-comments.pdf