Induced pluripotent stem cells derived from rabbits exhibit some characteristics of naïve pluripotency

Pluripotent stem cells (PSCs) can exist in two morphologically, molecularly and functionally distinct pluripotent states, designated as the naïve and primed states (Nichols and Smith, 2009). The first rabbit pluripotent stem cell (PSC) lines were produced from New Zealand White rabbit blastocysts in the form of embryonic stem cells [rabbit ESCs (rbESCs)] (Fang et al., 2006; Wang et al., 2007). In their undifferentiated state, rbESCs required both fibroblast growth factor 2 (FGF2) and growth factors of the transforming growth factor (TGF) β family (activin, nodal) for self-renewal (Wang et al., 2008). These results were confirmed in another study, which showed that activin/nodal signaling through Smad2/3 activation was necessary for maintaining the pluripotent status of rbESCs (Honda et al., 2009). These rbESCs did not colonize rabbit embryos after being injected into blastocysts (Honda et al., 2008). Honda et al. also established rabbit induced PSC lines [rabbit iPSCs (rbiPSCs)] from adult liver and stomach cells (Honda et al., 2010). Similar to rbESCs, rbiPSCs are dependent on FGF2 for self-renewal. Therefore, rbESCs and rbiPSCs seemed to exhibit the cardinal feature of primed pluripotency. However, this conclusion should be moderated because none of these studies evaluated these cell lines for other criteria that demarcate naïve and primed pluripotency; moreover, these studies did not compared these cell lines with rabbit inner cell mass (ICM) cells and PSCs from other species. Therefore, the status of rbESCs and rbiPSCs is yet to be thoroughly evaluated. In our study, we compared rbESCs and rbiPSCs at both molecular and functional levels and concluded that (i) rbESCs do not exhibit all the characteristic features of pluripotency and (ii) rbiPSCs exhibit some characteristic features of naïve pluripotency.

Twenty-four ICMs isolated from 52 blastocysts by immunosurgery were plated on feeder cells in rbESC medium supplemented with FGF2. Twelve of these ICMs could form outgrowths, of which 4 produced a population of highly proliferating cells that could be routinely passaged more than 40 times. These rbESC lines – designated rbES-4, rbES-8, rbES-18 and rbES-19 – formed flattened colonies of compact cells (Fig. 1A; supplementary material Fig. S1A). They were passaged every 3 days after collagenase II treatment of the colonies, followed by gentle dissociation into small clumps. All attempts to passage these cell lines after trypsinization into single-cell suspensions resulted in extensive differentiation and apoptosis (data not shown). All the 4 lines were positive for alkaline phosphatase activity (Fig. 1B; supplementary material Fig. S1A) and were strictly dependent on FGF2 for self-renewal, as determined by their altered morphology after FGF2 withdrawal for 24 h (Fig. 1C). Conventional Giemsa staining showed a normal chromosome number in all the selected lines (>50 metaphase spreads analyzed: 42XY) (Fig. 1D; supplementary material Fig. S1A). For rbES-18 line, the cells were karyotyped using G-banding, and normal chromosome complements were observed in the 4 analyzed metaphase spreads (Fig. 1E). Immunolabeling showed that all the 4 lines expressed Oct4 in virtually every cell in the respective populations (Fig. 1F; supplementary material Fig. S3). All the 4 cell lines could undergo in vitro differentiation into embryoid bodies (EBs) (Fig. 1H). Differentiation was accompanied by the rapid loss of Nanog and Oct4 expression and upregulation of the ectodermal marker Nestin and endodermal marker Hnf3β. The mesodermal marker Flk1 was expressed at robust levels in undifferentiated cells (Fig. 1G). All the 4 rbESC lines could induce teratomas after being injected under the kidney capsules in SCID mice. All teratomas contained derivatives of the 3 embryonic germ layers (Fig. 1I,J; supplementary material Fig. S2A). Therefore, these 4 rbESC lines appeared to exhibit the characteristics of PSCs.

Adult rabbit fibroblasts (rbF) were infected twice at a 48-h interval with a mixture of 4 retroviral vectors that expressed human Oct4, Sox2, Klf4 and c-Myc transcription factors. Three days after the second infection, the rbFs were dissociated and replated on feeder cells in an rbESC medium supplemented with 10 ng/ml FGF2. The culture medium was changed every day thereafter until some colonies with compact morphologies appeared (Fig. 2A). One hundred and eleven clones were selected between the 15^th and 29^th day after fibroblast infection. These 111 clones were passaged by trypsin dissociation into single-cell suspensions. Twelve clones expressing alkaline phosphatase and capable of rapid proliferation were selected for further studies. At passage 6, 4 of these clones expressed rabbit Nanog and Oct4 at the same level as that in rbESCs (Fig. 2B). For 3 clones – designated rbiPS-B19, rbiPS-B24 and rbiPS-B29 – silencing of all the 4 transgenes was achieved at passage 25 (Fig. 2C). Line rbiPS-B25 did not show suppressed hc-MYC transgene expression even after 25 passages and was excluded from all subsequent analysis. Withdrawal of FGF2 induced dramatic morphological changes within 48 h, suggesting differentiation (Fig. 2D). A normal chromosome number was observed in 2 lines – rbiPS-B19 and rbiPS-B29 (Fig. 2E; supplementary material Fig. S1B). For the rbiPS-B19 cell line, cells were karyotyped using G-banding, and normal chromosome complements were observed in the 20 analyzed metaphase spreads (Fig. 2F). The rbiPS-B24 cell line showed an abnormal 43XX karyotype. Immunolabelling showed that all the 3 cell lines expressed Oct4 in virtually every cell in the respective populations (Fig. 2G; supplementary material Fig. S3). Differentiation induced by suspension culture resulted in the rapid loss of Oct4 and Nanog expression (Fig. 2H). HNF3β, Flk1 and Nestin were already expressed at robust levels in undifferentiated rbiPS-B19 cells. All the 3 lines induced teratomas containing cells that had differentiated into the 3 germ layers after being injected into SCID mice (Fig. 2I; supplementary material Fig. S2B). Most importantly, all iPSC lines could be routinely passaged every 2 days by trypsin dissociation into single-cell suspensions.

The patterns of stage-specific embryonic antigen (SSEA) expression differed between cell types, and all showed considerable heterogeneity. All the 3 rbiPSC lines showed heterogeneous expression of SSEA1, SSEA4, and TRA-1-60, whereas all the 4 rbESC lines expressed only SSEA1 (supplementary material Fig. S3). Distribution of SSEA1 and SSEA4 expression in various cell populations was quantified using flow cytometry (Fig. 3A). SSEA1 was expressed in 30–60% rbESCs and rbiPSCs, compared with 94% mESCs (Fig. 3B). SSEA4 was expressed in 10–40% rbiPSCs and the pattern was complex with these lines, which comprised 4 cell subpopulations: SSEA1⁺ (32–49%), SSEA4⁺ (5–12%), double-positive (7–26%), and double-negative (29–37%) cells (Fig. 3C). The distributions between these 4 sub-populations differed between lines and were constant over time in culture (data not shown). We determined whether these 4 populations represented 4 distinct cell types or whether they were capable of converting into each other in culture. To address this each subpopulation was sorted by FACS and subsequently cultured for 1 week (secondary cultures) before analysis by flow cytometry. The heterogeneity originally observed in the unsorted cell population was restored from each subpopulation (Fig. 3C). These results indicate that rbiPSCs existed in 4 interchangeable states (SSEA1⁺/SSEA4⁺, SSEA1⁻/SSEA4⁺, SSEA1⁺/SSEA4⁻, and SSEA1⁻/SSEA4⁻), and that the distribution of the 4 populations was cell line specific.

E-cadherin expression is observed in all PSCs regardless of species, whereas N-cadherin is expressed only in mouse Epiblast Stem Cells (EpiSCs) (Tesar et al., 2007). rbESC and rbiPSC lines had different patterns of E- and N-cadherin expression. Immunolabeling showed that all rbESC lines expressed both E- and N-cadherins similar to rhesus ESCs. In contrast, rbiPSC lines expressed only E-cadherin similar to mouse ESCs (Fig. 3D). Each SSEA subpopulation was sorted by FACS, and the expression of E-and N-cadherins was analyzed by qPCR. All subpopulations of rbiPSCs expressed Cdh1 (E-cadherin) at levels 10- to 100-fold higher than those measured in rbESCs. In contrast, the same subpopulations expressed Cdh2 (N-cadherin) at levels 10- to 100-fold lower than those measured in rbESCs (Fig. 3E). Therefore, differential expression of E- and N-cadherins was associated with the ESC vs iPSC status and not with the cell population heterogeneity observed as the SSEA expression.

PSCs show an unusual cell cycle structure that is characterized by a short G1 phase and a high proportion of cells in the S phase. In addition, they do not undergo cell-cycle arrest at the G1 checkpoint in response to DNA damage. Instead, they undergo growth arrest only at the G2 checkpoint (Aladjem et al., 1998; Fluckiger et al., 2006; Filipczyk et al., 2007). Therefore, we investigated whether rabbit PSCs displayed the same cell cycle characteristics. For the 4 rbESC lines, the G1, S, and G2/M phases were observed in 74±2.7%, 20±2.8%, and 6±2.9% cell populations, respectively. For the 3 rbiPSC lines, the G1, S, and G2/M phases were observed in 43±2.4%, 38±3.5%, and 19±1.9% cell populations, respectively (supplementary material Table S3). Thus, the percentage of rbPSCs in the S phase (20–38%) was higher than that in rbFs (7%) but lower than that in mESCs (68%), mEpiSCs (56%), and rhESCs (41%). The cell cycle distributions of the SSEA1⁺ and SSEA4⁺ subpopulations were also analyzed. Of all the analyzed cell lines, the SSEA1⁺ subpopulations showed the lowest G1 and the highest S fractions than the SSEA1⁻ and SSEA4⁺ subpopulations. However, no SSEA1⁺ subpopulation from any line analyzed showed a cell cycle distribution similar to that observed for SSEA1⁺ mESCs (supplementary material Table S3).

To study the cell cycle response of rbPSCs to DNA damage, we analyzed the effects of doxorubicin, a DNA intercalator that induces double-strand breaks (Fig. 4). The G1 fraction of rbFs treated with doxorubicin for 24 h decreased from 80% to 41%, whereas the G2 fraction increased from 12% to 58%, indicating growth arrest in G1 and G2. In contrast, the G1 fraction of doxorubicin-treated mESCs decreased dramatically from 20% to 1%, whereas its G2 fraction increased from 11% to 99%, indicating a lack of G1 checkpoint in mouse PSCs and accumulation of cells at the G2 checkpoint. Similar results were obtained with rhESCs and mEpiSCs for which the G1 and S fractions virtually disappeared after doxorubicin treatment (Fig. 4A). The response of rbPSCs to doxorubicin treatment was more heterogeneous. They all showed moderate decreases in their G1 fractions (rbESCs, from 74% to 70%; rbiPSCs, from 43% to 20%) and significant increase in their G2 fractions (rbESCs, from 6% to 30%; rbiPSCs, from 19% to 80%). We assumed that this might be a reflection of the heterogeneity of rbPSC populations. Therefore, we analyzed the responses of different SSEA fractions to doxorubicin treatment. Only SSEA1⁺/SSEA4⁻ and SSEA1⁺/SSEA4⁺ subpopulations of rbiPSCs showed complete growth arrest in the G2 phase and did not show cell accumulation in the G1 phase; similar results were observed for mESCs, mEpiSCs, and rhESCs (Fig. 4A,B). All other fractions showed variable accumulation of cells in the G1 and G2 phases, as observed previously in the entire population. Of note, the SSEA1⁺ subpopulation of rbESCs still showed accumulation in the G1 and G2 phases after doxorubicin treatment. Thus, these results indicate that only the SSEA1⁺ subpopulation of rbiPSCs showed the cell cycle features of mESCs, mEpiSCs, and rhESCs (i.e. a shortened G1 phase, a high proportion of cells in the S phase, and absence of DNA damage checkpoint in G1). Neither the SSEA1⁺ subpopulation of rbESCs nor the other subpopulations of rbiPSCs showed the cell cycle features of PSCs.

The distal element (DE) of Oct4 enhancer showed robust activity in ICM cells of mouse blastocyst, primordial germ cells, and mESCs. In contrast, it is minimally active in the epiblast of mouse postimplantation embryo and its EpiSC derivatives (Yeom et al., 1996; Tesar et al., 2007). Thus, DE activity can be used to discriminate between the naïve and primed states of pluripotency (Nichols and Smith, 2009). To evaluate the DE activity in rabbit preimplantation embryos, we used the L-SIN-EOS-C(3)-EiP lentiviral vector (designated EOS) that expressed GFP and the puromycin resistance gene under the transcriptional control of the minimal early transposon promoter and trimer of the CR4 enhancer (Hotta et al., 2009). The CR4 enhancer is a 146-bp fragment of the mouse DE element that recruits the Oct4, Sox2, and Sp1 transcription factors and recapitulates the regulation of the DE element in pluripotent and differentiated cell lines (Okumura-Nakanishi et al., 2005). A vector expressing GFP from the ubiquitous CAG promoter [GAE-CAG-eGFP/wpre, (Mangeot et al., 2002)] – designated GAE – was used as a control. Both lentiviral vectors were concentrated and injected under the zona pellucida of 8-cell-stage rabbit embryos. The embryos were then cultured to the mid-blastocyst stage before analysis with confocal microscopy. Fluorescent cells were observed in the trophoblast and ICM, with both control and EOS lentivectors, indicating the activity of the trimer of the CR4 enhancer in rabbit preimplantation embryos (Fig. 5A,B). DE activity in the trophoblast cells was consistent with previous observations, indicating that Oct4 was expressed in the trophoblast of the rabbit blastocyst (Kobolak et al., 2009).

To evaluate the DE activity in rabbit PSCs, rbESC, and rbiPSC lines were infected with EOS and GAE lentivectors and analyzed using flow cytometry after 4 days (Fig. 5C; supplementary material Fig. S4A). All cell lines infected with GAE showed green fluorescence ranging from 12% to 54% of positive cells according to the cell type (threshold set at 2.5×10²). Approximately 60–95% of these GFP-positive cells showed fluorescence levels of ≥10³. In contrast, fluorescence levels in EOS-infected cells varied dramatically between cell types. mESCs had the highest level (76% GFP-positive cells had fluorescence levels of >10³). mEpiSCs showed no expression while rhESCs showed low expression (only 13% GFP-positive cells had fluorescence levels of >10³). This indicates that rbESCs and iPSCs showed intermediate expression levels (rbES-4, 32%; rbES-8, 33%; rbES-19, 30%; rbiPS-B19, 42%; rbiPS-B24, 41%; and rbiPS-B29, 48%; Fig. 5C; supplementary material Fig. S4A). This suggested that the trimer of the CR4 enhancer was more active in rbPSCs than in mouse EpiSCs and primate ESCs.

The GFP expression level in EOS-infected cells was influenced by the infection rate. Therefore, to eliminate the confounding effect of the variations in infection rates between the cell lines, we calculated an expression index (I_E): ⁠. %GFP⁺≥10³ is the percentage of GFP⁺ cells with expression levels of ≥10³; %GFP⁺_T is the percentage of all GFP⁺ cells; and I_R is the percentage of infected cells. Comparison of I_E values between cell lines showed that mESCs had the highest index, whereas rhESCs and EpiSCs had the lowest index. As determined from previous calculations, I_E of all rabbit lines was between that of mESCs and rhESCs, with all the 3 iPSCs lines showing higher I_E indices than the 3 ESC lines (Fig. 5D). From these experiments, we concluded that compared with rbESCs, rbiPSC lines could activate the mOct4 DE element at a higher level.

This conclusion was strongly reinforced by the results of puromycin selection of EOS-infected cells. All ESC and iPSC lines infected with EOS were cultured for 7 days in a medium supplemented with puromycin. Puromycin-resistant cells could be readily selected and expanded from mESCs and the 3 rbiPSC lines (designated rbiPS-B19-EOS, rbiPS-B24-EOS, and rbiPS-B29-EOS) (Fig. 5E). In contrast, no puromycin-resistant cells could be isolated from any of the rbESC lines infected with the EOS lentivector. Therefore, we concluded that only iPSCs could activate the trimer of the CR4 enhancer at a level sufficient to confer resistance to puromycin. This conclusion was supported by the observation that the mean intensity of GFP fluorescence increased during puromycin selection of iPSCs infected with EOS, suggesting that only some rare cells expressing GFP at high levels were resistant to puromycin and expanded during selection. These puromycin-resistant cells had notably higher percentages of SSEA1⁺ and SSEA4⁺ cells than the original cells (i.e. before puromycin selection) (supplementary material Fig. S4B).

Regions CR1 and CR4 were chosen for the methylation study. Bisulfite genomic sequencing showed that both regions were 3–10 times more methylated in rabbit embryonic fibroblasts (REFs) than in rabbit PSCs (Fig. 5F; supplementary material Fig. S5). Regions CR1 and CR4 were less methylated in rbiPSCs than in rbESCs [CR1: 3.1% methylated CpG dinucleotides in rbiPSCs vs 8.6% in rbESCs (P<0.01); CR4: 2.3% methylated CpG dinucleotides in rbiPSCs vs 7.6% in rbESCs (P<0.001)] (Fig. 5G; supplementary material Fig. S5).

To further define the molecular properties of rbESCs and rbiPSCs, we purified the SSEA1⁺/SSEA4⁻, SSEA1⁺/SSEA4⁺, and SSEA1⁻/SSEA4⁻ subpopulations by FACS sorting in 4 rbESC lines and 2 rbiPSC lines and performed a global analysis of their expression profiles. For this, we used a rabbit-specific gene expression microarray containing approximately 13,000 independent genes. For each cell line (rbES-4, rbES-8, rbES-18, and rbES-19; rbiPS-B19 and rbiPS-B29) and SSEA subpopulations, 3 biological replicates were hybridized to the customized rabbit array. Hierarchical clustering of the normalized data showed that all subpopulations of rbESC lines were clustered on one side and all subpopulations of rbiPSC lines were clustered together on the other (Fig. 6A). This indicates that despite the heterogeneity observed in the cell lines, the variability is mainly associated with the cell status (rbESC vs rbiPSC) and not with the SSEA expression status.

To further examine the transcriptome of various rbESC and rbiPSC subpopulations, we examined the expression of rabbit homologs of 22 mouse genes [Blimp1, Cdx2, Cdh1 (E-cadherin), Cdh2 (N-cadherin), Cldn6, Dax1 (NrOb1), Dazl, Essrb, Fbxo15, Fgf4, Gbx2, Klf4, Lefty2, Nanog, Oct4 (Pou5f1), Otx2, Pecam1, Pitx2, Piwil2, Rex1 (Zfp42), Tbx3, and Tcfcp2l1], the mRNA levels of which have been used to evaluate stemness and demarcate the naïve and primed pluripotent states in rodent cells (Tesar et al., 2007; Bao et al., 2009; Hanna et al., 2010; Tang et al., 2010; Wang et al., 2012). Of note, 12 genes (Klf4, Cldn6, Otx2, Nanog, Oct4, Cdx2, Dazl, Rex1, Fbxo15, Pitx2, Tbx3, and Gbx2) were represented on the microarray. The expression levels of these 22 genes were determined by qPCR for the SSEA1⁺/SSEA4⁻, SSEA1⁺/SSEA4⁺, SSEA1⁻/SSEA4⁺, and SSEA1⁻/SSEA4⁻ subpopulations purified from the 4 rbESC and 3 rbiPS lines (supplementary material Table S4; Fig. S6; Fig. 3E). In agreement with the microarray data, correlation clustering of the ΔCt values showed that the rbiPSC subpopulations showed fewer differences between them than with the rbESC subpopulations (Fig. 6B). To determine which genes accounted for the variations between rbESCs and rbiPSCs, PCA was performed from the qPCR data obtained with the SSEA1⁺/SSEA4⁻ and SSEA1⁺/SSEA4⁺ subpopulations (Fig. 6C). The first 2 principal components had 78% of the total variation. The first axis that represented 56% of the total variability distinguished between rbESC and rbiPSC lines in accordance with correlation clustering analysis. The second axis represented 22% of the total variability and accounted for variability between cell lines. A graphical representation of gene dispersion in a trigonometric circle highlighted 2 clusters of genes: the first comprised Cdh2, Gbx2, and Dax1 and the second comprised Otx2, Pitx2, Cldn6, Cdh1, Essrb, Klf4, Dazl, Oct4, Piwil2, Nanog, and Pecam1 (Fig. 6D). Both clusters were clearly anticorrelated over the cell lines, indicating that their relative expression levels were mutually exclusive. Both clusters appeared to be representative because most of the distances between the origin and the projected variables were globally close to the maximum distance (i.e. 1). Thus, these 14 genes explained most of the variability observed between rbESCs and rbiPSCs (first axis in Fig. 6C). The other 8 genes (Fgf4, Cdx2, Rex1, Tcfcp2l1, Lefty2, Fbxo15, Tbx3, and Blimp1) mostly accounted for the variability observed between lines in each group (2^nd axis Fig. 6D).

We next examined which cell types (rbESCs vs. rbiPSCs) and which subpopulations (SSEA1⁺vs. SSEA4⁺vs. SSEA⁻) showed the closest proximity to ICM cells of the rabbit embryo. To this end, the expression of the 22 genes was measured in the ICM of third-day rabbit mid-blastocysts (supplementary material Table S4) and compared to the expression measured in the rbESC and rbiPSC subpopulations. Calculation of the Euclidean distances between the ICM sample and each subpopulations showed that SSEA1⁺ rbiPSCs were the closest to the ICM and the SSEA1⁻ rbiPSCs were the farthest (Fig. 7A). Interestingly, all rbESC subpopulations, irrespective of their SSEA1 status, were more distantly related to the ICM than the SSEA1⁺ rbiPSCs. The relative proximity of the SSEA1⁺ rbiPSCs with the ICM could be further evidenced in a rank histogram of the relative expression level of the 22 genes (Fig. 7B,C). The RΝA level of each gene in ICM cells were normalized to the mRNA level of Tbp, and the resultant ΔCt was ranked from the highest (Oct4: ΔCt = 37.5) to the lowest (Dax1: ΔCt = 0.0013). Most genes that were strongly expressed in rbiPSCs than in rbESCs (Oct4, Cdh1, Essrb, Cldn6, Nanog, Klf4, Otx2, Dazl, Pitx2, Piwil2, Pecam1; P<0.05) were highly expressed in the ICM (ΔCt>1). In contrast, genes that were strongly expressed in rbESCs (Cdh2, Gbx2, Dax1; P<0.05) were expressed at low levels in the ICM (ΔCt<0.1). Taken together, these results indicate that the gene expression profile of SSEA1⁺ rbiPSCs best matched that of the ICM cells of rabbit preimplantation embryo.

To evaluate the capacity of rbESCs to colonize pre-implantation embryos and contribute to the formation of ICM, 5–10 rbES-18 cells that expressed GFP were injected under the zona pellucida of 20 eight-cell-stage rabbit embryos. Of note rbES18 was derived from a GFP transgenic embryo (Al-Gubory and Houdebine, 2006) and GFP was expressed by all cells in this population (Fig. 8A). In a parallel experiment, rbES-18 SSEA1⁺ cells were sorted by FACS and then injected into 67 eight-cell-stage embryos. All embryos were cultured until they reached the mid-blastocyst stage prior to GFP immunostaining and analysis by confocal microscopy (supplementary material Table S5; Fig. 8B). No GFP-expressing cells were observed in the 74 embryos that had successfully developed to the blastocyst stage. In contrast, when 5–10 mES GFP⁺ cells (line CGR8-GFP) were injected into rabbit embryos at the same stage, all the resulting blastocysts (15) exhibited extensive colonization of ICM (Fig. 8C). Thus, rabbit early-cleavage-stage embryos were permissive to colonization by naïve mESCs but not by rbESCs.

We performed a similar experiment with rbiPS-B19 line. First, rbiPS-B19 cells were infected with the GAE lentivector that expressed GFP under the CAG promoter. GFP⁺ cells were sorted by FACS and then injected into 47 eight-cell-stage rabbit embryos. All developed to the blastocyst stage, but only one embryo exhibited GFP⁺ cells in the ICM and in the trophoblast. Similarly, SSEA1⁺ rbiPS-B19-GFP cells were sorted prior to injection into 65 eight-cell-stage embryos. Only one blastocyst had GFP⁺ cells in the ICM after GFP immunostaining (Fig. 8D).

We thought that isolated rbiPSC survival would be very poor after injecting them into eight-cell-stage embryos and that this could explain the low rate of chimaeric blastocysts that was observed. To bypass the dissociation step, we generated aggregation chimaeras with 8-cell-stage embryos. Of note, this technique requires removal of the mucous coat, which prevents embryo implantation into the uterus after transfer to surrogates (Murakami and Imai, 1996). For this goal, rbiPS-B19 cells were infected with ABP-RP-TLCLACZL, a lentiviral vector that expressed the lacZ gene under the control of the EF1α promoter. Viable β-galactosidase⁺ cells were isolated by iterative FACS sorting until a population of rbiPS-B19 cells that expressed β-galactosidase at a robust level could be stably expanded. The rbiPSC-B19-LacZ cells were subsequently infected with the EOS lentivector and selected with puromycin to generate a homogenous population of mOct4 DE-positive cells (not shown). Clumps of 20–25 resulting rbiPSC-B19-LacZ-EOS cells were aggregated in the zona pellucida- and mucous coat-free 8-cell-stage embryos and were subsequently cultured to the mid-blastocyst stage. We obtained 62 blastocysts, of which 10 had β-galactosidase⁺ cells in the ICM. In one case, the ICM was heavily colonized by β-galactosidase⁺ cells (Fig. 8E,F; supplementary material Table S5).

Using the molecular and functional criteria that characterize pluripotency in rodents and primates, we demonstrated that rbESCs and rbiPSCs show some key differences. First, similar to mECSs, rbiPSCs are resistant to trypsinization to single-cell suspensions, a procedure that induces extensive cell death and differentiation when applied to rbESCs, mouse EpiSCs (Tesar et al., 2007), and primate ESCs (Wianny et al., 2008). Second, similar to mESCs, rbiPSCs express elevated levels of E-cadherin, whereas rbESCs express both E- and N-cadherins similar to mouse EpiSCs and primate ESCs (Tesar et al., 2007; Hawkins et al., 2012). Third, the cell cycle of rbESCs had a longer G1 phase than the S and G2 phases and similar to somatic cells, they undergo growth arrest in the G1 phase after DNA damage. This is in sharp contrast to the SSEA1⁺ rbiPSCs as well as to all mouse and primate PSCs examined until date, which have a relatively short G1 phase and lack a DNA damage checkpoint in the G1 phase (Aladjem et al., 1998; Hong and Stambrook, 2004; Fluckiger et al., 2006; Filipczyk et al., 2007; Momcilovic et al., 2010). Fourth, the CR4 element of the mouse Oct4 DE exhibits a more robust activity in rbiPSCs than in rbESCs. High CR4 element activity is observed in mESCs, whereas low activity is observed in mouse EpiSCs (Tesar et al., 2007). Fifth, the Oct4 promoter is less methylated in the rbiPSCs than in the rbESCs. Sixth, rbESCs and rbiPSCs show distinct gene expression profiles. We identified 14 genes whose activity is different between rbiPSCs and rbESCs. Similar to that in mESCs, Essrb, Klf4, Piwil2, Cdh1, Dazl, and Pecam1 are strongly expressed in rbiPSCs; however, they are expressed at low levels in rbESCs. Last, based on the expression of the 22 selected genes, rbiPSC lines are the closest to rabbit ICM, whereas the profiles of all rbESC lines are the farthest. Taken together, it can be concluded that rbiPSCs show many characteristic features of naïve pluripotency as defined in rodents, i.e. resistance to single-cell dissociation, no expression of N-cadherin, robust activity of the distal enhancer of mouse Oct4, and robust expression of ICM-specific markers. However, rbiPSCs do not show all the attributes of naïve pluripotency. For instance, they exploit FGF2 signaling for self-renewal, an attribute of primed pluripotency (Nichols and Smith, 2009). In addition, they do not express all the molecular markers of naïve pluripotency, including Rex1, Tbx3, Gbx2, Fgf4, and Dax1 (Tesar et al., 2007). Furthermore, they show a reduced capacity to colonize the ICM after injection into rabbit blastocysts. Based on these observations, we propose that rbiPSCs self-renew in an intermediate state between naïve and primed pluripotency. Such an intermediate state may resemble the intermediate epiblast stem cells (IESCs) recently described in the mouse (Chang and Li, 2013). Like IESCs, rbiPSCs express Rex1 (naïve pluripotency marker) and Otx2 (primed pluripotency marker) at low levels, and express Pecam1 (naïve pluripotency marker) and Pitx2 (primed pluripotency marker) at high levels. Moreover, both rbiPSCs and IESCs can colonize the ICM after aggregation with morulas. In the case of IESCs, the IESC-derived ICM cells failed to contribute to post-implantation development. In the case of rbiPSCs, we could not address this question because the morula aggregation technique requires removal of the mucin coat, which prevents embryo implantation into the uterus after transfer to surrogates (Murakami and Imai, 1996). However, the growth requirements of IESCs and rbiPSCs seem different. IESCs, which grow on gelatin, respond either to LIF or to activin to self-renew in the undifferentiate state. By contrast rbiPSCs grow on feeders and are strictly dependent on FGF2 signalling. Nevertheless, the possibility remains that rbiPSCs could also be dependent on LIF produced by feeder cells.

We have shown that rbESCs do not show all the attributes of pluripotent stem cells. In particular, they have a DNA damage checkpoint in the G1 phase like somatic cells. It is generally assumed that shortening the G1 phase sustains the pluripotent state by restricting the window of opportunity for differentiation cues (Mummery et al., 1987; Burdon et al., 2002; Sela et al., 2012; Coronado et al., 2013). In bona fide pluripotent stem cells, the lack of a DNA damage checkpoint in G1 might further reduce this window by preventing growth arrest at a critical point for self-renewal. Thus, rbESCs lack some key features of the pluripotent cell cycle, which might explain why they exhibit such a low proliferation rate and high spontaneous differentiation rate than rbiPSCs and primate ESCs.

To conclude, we have reported the first PSC lines in a non-rodent species that show characteristics of naïve pluripotency. These cells represent a first step towards genome engineering and production of germline chimeras in rabbits.

Sexually mature New Zealand White rabbits were purchased from Hycole (Marcoing, France) or HyPharm (Roussay, France). Superovulations of females were induced as described previously (Salvetti et al., 2010). Sixty hours after artificial insemination, the fertilized embryos were flushed from explanted oviducts using Euroflush® (IMV Technologies) and cultured in TCM199 medium (Sigma) supplemented with 10% new-born calf serum (Sigma) at 38°C in 5% CO₂.

Mouse embryonic fibroblasts (MEFs) were prepared from 12.5-day-old embryos from the OF1 strain (Charles River). Rabbit embryonic fibroblasts (REFs) were prepared from 12-day-old embryos from the New Zealand White (Hycole). Rabbit fibroblasts (rbFs) were prepared from the ear skin of a 5-month-old New Zealand White female. MEFs, REFs, rbFs, the 293T cell line (ATCC, CRL 11268), the rhesus monkey ESC (rhESC) line LyonES (Wianny et al., 2008), and mouse ESCs [mESCs, CGR8 line (Mountford et al., 1994)] were cultured as described elsewhere (Sandrin et al., 2002; Wianny et al., 2008; Savatier et al., 1996). EpiSCs were derived from epiblasts of 6.5-day-old mouse embryos of the OF1 strain and cultured as described previously (Brons et al., 2007). rbESCs and rbiPSCs were cultured on mitomycin C-treated MEFs (1.25×10⁴ MEF/cm²) in DMEM/F12 medium supplemented with 20% knockout serum replacement, 1% non-essential amino acids and 1% of a solution of 10,000 U/ml penicillin + 10,000 U/ml streptomycin + 29.2 mg/ml L-glutamine, 1 mM sodium pyruvate and 100 µM 2-mercaptoethanol (Invitrogen). For rbESC and rbiPSCs culture, the medium was supplemented with 13 ng/ml or 10 ng/ml FGF2, respectively. rbiPSCs were routinely dissociated into single-cell suspensions after treatment with 0.05% trypsin–EDTA. rbESCs were dissociated into small clumps of cells after treatment with 1 mg/ml collagenase II (Sigma).

The mucin coat and the zona pellucida of rabbit blastocysts (70 to 72 h post-insemination; mid-blastocyst stage) were mechanically removed after brief exposure (2 min) to 5 mg/ml pronase (Sigma). ICMs were separated from the trophectoderm by immunosurgery after incubation of hatched blastocysts in anti-rabbit whole goat serum (Sigma) at 37°C for 30 min and brief exposure (5 min) to guinea pig complement serum (Sigma). Isolated ICMs were transferred to four-well plates on mitomycin C-treated MEFs at a concentration of 3.5×10⁴ cells/cm² in rbESC medium. After 5 to 6 days, the developing outgrowths were mechanically dissociated into small clumps and transferred onto fresh feeder cells at 2.5×10⁴ cells/cm². All subsequent passages were performed by mechanical dissociation after treating colonies with collagenase II at 1 mg/ml, followed by plating cell clumps on feeder cells at a concentration of 1.25×10⁴ cells/cm².

Plasmids used for producing retroviral and lentiviral vectors were either purchased from Addgene [pMXs-hOCT3/4 (17217), pMXs-hSOX2 (17218), pMXs-hKLF4 (17219), pMXs-hc-MYC (17220) (Takahashi et al., 2007) and pL-SIN-EOS-C(3+)-EiP (21313) (Hotta et al., 2009)] or obtained from Dr Cosset [pTG5349, phCMV-GP and pGAE-CAG-eGFP/wpre (Sandrin et al., 2002; Mangeot et al., 2002)] or Dr Suter [psPax2 and pMD2G (Suter et al., 2006)]. Virus production, concentration and titration were performed as previously described by transfection of 293T cells (Sandrin et al., 2002). Concentrated LacZ-expressing lentiviral particles were purchased from Allele Biotech (ABP-RP-TLCLACZL).

To generate rbiPSCs, 5×10⁵ rbFs were infected twice at a two-day interval with 3×10⁶ of viral particles, which contained equal parts of the four freshly produced hOct4, hSox2, hKlf4 and hMyc retroviruses [multiplicity of infection (m.o.i.) = 8], in the presence of 4 µg/ml polybrene (Sigma). Three days after the second infection, rbFs were trypsinized into single cells and replated at low density (6×10² cells/cm²) onto growth-inactivated MEFs in rbiPSC medium. The medium was changed every other day for 2 to 3 weeks until iPSC-like colonies appeared. To infect ESCs and iPSCs, the cells were plated on Matrigel (Beckton Dickinson) prior to infection with L-SIN-EOS-C(3+)-EiP and GAE-CAG-eGFP/wpre lentivectors at an m.o.i. of 5–10 in appropriate media supplemented with 6 µg/ml of polybrene. Rabbit eight-cell-stage embryos were infected by microinjection of 1 nl L-SIN-EOS-C(3+)-EiP and GAE-CAG-eGFP/wpre lentivectors (2×10⁷ infectious particles/ml) under the mucous coat and zona pellucida.

Alkaline phosphatase and β-galactosidase activities were assessed using a PAL 86R-1KT kit (Sigma) and a GALS staining kit (Sigma), respectively, according to the manufacturer's instructions. To detect β-galactosidase activity in viable cells, LacZ-expressing cells were analysed by flow cytometry using a FluoReporter® lacZ flow cytometry kit (Molecular Probes) according to the manufacturer's instructions.

For immunolabelling, the cells were grown on cover slips pre-treated with gelatin and FBS and coated with feeder cells. The cells were fixed with 2% paraformaldehyde (PFA; Sigma) in phosphate-buffered saline (PBS) at 4°C for 30 min, permeabilized in Tris-buffered saline (TBS, 50 mM Tris pH 7.6, 0.9% NaCl, 0.2% Triton X-100) and washed three times (10 min each) in TBS. Non-specific binding sites were blocked with TBS supplemented with 10% FBS and 1% bovine serum albumin (BSA) for 1 h and subsequently incubated overnight at 4°C with primary antibodies (supplementary material Table S1). After three rinses (10 min each) with TBS, the cells were incubated with fluorochrome-conjugated secondary antibodies (supplementary material Table S1) at room temperature (RT) for 1 h, followed by DNA staining with 0.5 µg/ml Hoechst 33342 for 3 min. Cover slips were mounted with mounting medium M1289 (Sigma), fixed with DPX medium (Cell Path), examined under a conventional fluorescence microscope (Leica DMRE) and analysed with Mercator software (Exploranova). For flow cytometry analysis, the cells were immunolabelled using the same protocols and the same antibodies used for immunostaining on coverslips, except that incubation with antibodies was at RT for 30 min. The cells were analysed using a FACS Canto II cytometer and FACSDiva software (Becton Dickinson). The same protocol was used to detect GFP expression in rabbit blastocysts by immunolabelling, except that embryos were fixed with 2% PFA and permeabilized with 0.5% Triton X-100 in TBS. After DNA staining with 3 µg/ml of propidium iodide in 1% BSA, embryos were incubated sequentially in baths (30 min each) of 10%, 20% and 40% glycerol in TBS and then mounted on slides with mounting medium M1289 (Sigma) containing 0.1% p-phenylene-diamine (w/v) (Sigma). The embryos were observed by confocal microscopy (Leica TCS SP1 microscope with LSM software) at 488 nm and 555 nm (serial slices of 4 mm).

Total RNA was isolated from cell clumps, cell pellets and embryoid bodies (EBs) using an RNeasy mini kit (Qiagen) according to the manufacturer's protocol and reverse-transcribed using 200 U of M-MLV retrotranscriptase (Promega) and random primer mix (MWG). Real-time PCR used the StepOnePlus real-time PCR system (Applied Biosystems) and Fast SBYR® Green Master Mix (Applied Biosystems) according to the manufacturer's instructions. Following 40 amplification cycles, melt-curve analysis was used to verify that only the desired PCR products had been amplified. PCR efficiency for both target and reference genes was determined from the relative quantitative values for calibrator-normalized target gene expression using StepOnePlus Software V2.1 (Applied Biosystems). In all cases, expression of the target genes was normalized to those of rabbit TATA-box binding protein (Tbp) gene. But the same results were obtained with a normalization with the mRNA level of rabbit Gapdh gene. All primers used for RT-PCR and qPCR are shown in supplementary material Table S2.

Chromosome spreads were prepared according to standard procedures and stained with Giemsa solution (Sigma). Thirty to fifty metaphase spreads were counted for each of the rbESC and rbiPSC lines. For karyotyping, G-banding was performed using 0.045% trypsin in Ca²⁺- and Mg²⁺-free Tyrode (Sigma). For some rbESC and rbiPSC line analysed, 4 to 20 banded karyotypes were evaluated for chromosomal rearrangements.

For EB formation, rbESCs and rbiPSCs were dissociated into single cells using 0.05% trypsin, resuspended in rbESC medium without FGF2 and cultured in hanging drops. For teratoma formation, 2–5×10⁶ rbESCs or rbiPSCs were injected under the kidney capsules of 8-week-old severe combined immunodeficient (SCID) mice (CB17/SCID, Charles River). After 4–7 weeks, the mice were euthanized and lesions were surgically removed and fixed in 10% neutral buffered formalin (Curtin Matheson Scientific, Inc.) for paraffin sections, or in 4% PFA for cryosections.

Cells were cultured for 24 h with 50 ng/ml of doxorubicin hydrochloride (Adriamycin®, Sigma). Single-cell suspensions were fixed with 70% ethanol at −20°C, rehydrated with PBS at 4°C for 30 min, rinsed twice with PBS and incubated with 0.1 mg/ml of RNAse A at RT for 30 min. Propidium iodide (30 mg/ml) was added to the cells for 1 min before analysis with a FACS Canto II cytometer (Becton Dickinson). Data were acquired using FACSDiva software (Becton Dickinson). The percentages of G1-, S- and G2-phase cells were determined using ModFitLT software (Verity Software House).

The methylation profile of the rabbit Oct4 promoter (GenBank: AC235550.2) was determined by bisulfite mutagenesis and sequencing as previously described (Borghol et al., 2008). Two regions of the Oct4 promoter were subjected to duplex-nested PCR. Region 1 is situated between positions −307 and +96 (25 CpG), and encompasses Conserved Region 1 (CR1) as defined by Kobolak et al. (Kobolak et al., 2009). Region 2 is situated between positions −2060 and −1626 (17 CpG sites), and encompasses Conserved Region 4 (CR4). Sequence of primers specific for bisulfite-converted DNA are as follows: Region 1, external forward: 5′-GTTTTTTTAGGGAGGGGGTAGAG-3′, external reverse: 5′-AAAACCTTAAAAACTCAACCAAATCC-3′, internal forward: 5′-ATGGGGTGGAAGGGATTTTAG-3′, internal reverse: 5′-AAAATCCACCCAACCTAACTCC-3′; Region 2, external forward: 5′-GTTGGTTGGGTAGGAGTTTAT-3′, external reverse: 5′-TAACCCTATCAAACTTCTAAAAAACT-3′, internal forward: 5′-ATAAGTTAAAGAGTTTTGTTTTTGG-3′, internal reverse: 5′-AACTTCTAAAAAACT AAATAACCTAACTCT-3′. The PCR products were cloned into the pNZY28-A plasmid (Nzyteh; Portugal) and 24 PCR-positive colonies were sequenced for each PCR product (Biofidal, Lyon, France).

A customized rabbit microarray was designed starting from the commercially available Agilent rabbit microarray (G2519F), and including genes expressed during preimplantation development in the rabbit (B.S.-P., L.J. and V.D., unpublished data). Briefly, rabbit embryonic ESTs obtained from four suppressive subtractive cDNA hybridizations (SSH) were included in the custom-made microarray. ESTs expressed during the first cleavage stages, blastocyst stage, and early gastrulation embryos were sequenced (Léandri et al., 2009), assembled into contigs and annotated through homology searches against nucleic and proteic databases. Furthermore, candidate genes sourced from Ensembl Genome Browser (oryCun2.0), and associated with pluripotency or embryonic development, were added to the design. The best quality probes from the commercially available Agilent transcriptome-wide rabbit microarray was added, resulting in a microarray with 62, 976 probes including positive and negative controls. The 60-mer oligosequence probes were designed by eArray (https://earray.chem.agilent.com/earray), and were used for microarrays synthesis in situ using the Agilent SurePrint™ technology (Agilent Technologies, Mississauga, On, Canada) with a 8×60 K format. The custom microarray design of the platform has been submitted to the NCBI GEO (Gene Expression Omnibus). The accession number of the platform is GSE43403.

For hybridization, cyanine-3 (Cy3) labeled cRNAs were prepared from 0.2 µg RNA using the One-Color Microarray-Based Gene Expression Analysis Low Input Quick Amp Labeling kit (Agilent) according to the manufacturer's instructions. Dye incorporation and cRNA yield were checked with the NanoDrop ND-1000 Spectrophotometer. 0.6 µg of Cy3-labelled cRNA (specific activity >6.0 pmol Cy3/ug cRNA) was fragmented at 60°C for 30 minutes following the manufacturer's instructions. Hybridizations were performed for 17 hours at 65°C in a rotating Agilent hybridization oven. After hybridization, microarrays were washed 1 minute at room temperature with GE Wash Buffer 1 (Agilent) and 1 minute with 37°C GE Wash buffer 2 (Agilent), then dried immediately. Slides were scanned immediately after washing on the Agilent DNA Microarray Scanner using one color scan setting for 8×60k array slides (Scan Area 61×21.6 mm). The scanned images were analyzed with Feature Extraction Software (Agilent, 7.10.3.1). Features flagged in Feature Extraction as Feature Non-uniform outliers were excluded. Normalization procedures were performed using R statistical software (R: A language and environment for statistical computing (R Foundation for Statistical Computing), Vienna, Austria; 2012, R-Core-Team; http://www.r-project.org) and data were normalized using intra-array median substraction and log2 transformation.

Principal component analysis (PCA) was performed on qRT-PCR data (centered and reduced) using R statistical software (FactoMineR: Multivariate exploratory: data analysis and data mining with R; R package version 1.20, 2012, Husson et al.; http://cran.r-project.org/web/packages/FactoMineR/index.html). Other statistical analyses were done using the test of Kruskal–Wallis. Normalized transcriptome data were analyzed by hierarchical clustering using the Pearson correlation coefficient as a measure of distance.

For microinjection, ESCs and iPSCs were dissociated into single-cell suspensions. Five to 10 cells were microinjected under the mucous coat and zona pellucida of eight-cell-stage rabbit embryos. The embryos injected with GFP-expressing cells were cultured in vitro until they reached the mid-blastocyst stage, fixed in 2% PFA and immunostained with anti-GFP antibody.

For morula aggregation, fertilized eggs were harvested 20 hours after insemination, and cultured in RDH medium (Invitrogen), as described elsewhere (Jin et al., 2000). When they reached the eight-cell stage, embryos were treated with pronase to eliminate the zona pellucida and placed with clumps of 20–25 PSCs in microwells to induce aggregations, as described previously (Wood et al., 1993). Successful aggregates were cultured for 48 hours until they reached the blastocyst stage and subsequently analysed for β-galactosidase activity.

We thank Nicolas Gadot and Jean-Yves Scoazec (ANIPATH, SFR Santé Lyon-Est, Lyon, France) for teratoma histology, Thibault Andrieu and Sébastien Dussurgey (SFR BioSciences Gerland–Lyon Sud) for technical assistance with flow cytometry and cell sorting, the INRA SIGENAE team (http://www.sigenae.org) for their help in designing the customized rabbit array, and Jérome Lecardonnel for technical assistance in array hybridization. This work was supported by research grants from “Agence Nationale de la Recherche” (ANR) (project PLURABBIT no. PCS-09-GEM-08), European Cooperation in Science and Technology (COST) Action (Rabbit Genome Biology (RGB)-Net no. TD1101), the HyPharm Company (PhD Grant INRA/HyPharm no. 04012010) and the “région Rhône-Alpes” (ADR CIBLE 2010 project no. R10065CC). This work was performed within the framework of (i) the LABEX DEVWECAN (ANR-10-LABX-0061) and the LABEX CORTEX (ANR-11-LABX-0042) of Université de Lyon (“Investissements d’Avenir” no. ANR-11-IDEX-0007), (ii) the Institute CESAME (“Investissements d’Avenir” no. ANR-10-IBHU-003) and (iii) the “Infrastructure Nationale en Biologie et Santé” INGESTEM (“Investissements d’Avenir” no. ANR-11-INBS-0009), operated by the ANR.