High-quality single-cell transcriptomics from ovarian histological sections during folliculogenesis

Cell culture

mESCs were cultured in a medium containing 2iLIF (GSK inhibitor [CHIR99021; 3 μM], MEK inhibitor [PD0325901; 0.4 μM], and LIF [1 U/μl] [2iLIF]) as described previously (Hayashi et al, 2011).

Sectioning of frozen cell blocks

To make frozen cell blocks, the mESCs were pelleted and washed with PBS. The mESCs were then suspended in 1 × PBS containing 10% PVA (Sigma-Aldrich), embedded in a cryomold, and snap-frozen in liquid nitrogen (N₂). Subsequently, the frozen cell blocks were sectioned at a thickness of 15 μm, using a CM1860UV cryostat (Leica), and were mounted on MembraneSlide 1.0 PEN slides (Carl Zeiss). The sectioned cell blocks were dried overnight at room temperature and fixed with 80% ethanol. The fixed sections were then stained with 1% cresyl violet acetate (MP Biomedicals) dissolved in 80% ethanol, washed with 100% ethanol, and dried for more than 1 h. For the histological images, contrast and brightness were optimized linearly with Photoshop for the ease of visual inspection. These procedures are summarized in Fig 1A.

Sectioning, fixation, and staining of mouse ovaries

Mouse ovaries were excised from 8-wk-old female mice and snap-frozen in liquid N₂. Then, the frozen ovaries were sunk in 10% PVA in cryomolds and again frozen in liquid N₂. The embedded ovaries were sectioned using a CM1860UV cryostat at a thickness of 15 μm, dried at room temperature for more than 1 h, fixed and stained with 1% cresyl violet acetate dissolved in 50% isopropanol, washed with 100% isopropanol, and dried for more than 1 h. Formalin fixation was performed with 10% formalin neutral buffer solution (Wako) at room temperature for 24 h, and the fixed sections were stained with 1% cresyl violet acetate in 50% isopropanol and dried for more than 1 h.

Laser capture microdissection

The mESCs on the frozen sections were isolated using a laser microdissection system PALM MB4 (Zeiss) with a ×20 objective lens. The dissection laser setting, cutting speed, and laser pulse catapulting (LPC) were set as 46%, 50%, and 54%, respectively. The dissected cells were collected into the lysis buffer in the caps of single, flat-top 200-μl PCR tubes (Greiner Bio-One).

In mouse ovarian sections, we analyzed follicle morphology using PALM MB4 microscopy. Specifically, we selected normally growing primary follicles and secondary-to-early antral follicles in sections containing nuclei of oocytes for subsequent transcriptome analysis. To isolate oocytes and granulosa cells, we identified single granulosa cells and oocytes under microscopic inspection with a ×60 objective lens. To avoid collecting granulosa cells overlapped within a single section, we carefully inspected the morphology of the nuclei and drew lines to guide the cutting laser, encircling the single nuclei using the “Joint Cut” mode. Given that the thickness of the sections was 15 μm, which is similar to the size of a single granulosa cell, we consider it unlikely that two or more entire cells were isolated simultaneously with this procedure. The parameters used for isolating oocytes and granulosa cells were as follows: cutting energy at 44%, LPC energy at 30%, and speed at 40. The dissected cells were collected into the lysis buffer in the caps of single, flat-top 200-μl PCR tubes. These LCM procedures are summarized in Fig 1A and the results of LCM on ovarian sections are shown in Fig S6A and B.

DRaqL-SC3-seq cDNA amplification

Cells in alcohol-fixed sections were isolated with LCM in 6.4 μl of cell lysis buffer (0.8 μl of GeneAmp 10xPCR Buffer II [Thermo Fisher Scientific], 0.48 μl of 25 mM MgCl₂ [Thermo Fisher Scientific], 0.8 μl of 5% sodium deoxycholate [SDc] [Nacalai Tesque], 0.4 μl of 100 mM dithiothreitol [DTT] [Thermo Fisher Scientific], 0.64 μl of 40 U/μl RNaseOUT [Thermo Fisher Scientific], 0.08 μl of 40 U/μl porcine liver RNase inhibitor [Takara Bio], 0.16 μl of 2.5 mM dNTP [Takara Bio], 0.16 μl of 1:500,000 diluted ERCC RNA Spike-In Mix 1 [Thermo Fisher Scientific], 0.16 μl of 10 ng/μl V1[dT]₂₄ primer [Hokkaido System Science], and 2.72 μl of deionized distilled water [DDW, Gibco]) in the flat caps of 0.2-ml PCR tubes (Greiner Bio-One), and spun down into the tubes by brief centrifugation. Cells were then lysed at 70°C for 6 min, followed by the addition of 2.8 μl quenching buffer (0.7 μl of Triton X-100 [Nacalai Tesque], 0.7 μl of 2% BSA [Takara Bio]) and incubation at 70°C for 90 s to quench the denaturing effect of SDc.

The cDNA synthesis and amplification were performed as described previously (Kurimoto et al, 2006) with minor modifications. In brief, the first-strand cDNA was synthesized by adding 0.5 μl of reverse transcriptase mix (0.2 μl of 200 U/μl SuperScript III [Thermo Fisher Scientific], 0.033 μl of 40 U/μl porcine liver RNase inhibitor, 0.067 μl of 1–10 mg/ml T4 gene 32 product [Roche], and 0.2 μl of DDW [Gibco]) to the cell lysate (9.2 μl) and incubating at 50°C for 5 min followed by heat inactivation at 70°C for 10 min. The excess primers were degraded by adding 1 μl of exonuclease I mixture (1×exonuclease I buffer and 0.5 U exonuclease I [Takara Bio]) and incubating at 37°C for 30 min followed by heat inactivation at 80°C for 25 min. Then, the first-strand cDNA product was poly dA-tailed by adding 3 μl TdT mixture (0.6 μl of 10× PCR Buffer II, 0.36 μl of 25 mM MgCl₂, 0.18 μl of 100 mM dATP, 0.3 μl of 15 U/μl recombinant terminal deoxynucleotidyl transferase [Thermo Fisher Scientific], 0.3 μl of 2 U/μl RNaseH [Thermo Fisher Scientific], and 1.26 μl of DDW) and incubating at 37°C for 1 min followed by heat inactivation at 70°C for 10 min. The dA-tailed product (13.7 μl) was divided into four tubes, and second-strand cDNA was synthesized by adding 9.5 μl V3-PCR mixture (0.95 μl of 10× ExTaq Buffer [Takara Bio], 0.95 μl of 2.5 mM dNTP, 0.1875 μl of 1 μg/μl V3[dT]₂₄ primer, 0.095 μl of ExTaq Hot Start version [ExTaqHS, Takara Bio], and 7.3175 μl of DDW) and applying thermal cycling program 1. Then, cDNA amplification was performed by adding 9.5 μl of V1-PCR mixture (0.95 μl of 10× ExTaq Buffer [Takara Bio], 0.95 μl of 2.5 mM dNTP, 0.1875 μl of 1 μg/μl V1[dT]₂₄ primer, 0.095 μl of ExTaqHS, and 7.3175 μl of DDW) and applying thermal cycling program 2. For oocytes, we used a PCR program with 18 cycles. The amplified cDNA (∼89.7 μl) was purified with a 0.6× volume of AxyPrep MAG PCR clean-up reagent (Axygen) according to the manufacturer’s instructions. All oligonucleotides and thermal cycling programs used in this study are listed in Table S5. DRaqL and adapted cDNA amplification are schematically represented in Fig 1A and B.

DRaqL-SMART-Seq v4 cDNA amplification

A Smart-seq v4 3′ DE Kit (Takara Bio) was adapted to DRaqL as follows: cells in alcohol-fixed sections were isolated with LCM in 6.4 μl of cell lysis buffer (0.25 μl of 40 U/μl Takara RNase Inhibitor, 0.4 μl of 100 mM DTT, 0.16 μl of 1:500,000 diluted ERCC RNA Spike-In Mix 1, 0.4 μl of 5% SDc, and 5.19 μl of nuclease-free water) in the caps of PCR tubes, and spun down into the tubes by brief centrifugation. Cells were then lysed at 72°C for 3 min, followed by the addition of 2.8 μl quenching buffer (0.7 μl of Triton X-100, 0.7 μl of 2% BSA, and 1.4 μl of 5× SuperScript II buffer [Thermo Fisher Scientific]) to reduce the denaturing effect of SDc. Then, 2.8 μl of RT oligo (1 μl of 12 μM Oligo dT In-line Primer and 1.8 μl of nuclease-free water: kit components) was added, followed by incubation at 72°C for 90 s. Reverse transcription was performed by adding 7.5 μl of Master Mix (4 μl of 5× ultra-low first-strand buffer, 1 μl of 48 μM SMART-Seq V4 oligonucleotide, 0.5 μl of 40 U/μl RNase inhibitor, and 2 μl of SMARTScribe Reverse Transcriptase: kit component) supplemented with 0.25 μl of SuperScript II and 0.25 μl of SuperScript III, and by incubating at 42°C for 90 min and 70°C for 10 min. Then, cDNA amplification was performed by adding 30 μl of kit component PCR Mix (25 μl of 2× SeqAmp PCR Buffer, 1 μl of 12 μM Blocked PCR Primer II A, 1 μl of SeqAmp DNA polymerase, and 3 μl of nuclease-free water) and applying the thermal cycling program for SMART-seq v4. The amplified cDNA (50 μl) was purified with a 0.8× volume of AxyPrep MAG PCR clean-up reagent according to the manufacturer’s instructions.

DRaqL-smart-seq2 cDNA amplification

Cells in alcohol-fixed sections were isolated with LCM in 6.4 μl of cell lysis buffer (0.6 μl of 5× SuperScript II buffer, 0.1 μl of 100 μM Oligo dT VN, 0.8 μl of dNTP mix [25 mM each], 0.25 μl of 40 U/μl recombinant RNase inhibitor [Takara Bio], 0.5 μl of 100 mM DTT, 0.06 μl of 1 M MgCl₂ [Sigma-Aldrich], 2 μl of 5 M Betaine [Sigma-Aldrich], 0.16 μl of 1:500,000 diluted ERCC RNA Spike-In Mix 1, 0.4 μl of 5% SDc, and 1.53 μl of DDW) in the caps of PCR tubes, and spun down into the tubes by brief centrifugation. The isolated cells were lysed at 72°C for 6 min, followed by addition of 2.8 μl quenching buffer (0.7 μl of Triton X-100, 0.7 μl of 2% BSA, and 1.4 μl of 5× SuperScript II buffer [Thermo Fisher Scientific]). Then, 1.6 μl of a template-switching mixture (0.1 μl of 100 μM N-template-switching oligo, 0.25 μl of SuperScript II, 0.25 μl of SuperScript III, 0.2 μl of recombinant RNase inhibitor, and 0.8 μl of DDW) was added, followed by the cycling RT program for DRaqL-Smart-seq2. Then, cDNA amplification was performed by adding 15 μl of Seq amp PCR mixture (12.5 μl of 2x SeqAmp buffer [Takara Bio], 0.05 μl of 100 μM N-IS PCR primer, 0.5 μl of SeqAmp DNA polymerase [Takara Bio], and 1.95 μl of DDW) and applying the thermal cycling program for SeqAmp. For oocytes, we used a PCR program of 18 cycles. The amplified cDNA (25.8 μl) was purified with a 0.8× volume of AxyPrep MAG PCR clean-up reagent according to the manufacturer’s instructions.

DRaqL-Protease-smart-seq2 cDNA amplification

Cells in formalin-fixed sections were isolated with LCM in 6.4 μl of cell lysis buffer (0.6 μl of 5× SuperScript II buffer, 0.1 μl of 100 μM Oligo dT VN, 0.8 μl of dNTP mix (25 mM each), 0.25 μl of 40 U/μl recombinant RNase inhibitor, 0.5 μl of 100 mM DTT, 0.06 μl of 1 M MgCl₂, 2 μl of 5 M Betaine, 0.4 μl of 5% SDc, 0.32 μl of 900 mAU/ml QIAGEN Protease, and 1.37 μl of DDW) in the caps of PCR tubes, and spun down into the tubes by brief centrifugation. The isolated cells were lysed by protease digestion at 50°C for 10 min followed by heat inactivation at 80°C for 15 min. The denaturing effect of SDc was quenched by addition of 0.8 μl of 1:2,500,000 diluted ERCC RNA Spike-In Mix 1 and 2.8 μl of quenching buffer (0.7 μl of Triton X-100, 0.7 μl of 2% BSA, and 1.4 μl of 5× SuperScript II buffer [Thermo Fisher Scientific]), followed by incubation at 72°C for 90 s. Then, 0.8 μl of a template-switching mixture (0.1 μl of 100 μM N-template-switching oligo, 0.25 μl of SuperScript II, 0.25 μl of SuperScript III, and 0.2 μl of recombinant RNase inhibitor) was added, followed by the cycling RT program for DRaqL-Smart-seq2. Then, cDNA amplification was performed by adding 15 μl of the SeqAmp PCR mixture, and applying the thermal cycling program for SeqAmp. The amplified cDNA (25.8 μl) was purified with a 0.8× volume of AxyPrep MAG PCR clean-up reagent according to the manufacturer’s instructions.

SC3-seq cDNA amplification with Triton X-100

To evaluate cell lysis efficiency with DRaqL, the cell lysis step of DRaqL-SC3-seq was replaced with cell lysis with Triton X-100 only (SC3-seq cDNA amplification [Triton X-100]) as follows. Cells were isolated with LCM in 6.4 μl of cell lysis buffer (0.8 μl of GeneAmp 10xPCR Buffer II, 0.48 μl of 25 mM MgCl₂, 0.8 μl of 5% Triton X-100, 0.4 μl of 100 mM DTT, 0.64 μl of 40 U/μl RNaseOUT, 0.08 μl of 40 U/μl porcine liver RNase inhibitor, 0.16 μl of 2.5 mM dNTP, 0.16 μl 1:500,000 diluted ERCC RNA Spike-In Mix 1, 0.16 μl of 10 ng/μl V1[dT]₂₄ primer, 2.72 μl of DDW) in the flat caps of PCR tubes, and spun down into the tubes by brief centrifugation. Cells were then lysed at 70°C for 6 min, followed by the addition of 2.8 μl quenching buffer (0.7 μl of Triton X-100, 0.7 μl of 2% BSA) and incubation at 70°C for 90 s to quench the denaturing effect of SDc. The following steps were performed as in DRaqL-SC3-seq.

SC3-seq cDNA amplification with NP-40

As a positive control of the cDNA amplification of 10-pg total RNA purified from 2i-LIF mESCs, we performed the original protocol of the SC3-seq cDNA amplification with the same procedure as described in Kurimoto et al (2006), with minor modifications. Briefly, 0.5 μl of 25 pg/μl total RNA were added to 4 μl of cell lysis buffer (0.8 μl of GeneAmp 10xPCR Buffer II, 0.48 μl of 25 mM MgCl₂, 0.8 μl of 5% NP-40, 0.4 μl of 100 mM DTT, 0.64 μl of 40 U/μl RNaseOUT, 0.08 μl of 40 U/μl Porcine Liver RNase Inhibitor, 0.16 μl of 2.5 mM dNTP, 0.16 μl 1:500,000 diluted ERCC RNA Spike-In Mix 1, 0.16 μl of 10 ng/μl V1[dT]₂₄ primer, 2.72 μl of DDW). The mixture was incubated at 70°C for 6 min, followed by the addition of 2.8 μl DDW, and then was subjected to first-strand cDNA synthesis. The following steps were performed as in DRaqL-SC3-seq.

SMART-Seq v4 3′DE Kit cDNA amplification

The cDNA amplification using SMART-Seq v4 3′DE Kit was performed according the manufacturer’s instruction (Takara Bio), except for the step of cell capture. 20 μl of 10× reaction buffer was prepared (19 μl of 10×Lysis Buffer, 1 μl of RNase Inhibitor), of which 1 μl was mixed with 10.5 μl DDW, resulting in 11.5 μl of cell lysis buffer. Then, 6.4 μl of cell lysis buffer was added the flat caps of PCR tube, and cells were isolated with LCM into the reaction buffer, followed by the brief centrifugation. Then, another 5.1 μl of cell lysis buffer was added to the tubes (11.5 μl in total). The following steps were performed exactly according to the manufacturer’s instruction. Briefly, 1 μl of 12 μM Oligo dT In-line Primer was added to the mixture and incubated at 72°C for 3 min. Then, reverse transcription was performed by adding 7.5 μl of Master Mix (4 μl of 5× ultra-low first-strand buffer, 1 μl of 48 μM SMART-Seq V4 Oligonucleotide, 0.5 μl of 40 U/μl RNase Inhibitor, 2 μl of SMARTScribe Reverse Transcriptase) and by incubating at 42°C for 90 min and 70°C for 10 min. Then, cDNA amplification was performed by adding 30 μl of kit component PCR Mix (25 μl of 2× SeqAmp PCR Buffer, 1 μl of 12 μM Blocked PCR Primer II A, 1 μl of SeqAmp DNA Polymerase, and 3 μl of nuclease-free water) and applying the thermal cycling program for SMART-seq v4.

Smart-seq2 cDNA amplification

Smart-seq2 cDNA amplification was performed as described (Nichterwitz et al, 2016), except for the volume of cell lysis solution. Briefly, cells were isolated with LCM in 6.4 μl cell lysis solution (1 μl of 10 μM Oligo dT VN, 1 μl of dNTPs mix [10 mM each], 0.21 μl of 40 U/μl Recombinant RNase inhibitor, 0.17 μl of 5% Triton X-100, 0.16 μl of 1:500,000 diluted ERCC RNA Spike-In Mix 1, and 3.86 μl of DDW) in the caps of PCR tubes, and spun down into the tubes by brief centrifugation. The isolated cells were lysed at 72°C for 3 min, and were put on ice immediately after the cell lysis. Then, 5.45 μl of the reverse transcription mixture (2 μl of 5× SuperScript II buffer, 0.5 μl of 100 mM DTT, 2 μl of 5 M Betaine, 0.1 μl of 1 M MgCl₂, 0.25 μl of 40 U/μl Recombinant RNase inhibitor, 0.1 μl of 100 μM template-switching oligo, and 0.5 μl of SuperScript II) was added to the tube, followed by cycling RT program for Smart-seq2. Then, cDNA amplification was performed by adding 15 μl of KAPA PCR mixture (12.5 μl of 2× KAPA HiFi HotStart ReadyMix, 0.2 μl of 10 μM IS PCR primer, 2.3 μl of DDW), and applying the thermal cycling program for HiFi.

Protease-Smart-seq2 cDNA amplification

Cells in formalin-fixed sections were isolated with LCM in 6.4 μl of cell lysis buffer (0.6 μl of 5× SuperScript II buffer, 0.1 μl of 100 μM Oligo dT VN, 0.8 μl of dNTP mix [25 mM each], 0.25 μl of 40 U/μl recombinant RNase inhibitor, 0.5 μl of 100 mM DTT, 0.06 μl of 1 M MgCl₂, 2 μl of 5 M Betaine, 0.4 μl of 5% Triton X-100, 0.32 μl of 900 mAU/ml QIAGEN Protease, and 1.37 μl of DDW) in the caps of PCR tubes, and spun down into the tubes by brief centrifugation. The isolated cells were lysed by protease digestion at 50°C for 10 min and then heat inactivated at 80°C for 15 min, followed by the addition of 0.8 μl of 1:2,500,000 diluted ERCC RNA Spike-In Mix 1 and incubation at 72°C for 90 s. Then, 3.6 μl of a template-switching mixture (1.4 μl of 5× SuperScript II buffer, 0.1 μl of 100 μM N-template-switching oligo, 0.25 μl of SuperScript II, 0.25 μl of SuperScript III, 0.2 μl of recombinant RNase inhibitor, and 1.4 μl DDW) was added, followed by the cycling RT program for DRaqL-Smart-seq2. Then, cDNA amplification was performed by adding 15 μl of the SeqAmp PCR mixture, and applying the thermal cycling program for SeqAmp. The amplified cDNA (25.8 μl) was purified with a 0.8× volume of AxyPrep MAG PCR clean-up reagent according to the manufacturer’s instructions.

Pooled and bulk RNA-seq of granulosa cells using DRaqL-SC3-seq

For pooled cell analysis, about 10 granulosa cells were isolated with LCM from the regions neighboring oocytes and non-neighboring regions in early antral follicles in the alcohol-fixed, mouse-frozen ovarian sections (n = 7 each). Cell lysis and cDNA amplification were performed using the DRaqL-SC3-seq method with a reduced number of PCR cycles (18 cycles).

For bulk RNA-seq of granulosa, the whole granulosa of individual sectioned follicles were isolated using LCM without considering the histological affiliation of individual cells. Early antral follicles were identified in the sections that contained oocytes (n = 4), for which oocytes were first removed by LCM and then whole granulosa were collected. The developmental stages of follicles that did not contain oocytes within the investigated sections were not identified. Cell lysis and cDNA amplification were performed using the DRaqL-SC3-seq method with a reduced number of PCR cycles (18 cycles).

Quantification of amplified cDNA

The amount of the amplified cDNAs was measured with Qubit 4 Fluorometer using 1X dsDNA High-Sensitivity Kit (Thermo Fisher Scientific) according to the manufacturer’s instruction. Briefly, 1 μl cDNAs and 10 μl standard solutions were added to 199 μl and 190 μl of 1X working solution, respectively, in a 0.5 ml thin-wall, clear PCR tube (Axygen). The mixtures were vortex and incubated for 2 min at room temperature.

Gene expression analysis with real-time PCR

Real-time PCR was performed using a CFX384 real-time PCR system (Bio-Rad) with Power SYBR Green PCR Master Mix (Thermo Fisher Scientific), according to the manufacturer’s instructions. The primers used for real-time PCR were as previously described (Nakamura et al, 2015) and are listed in Table S5.

Quality control of amplified cDNAs

The success rate of the cDNA amplification from single cells isolated with LCM was evaluated as follows. The expression level of a housekeeping gene, Arbp, was quantified with real-time PCR. Then, the cDNAs for which the Smirnov–Grubbs test for outliers yielded a P-value < 0.01 were considered to have been unsuccessfully amplified and were removed from the subsequent analyses (Fig S1). The size of the amplified cDNA was analyzed using Bioanalyzer with High-Sensitivity DNA kit (Agilent Technologies).

Library preparation for 3′-sequencing with DRaqL-SC3-seq

Preparation of the libraries for SC3-seq was performed as previously described (Nakamura et al, 2015, 2017) with minor modifications. Before the library preparation, ∼5 ng of amplified single-cell cDNAs were pre-amplified with thermal cycling program 3 in 15 μl V1-NV3 pre-amplification mixture (containing 1.5 μl of 10× ExTaq Buffer, 1.2 μl of 2.5 mM dNTP, 0.15 μl of 1 μg/μl N-V3[dT]₂₄, 0.15 μl of 1 μg/μl V1[dT]₂₄, and 0.075 μl of ExTaqHS). The pre-amplified cDNAs were subjected to a primer–dimer removal by 6–9 cycles of DNA purification with a 0.6× volume of AxyPrep reagent, and elution in 20 μl DDW. The size-selected cDNAs were then fragmented with Picoruptor (Diagenode) with 10 cycles of a 30-s sonication and a 30-s interval at 4°C, followed by end-polishing in 25 μl End polish mix (20 μl of the sonicated DNA, 2.5 μl of 10× NEB Next End Repair Reaction Buffer [NEB], 0.085 μl of T4 polynucleotide kinase [NEB], 0.085 μl of T4 DNA polymerase [NEB], and 2.33 μl of DDW), and incubating at 20°C for 30 min. The fragmented, end-polished cDNAs were then subjected to size selection, with 0.7× volume AxyPrep reagent for cutting off the larger molecular size, and the supernatant was recovered by supplementing with an additional 0.2× volume of AxyPrep reagent (0.9× volume total) and eluted in an 8 μl DDW. Next, internal adapter extension was performed in a 10-μl internal-adapter mixture (1 μl of 10× ExTaq Buffer, 0.94 μl of 2.5 mM dNTP, 0.67 μl of 10 μM tRd2SPV1[dT]₂₀, 0.067 μl of ExTaqHS, and 7.3 μl of the eluted DNA) with thermal cycling program 4. Then, the T-adapter ligation was performed by adding 6.63 μl of a ligation mixture (3.3 μl of 5× NEB Next Quick Ligation Reaction Buffer [NEB], 0.2 μl of 10 μM tRd1SP adapter [double-stranded DNA formed by annealing tRd1SPTs and tRd1SPTas oligos], 0.33 μl of 2,000 U/μl T4 DNA ligase [NEB], and 2.8 μl of DDW), and incubating at 20°C for 15 min and 72°C for 20 min, and the ligation product was purified with a 0.8× volume of AxyPrep reagent. Finally, an additional PCR was performed in a 20-μl S5-N7 PCR mixture (2 μl of 10×ExTaq Buffer, 1.6 μl of 2.5 mM dNTP, 2 μl of 10 μM S5-index primer, 2 μl of 10 μM N7-index primer, 0.1 μl of ExTaqHS, 2.3 μl of DDW, and 10 μl of the eluted DNA) with thermal cycling program 5, followed by DNA purification with a 0.9× volume of AxyPrep reagent. Then, libraries with different index sequences were mixed and subjected to sequencing from the Read1 primer (75-bp) using a NextSeq 500/550 High-Output Kit v2.5 (75 cycles) (Illumina).

Library preparation for Y-shaped adapter sequencing for cDNAs prepared with the DRaqL-SC3-seq cDNA amplification method

Single-cell cDNAs were pre-amplified in a 15 μl V1–V3 pre-amplification mixture (1.5 μl of 10×ExTaq Buffer, 1.2 μl of 2.5 mM dNTP, 0.15 μl of 1 μg/μl V1[dT]₂₄, 0.15 μl of 1 μg/μl V3[dT]₂₄, 0.075 μl of ExTaqHS, 7 μl of DDW, and 5 μl of the cDNA) with thermal cycling program 3. Pre-amplified cDNAs were then subjected to the primer–dimer removal, sonication, end-polishing, and size selection as in the 3′-sequencing with DRaqL-SC3-seq described above. Then, purified cDNA was dA-tailed in a 10-μl dA tailing mixture (1 μl of 10× ExTaq Buffer, 0.94 μl of 2.5 mM dNTP, 0.066 μl of ExTaqHS, 8 μl of the purified DNA) with incubation at 72°C for 10 min. The dA-tailed product (10 μl) was subjected to the Y-shaped adapter ligation by adding a 6.6-μl Y-adapter mixture (3.3 μl of 5×NEB Next Quick Ligation Reaction Buffer, 0.2 μl of 10 μM Y-shaped S5-index adapter [Y-shaped DNA formed by annealing rP7-Ycore and P5-S5-index oligos, 0.33 μl of 2,000 U/μl T4 DNA Ligase, and 2.8 μl of DDW) and incubating at 20°C for 15 min and at 72°C for 20 min. The ligation product was purified with a 0.9× volume of AxyPrep reagent, and subjected to additional PCR amplification in 20 μl of the S5-N7 PCR mixture with the thermal cycling program 5, followed by DNA purification with a 0.9× volume of AxyPrep reagent. The libraries were mixed and subjected to sequencing from the Read1 and Read2 primers (75-bp each) using a NextSeq 500/550 High-Output Kit, ver. 2.5 (150 Cycles) (Illumina). We applied 10 samples of single-cell cDNAs from each cell type to the Y-shaped adapter sequencing, and used samples of the top 5 uniquely mapped reads (1.8–4 M reads) for subsequent analyses, to remove samples with small numbers of mapped reads (<1.5 M reads).

Library preparation for DRaqL-SMART-Seq v4, DRaqL-smart-seq2 and DRaqL-Protease-Smart-seq2

Libraries of the cDNAs amplified with the Smart-seq2-based methods (DRaqL-Smart-seq v4, DRaqL-Smart-seq2 and DRaqL-Protease-Smart-seq2) were prepared using the library construction module of the SMART-Seq v4 3′ DE Kit, according to the manufacturer’s instructions. Briefly, 10–12 purified cDNA samples were pooled and purified with a 0.8× volume of AxyPrep reagent. 300–400 pg cDNAs were tagmented using the Amplicon Tagment Mix of the Illumina Nextera XT DNA Sample Preparation Kit (Illumina). 3′-ends of the tagmented cDNAs were amplified with Nextera PCR Master Mix and dual index primers (Reverse PCR primer HT Index and Tn Forward PCR Primer HT Index in the SMART-Seq v4 3′ DE Kit) designed for TnRP1 and TnRP2 sequences. All libraries were mixed and sequenced using the NextSeq 500/550 High-Output Kit v2.5 (75 Cycles) with paired-end mode using Read2 (8 bps) for the in-line index and with dual-index mode (8 bps).

Single-cell RNA-seq of mouse ovarian oocytes and granulosa cells with DRaqL-SC3seq

Single oocytes (72 cells) in the primary-to-early antral follicles and the single granulosa cells associated with these oocytes in secondary-to-early antral follicles (140 cells) were histologically identified in the frozen ovarian sections stained with cresyl violet, and isolated with LCM. cDNAs were amplified with DRaqL from the isolated single cells, and 45 oocytes and 94 granulosa cells were passed through the Smirnov–Grubbs test for outliers (P > 0.01). We removed an oocyte and granulosa cells with low mapping rates (<3% and <13%, respectively), and subjected 44 oocytes and 61 granulosa cells (105 samples in total) to the 3′-sequencing (Nakamura et al, 2015), resulting in ∼2.2 M and ∼3.2 M uniquely mapped reads per single oocyte and granulosa cell, respectively (Fig S5). One cell showed a small number of detectable protein-coding genes (3,138) and low correlation coefficients of gene expression levels with all other cells in this study (<0.3), and was removed from subsequent analyses. In addition, after the sequence data processing described below, granulosa cells that showed mixed expression profiles with oocytes (3/60: 5%) were excluded from subsequent analyses. Details of evaluation of the mixed profiles are described below. Therefore, a total of 44 oocytes and 57 granulosa cells were finally used for the 3′-sequencing analyses. To examine whether known splice isoforms in oogenesis were detected, we subjected the cDNAs of 10 oocytes to the Y-shaped adapter sequencing and applied the top 5 samples with respect to the number of mapped reads to the subsequent analyses according to the above-described criteria.

Data processing and transcriptome analysis for 3′-sequencing

Data processing for 3′ sequencing was performed according to the previously reported method with minor modifications. FASTQ files were generated using the bcl2fastq tool (v. 2.20.0.422) and were trimmed to remove polyA/T sequences, adapter sequences, and poor-quality reads, using fastp (v. 0.20.0) (Chen et al, 2018) and cutadapt (ver. 1.16) (Martin, 2011). Then, reads were aligned to the GRCm38 mouse genome assembly using the HiSAT2 tool (v. 2.1.0) (Kim et al, 2015). Best-matched reads on annotated genes were counted using HTSeq 0.10.0 (Anders et al, 2015; Putri et al, 2021 Preprint). Gene annotation was performed as described previously (Nakamura et al, 2015), and reads mapped within 10 kb from the 3′-ends of genes annotated in the Ensembl database were included in the transcript counts of the genes.

For transcriptome analyses, read counts RPM matched to each gene were converted to log₂ values (log₂ [RPM+1]). Protein-coding genes with at least one mapped read were considered to be detectable genes. The detection rate of a gene in a group of cells was defined as the frequency of samples wherein the gene was detected. The two proportion z-test (two-sided) was performed to analyze the statistical significance of the difference of detection rates. The t test (two sided) for differences in average expression levels was performed using log₂ (RPM+1) values. False discovery rates of multiple testing were calculated using the Benjamini–Hochberg method (Benjamini and Hochberg, 1995). DEGs in granulosa cells were identified by the following criteria: difference of averaged log₂ (RPM+1) >1.8 (i.e., >3.5-fold), which was employed according to the previous study using the same cDNA amplification method (Kurimoto et al, 2008); log₂ (RPM+1) >4 in at least one sample; P-value by t test < 0.05 and FDR <0.05. Enrichment of gene ontology terms in DEGs was calculated using the DAVID tool (Huang da et al, 2009a; Huang da et al, 2009b). PCA was performed for genes showing log₂ (RPM+1) >4 in at least one sample, using the prcmp tool in the R program suite (R Core Team, 2020) with default settings, respectively. 97% confidence interval ellipses were calculated using the stat_ellipse tool in the R program suite.

Data processing and exon profiling for Y-shaped adapter sequencing

Sequence trimming and genomic mapping were performed as described above. Reads mapped to exons and exon–exon junctions were counted for their expression profiling by featureCounts (v2.0.1) (Liao et al, 2014) with -J -M -t exon -g gene_id options. Taking account of the size variations among exons, we quantified the expression levels of exons and exon–exon junctions using TMM values, rather than RPM. The read counts were normalized using EdgeR (3.36.0) (Robinson et al, 2010), and transformed to TMM for expression level analyses. For expression profiling of exon–exon junctions, only genes with detectable expression levels in the 3′-sequencing were used, and exons and their junctions showing >2 TMM values were defined as detectable. To analyze expression levels, TMMs were transformed to log₂ values (log₂ [TMM+1]). To define the number of exons for each gene, we calculated the detection rates of exons for every transcript from each protein-coding gene registered in the Ensembl database, among all samples, and selected the transcript that showed the largest detection rate. The sequencing statistics for mESCs and cells of mouse ovaries (oocytes and granulosa cells) analyzed in this study are shown in Tables S6 and S7, respectively.

Data processing for DRaqL-Smart-Seq v4, DRaqL-Smart-seq2 and DRaqL-Protease-Smart-seq2

The binary base call sequence files were converted to fastq files by the bcl2fastq tool with --minimum-trimmed-read-length 0 --mask-short-adapter-reads 0 --no-lane-splitting options. The fastq files were demultiplexed by using a SMART-Seq DE3 Demultiplexer. Quality control, read mapping, and read counting were performed for Read1 sequences by the same method as the data processing of the DRaqL-SC3-seq datasets.

Statistical model of the oocyte transcriptome based on diameter

To reconstruct the transcriptome of oocytes, we employed simple linear regression analyses. PC1 and PC2 of oocytes served as objective variables, whereas their corresponding diameter was taken as an explanatory variable. We calculated PC1 and PC2 for each diameter bin, ranging from 10 μm to 100 μm at 10-μm intervals. The gene expression level of each gene was then reconstructed by adding the products of the calculated PC values and factor loadings to the mean log₂-expression level. This process allowed us to create the reconstructed transcriptome for oocyte corresponding to each diameter bin. Subsequently, we calculated the Spearman’s rank correlation coefficient (rs) between the reconstructed transcriptomes and the transcriptome of each individual oocyte. The transcriptome with the highest rs value was considered the best-matched reconstructed transcriptome for that specific oocyte. We annotated oocytes with diameters significantly different from those of their best-matched models (>20 μm difference). To assess the similarity between the reconstructed transcriptomes and transcriptomes of real oocytes, we selected the top 500 genes with the largest positive and negative contributions to PC2, and performed hierarchical clustering using the Ward.D2 method in the R program suite.

Estimation of RNA copy number

The mRNA copy number was estimated for each gene using the spike-in RNAs (ERCC). The read number of each gene was fitted to the linear regression model of log₂-transformed read numbers and known copy numbers of ERCC RNAs, using the lm tool of the R program suite. The read numbers of the ERCC RNAs were obtained from the mapped sequence data using the idxstats function in samtools. ERCC RNAs that showed >0 reads were used for the calculation of linear regression models.

Single-cell transcriptome analysis of freshly dissociated oocytes and granulosa cells, and those isolated from frozen sections through LCM, from ovaries in the same mouse

The estrous cycle of an 8-wk-old mouse was determined using vaginal smear analysis. The mouse was euthanized during proestrus, and its ovaries were isolated in L15 medium (L5520; Sigma-Aldrich). One ovary was immediately frozen in liquid N₂ for subsequent LCM-based single-cell transcriptome analysis using DRaqL-Smart-seq2.

From the other ovary, COCs were isolated from fully grown follicles using tweezers and a needle into L15 medium containing 240 μM Dibutyryl-cAMP (#D0627; Sigma-Aldrich), an analog of cAMP. This COC isolation was conducted as previously reported (Takashima et al, 2021), with a minor modification in that natural estrus cycle was investigated without superovulation. The oocytes and their attached granulosa cells, and other granulosa cells within COCs, were dissociated by pipetting and then isolated using a glass capillary under visual inspection with a stereomicroscope. Each cell was washed twice in 200 μl of a 0.1% PVA solution in PBS. These freshly dissociated oocytes and granulosa cells were collected into cell lysis buffer, and cDNAs were amplified using Smart-seq2.

RNA sequencing and data processing were performed using the same method as described above, with slight modifications in the quality control process: granulosa cells with less than 100,000 reads and <3,000 detectable protein-coding genes, and oocytes with less than 100,000 reads and <10,000 detectable protein-coding genes, were excluded from further analysis. Consequently, a total of 14 out of 16 freshly dissociated oocytes, 8 out of 8 LCM-isolated oocytes from frozen sections, 35 out of 80 freshly dissociated granulosa cells, and 70 out of 80 LCM-isolated granulosa cells from frozen sections were included in the subsequent analysis.

Evaluation of mixed expression profiles of oocytes and granulosa cells

To assess the influence of oocyte transcriptome on adjacent granulosa cells in frozen sections analyzed using DRaqL-SC3-seq, we compared transcriptome of granulosa cells with the average expression profile of oocytes using scatterplots of log₂ gene expression levels with contour profiles. We observed a subset of granulosa cells (5%, 3/60) displaying bimodal patterns in the contour plots, indicating the presence of mixed expression profiles of granulosa cells and oocytes (see Fig S7A and B).

To further investigate the existence of mixed expression profiles, we utilized the single-cell RNA-seq data of freshly dissociated oocytes and granulosa cells described above, ensuring the analysis of genuine single cells without contamination. Consistent with the findings from the frozen sections, the contour plots of log₂ expression levels between granulosa cells and oocytes also exhibited bimodal patterns (Fig S8). Notably, granulosa cells directly attached to oocytes showed higher frequency of these mixed profiles (42% [6/14]) compared with other granulosa cells (5% [1/21]). Thus, we consider that these profiles are authentic expression profiles rather than a result of contamination. However, to ensure a focused investigation of the pure transcriptome of granulosa cells, we excluded the cells exhibiting mixed profiles from subsequent analyses using DRaqL-SC3-seq, resulting in the detailed exploration of the remaining 57 granulosa cells (Figs 4, 5, and 6).

Analysis of published datasets for oocytes and granulosa cells

A previously published single-cell RNA-seq dataset for mouse ovarian somatic cells (CRA003928) (Li et al, 2021) was downloaded and re-analyzed using the Scanpy tool (Wolf et al, 2018) in the Seurat v4.1.1 program suite (Hao et al, 2021) with the previously reported settings (Li et al, 2021). The features were filtered as follows: 200 < nFeature < 5,000, 1 < percent.mt < 5 and nCount_RNA > 1,000. The filtered feature counts were normalized by the total number of tags, and the counts per 10,000 reads were calculated using the NormalizedData function with LogNormalize method. Highly variable genes for subsequent analyses (6,883 genes) were identified by the FindVariableFeatures tool with the dispersion method and variable nfeatures. UMAP was calculated using the RunUMAP tool with dims = 1:50, and identified 12 clusters, which were then annotated using the marker genes as reported previously (Li et al, 2021). Cluster 1 was identified as cumulus cells based on the expression of Fabp5 and Ldha, and low/no detectable expression of corpus luteal marker (Hsd3b1, Inhba) and early granulosa markers in pre-antral follicles (Birc5, Cdc8). Among all genes annotated, the DEGs between the neighboring and non-neighboring granulosa cells identified in the DRaqL-SC3-seq analysis (Fig S16) were extracted for further analysis. The average expression levels of the genes up-regulated in the neighboring and non-neighboring cells were calculated, and cells were ordered according to the rank of the difference between these expression levels. The expression level differences between the top and bottom 25% of cells were examined with the Wilcoxon rank test.

For the analysis of single-cell RNA-seq of human oocytes and granulosa cells (GSE186504) (Fan et al, 2021), expression data were downloaded, and 130 granulosa cells were identified based on the information in SRA Run Selector and their gene expression profiles. The average expression levels (RPM values) of the above DEGs were calculated and ordered according to the rank of the expression level difference followed by the Wilcoxon rank test as described above.

For the comparison with transcriptome data during mouse oogenesis from nongrowing to germinal vesicle (GV) oocytes (Gahurova et al, 2017), we downloaded RNA-seq read count data from GSE86297 and calculated log₂ (sample read counts per million reads + 1) values. To identify DEGs among nongrowing and growing oocytes, we performed analysis of valiance and identified 1,372 genes that exhibited P < 0.05, a minimum log₂ expression level >4 in at least one sample, and a minimum log₂ expression level difference > 1 in at least one pair of samples. Among these genes, we identified 1,286 whose gene symbols were found in the gene table used in our study (GRCm39). For these genes, we calculated correlation coefficients with our transcriptome data of single oocytes.

For evaluation of the activation phase of granulosa cells, we investigated genes identified by Morris et al (2022) (Supplementary File 2: top 10 markers expressed in each ovary cluster) (Morris et al, 2022). We analyzed log₂ RPM+1 values in the single granulosa cells and oocytes in our dataset for “Granulosa,” “Preantral-Cumulus,” “Antral-Mural,” “Atretic,” “Mitotic,” “Luteinizing mural,” “Active CL,” “Mitotic-Antral,” and “Regressing CL” (CL, corpus luteum).

For the comparison with transcriptome data of human oocytes published by Ernst et al (2017), we downloaded data of DEGs between oocytes in primordial and primary follicles (Table S1 in Ernst et al [2017]: DEGs). Human Ensembl gene IDs were converted to mouse gene ID using g:Profiler (Raudvere et al, 2019). 249 and 186 genes were identified as up- and down-regulated in human primary-follicle oocytes, respectively, sharing orthologues with mice. Z scores were calculated using the primary oocyte mean (FPKM) values of these genes. We also calculated Z scores for average log₂ RPM+1 values in the primary-follicle oocytes in our dataset for mouse orthologues. These Z scores were compared between human and mouse, and genes that exhibited more than twofold differences were identified as differentially expressed. Gene ontology analysis was performed for DEGs between humans and mice using the DAVID tool.

Immunofluorescence

Snap-frozen mouse ovary embedded in optimal cutting temperature compound (Tissue-Tek) were sectioned with a 6-μm thickness and mounted on MAS-GP typeA glass slides (MATSUNAMI). Ovarian sections were fixed with 50% iso-propanol for 30 s at RT followed by permeabilization in 1×PBS containing 1% Triton X-100 (12967-32; nacalai tesque) for 10 min and three times of 5-min washing in 1×PBS. Then, the ovarian sections were incubated with 5% BSA (Sigma-Aldrich) and 10% in 1×PBS donkey serum (Sigma-Aldrich) for 60 min at RT for blocking. Then, sections were incubated with primary antibodies in 1×PBS containing 1% BSA in a dark humidity chamber for 1 h at RT, followed by three 5-min washes in 1×PBS. The following primary antibodies were used at a dilution rate of 1:00 in this study: rabbit anti-PBX1 antibody (18204–1-AP; Proteintech): rabbit anti-SUZ12 antibody (3737T; Cell Signaling Technology). Then, the sections were incubated with secondary antibody (donkey anti-rabbit Alexa Fluor 594 #A21207; Life Technologies) at a dilution rate of 1:300 in 1× PBS containing 1% BSA for 1 h at RT, followed by three 5-min washes in 1×PBS, and were mounted in VECTASHIELD Vibrance Antifade Mounting Medium with DAPI (VECTOR LABORATORIES) with a cover glass (MATSUNAMI).

These ovarian sections were analyzed with confocal laser microscopy using Olympus FV3000 with a ×20 objective lens. For quantitative analysis of fluorescent signals, regions of nuclei were identified with DAPI signals using ImageJ ROI manager, and signal densities of the target proteins in nuclei were quantified using the defined ROIs. Four follicles were analyzed for the measurement of fluorescence signal intensities for both proteins. Robust z-score of the signal intensities were calculated and Wilcoxon test were performed by the R program suit.