Small RNA sequencing reveals distinct nuclear microRNAs in pig granulosa cells during ovarian follicle growth
Journal of Ovarian Research volume 14, Article number: 54 (2021)
Nuclear small RNAs have emerged as an important subset of non-coding RNA species that are capable of regulating gene expression. A type of small RNA, microRNA (miRNA) have been shown to regulate development of the ovarian follicle via canonical targeting and translational repression. Little has been done to study these molecules at a subcellular level. Using cell fractionation and high throughput sequencing, we surveyed cytoplasmic and nuclear small RNA found in the granulosa cells of the pig ovarian antral preovulatory follicle. Bioinformatics analysis revealed a diverse network of small RNA that differ in their subcellular distribution and implied function. We identified predicted genomic DNA binding sites for nucleus-enriched miRNAs that may potentially be involved in transcriptional regulation. The small nucleolar RNA (snoRNA) SNORA73, known to be involved in steroid synthesis, was also found to be highly enriched in the cytoplasm, suggesting a role for snoRNA species in ovarian function. Taken together, these data provide an important resource to study the small RNAome in ovarian follicles and how they may impact fertility.
Increasing ovarian follicle size has long been an important indicator of oocyte quality, with oocytes obtained from large follicles showing consistently higher in vitro and in vivo developmental potential than their counterparts from small follicles [3, 51]. In addition to increased expression of key growth factors like GDF9 and BMP4, we and others have shown differences in microRNA (miRNA) expression in the somatic cells of these follicles [44,45,46, 54, 55, 57]. Recent investigations have looked at the distinct populations of small RNA that perform key biological functions in the cytoplasmic and nuclear compartments of cells [19, 41, 47, 60]. All of these RNA species are transcribed from chromosomal DNA, and most are processed in both the nucleus and cytoplasm. Accordingly, a complex system of transport exists to dynamically localize small RNA within the cell. For example, after processing into mature miRNA by Dicer in the cytoplasm, strands loaded onto Ago2 may be subsequently imported back into the nucleus via Importin 8 . The availability of complex secondary structures between RNA and DNA species permit these RNA to function in a myriad of ways. In the nucleus, miRNA can form triple helices with chromosomal DNA, to regulate gene expression . Even degraded fragments of miscoded messenger RNA have been shown to participate in the assembly of such scaffolds . Although small nucleolar RNAs (snoRNAs) are principally involved with RNA modification in the nucleus, they are shuttled to the cytoplasm during periods of stress [20, 31]. Despite evidence that such RNA species have functional roles in both the nuclear and cytoplasmic compartments of the cell, little systematic analysis has been done to study them at a population level. Here, we present a survey of small RNA species present in gilt granulosa cells obtained from small and large preovulatory follicles, herein referred to as SGCs and LGCs. While we found a few significant differences between SGCs and LGCs, we did reveal a diverse network of small RNA that showed distinct subcellular localization. This small RNAome of ovarian granulosa cells will provide an important resource for studying their subcellular function during follicle growth.
Validation of subcellular fractionation
Following cell fractionation, intact nuclei were visualized by phase-contrast and fluorescence microscopy, after chromosome staining with Hoescht 33342 (Fig. 1a). Before proceeding with deep sequencing of small RNAs, we confirmed the quality and purity of the nuclear and cytoplasmic fractions. The RNA integrity number (RIN) for each sample was determined by capillary electrophoresis: LGC cytosol, 8.62 ± 0.40; LGC nucleus, 6.92 ± 0.17; SGC cytosol, 9.02 ± 0.22; SGC nucleus, 6.78 ± 0.46. Western blotting revealed the purity of each fraction by the exclusion of glyceraldehyde 3-phosphate dehydrogenase from the nucleus and lamin B from the cytosol (Fig. 1b). RT-qPCR analysis of the nuclear spliceosomal RNA U6 revealed 1000-fold higher expression in the nucleus, as expected (Fig. 1c), confirming clear separation of the nucleus and cytosol fractions.
Overview of small RNA sequencing data
Four nuclear and four cytoplasmic RNA pools were sequenced from both LGCs and SGCs. After trimming adapter sequences, reads with a sequence length less than 18 nt were discarded. Remaining sequencing reads had lengths between 18 and 50 nt (Fig. 2a). The majority of reads (97%) were between 22 and 25 nt in length. Twenty percent of these reads could not be mapped to the pig genome and 7% mapped to coding regions.
We mapped the remaining reads to S. scrofa small RNA databases and looked broadly at RNA species breakdown in nuclear and cytoplasmic fractions (Fig. 2b). MicroRNA, both mature sequences and hairpin precursors, made up the majority of the mapped reads in both cell fractions (55% of the nuclear reads; 61% of cytoplasmic reads) followed by 30.4% of nuclear reads mapping to snoRNA, while 16.3% of cytoplasmic reads remained uncharacterized as “miscellaneous RNA”. As expected, considerably fewer cytoplasmic reads mapped to snoRNA, while reads mapping to ribozyme, transfer RNA and mitochondrial RNA were five- to ten-fold lower in the nuclear fraction compared to cytoplasmic fractions.
MicroRNAs have distinct subcellular localization patterns
After mapping reads to annotated swine miRNA, we looked at the relationship between samples using principal component analysis and unsupervised clustering. The first principle component comprised 54% of the variance and separated nuclear samples from cytoplasmic samples. The second principal component that comprised 11% of the variance in the data was related to batch effect, although interestingly this primarily affected the variance of nuclear samples. This batch effect was corrected for in subsequent statistical analysis. Expression of miRNA from distinct subcellular components clustered together without any obvious differences between miRNA from SGC and LGC (Fig. 3a). We next examined expression of miRNA using a generalized linear model (see Methods) comparing all four groups (SGC nuclei, LGC nuclei, SGC cytoplasm, and LGC cytoplasm). While most miRNA that showed differences in expression between the nucleus and cytoplasm were common between SGC and LGC, several miRNA differed in their subcellular localization between granulosa cells from the two stages of follicular development. Further examination of nuclear miRNA revealed seven differentially expressed between SGC and LGC (Fig. 3b).
We then ranked miRNA based on their fold change between the nucleus and cytoplasm (Supplemental Tables 1 and 2). SGC showed more significantly different miRNAs between these subcellular compartments than LGC (101 versus 83, respectively).
To confirm our sequencing results, we analyzed expression of 12 miRNA by digital droplet RT-PCR in our fractionated samples. All analyzed miRNA showed identical enrichment, either cytoplasmic or nuclear, and a significant correlation of expression ratios (R = 0.50, P = 0.013) between the two technologies (Fig. 4).
Some nuclear miRNAs are predicted to target promoter regions
With specific miRNAs observed to be enriched in the nucleus of granulosa cells and previous work demonstrating that small RNA in the nucleus can regulate transcription, we examined the potential targeting of these nuclear miRNAs to genomic promoter regions. Because of a dearth of annotation for pig promoter regions, we used the fact that mature miRNAs are highly homologous  and looked first at binding sites in human promoter regions that were then mapped to homologous regions in the pig. Our analysis revealed 417 potential binding sites at the promoter regions across the human genome, of which 55 had identical matches to homologous regions in the pig genome (Supplemental Table 3 and Supplemental File 1). We then used all putative targets to conduct a gene set analysis  to check whether these miRNA were targeting common pathways. DNA-binding genes that positively regulate patterning, development and locomotion were all significantly enriched (Table 1).
Putative novel miRNA identified in pig granulosa cells
We took advantage of the miRDeep2 algorithm to predict novel miRNA based on precursor structure and read depth . Using all reads from all groups, 39 miRNA were predicted with a score above 10, of which 34 ± 2 were estimated to be true positives. From this list, we checked ten with the highest number of mature read counts (all greater than 500; median of 2841). Six sequences aligned to other RNA species (e.g. Y RNA, rRNA, tRNA), while one appeared to be a single nucleotide mismatch to ssc-miR-26a. The three remaining sequences shown in Fig. 5 likely represent novel pig miRNA based on structure and seed sequence similarity to known human miRNAs: hsa-miR-130a-3p, hsa-miR-3193 and hsa-miR-5693.
Other small RNA species are also present in the cytoplasm and nucleus
As we had observed significant numbers of reads mapping to snoRNA and piRNA (Fig. 2b), we looked further into their subcellular distributions in granulosa cells. We found that piRNA were more abundant in the nucleus and showed a greater variation between SGC and LGC here than in the cytoplasm (data not shown). As expected, no reads mapped to piRNA found on the Y chromosome.
Similarly, we looked at the distribution of snoRNA, which were also found throughout the cell. Most significantly different snoRNA were those found in the cytoplasm, made up exclusively of box H/ACA type snoRNA, while those found in the nucleus belonged to the box C/D class snoRNA (Fig. 6). One snoRNA, SNORD86 (accession RF00594), showed significant nuclear enrichment, but only in small granulosa cells.
Our expression profile in porcine granulosa cells compared small non-coding RNA at two stages of ovarian follicle development and looked at their distribution in the nucleus and cytoplasm. Read lengths were similarly distributed between all samples and mapping to specific RNA species followed an expected subcellular distribution. The large fraction of uncharacterized “miscellaneous reads” found in the cytoplasm (Fig. 2b) may be due to degraded RNA species. While most of these will be degraded (reviewed in ) recent evidence has suggested that some short RNAs have functional roles within the cell . With continued investigation in this field, a much broader role for these small RNA is likely to be uncovered. We observed a substantial batch effect on nuclear miRNA and snoRNA, suggesting an unstable nuclear component that could be sensitive to processing. Our statistical models take this into account but it would be interesting to further investigate whether this source of variation is technical or biological in nature.
Broadly, miRNA were not substantially different between SGCs and LGCs, nor was the expression of cytoplasmic miRNAs. We did, however, identify seven nuclear miRNAs that were differentially expressed between SGC and LGC (Fig. 3b). Of note, these microRNA have been shown to play roles in granulosa cell function, including proliferation and apoptosis. The miR-182 cluster has been shown to promote proliferation in bovine granulosa cells by directly targeting the forkhead box protein O1 (FOXO1) . Our data shows higher levels of miR-182 in LGC, consistent with a proliferative phenotype. Conversely, miR-181a showed a decrease in expression from SGC to LGC. This miRNA was demonstrated to increase FOXO1 acetylation by targeting the deacetylase sirtuin 1, ultimately promoting granulosa cell apoptosis . Another study found that this same miRNA targeted the activin receptor IIA (ACVR2A) in mouse granulosa cells and its expression promoted an anti-proliferative phenotype . Interestingly, this group also found that whole-cell expression of miR-181a decreased during follicle development, which is in line with our findings. In rat ovaries, miR-122 has been shown to target the luteinizing hormone receptor  although we found no change in cytoplasmic levels of this miRNA during follicular growth. We did see a considerable increase in nuclear expression of miR-122 in LGC. This miRNA has been shown to regulate the processing of pre-miR-21 in hepatocellular carcinoma cells and promoted apoptosis . Whether this same miRNA processing control is exerted in granulosa cells is intriguing, especially as we have recently identified miR-21 as playing an important role in cumulus-oocyte complex maturation . MiR-126-5p showed a similar trend as miR-122, namely increased nuclear expression in LGCs. This miRNA has been shown to target the follicle-stimulating receptor (FSHR) in pig granulosa cells, suppress cytochrome 19A1 (aromatase) expression, and promote apoptosis [8, 32]. Aromatase, and consequently estradiol production, is also negatively regulated by miR-378, as we have previously reported [45, 57, 64]. Other aspects of follicle development have also been reported to be regulated by miR-378 [23, 54, 59], although to our knowledge this is the first time this analysis has been undertaken at a subcellular level. Here it appears that the known decrease in miR-378 in porcine granulosa cells taking place during follicle growth is derived primarily from changes in nuclear expression. Further studies should reveal how nuclear miR-378 regulates granulosa function. Less well-studied, nuclear expression of miR-142-5p and miR-192 was higher in LGC; these two miRNA have been found to be increased in patients with polycystic ovarian syndrome and diminished ovarian reserve, respectively [33, 63], but their role in ovarian function is not known. Taken together, our results are supported by other studies of granulosa cell miRNAs. We did observe very significant differences in miRNA expression between the nucleus and cytoplasm when considering both groups together; a larger sample size may have had sufficient power to detect additional differentially expressed miRNA. Combining samples of both sizes (LGC and SGC) also allowed us to increase the reads available to miRDeep2 to predict novel miRNA. We identified three high-confidence novel miRNAs present in pig granulosa cells on chromosomes 2, 12, and 13. None of these showed a strong bias towards either subcellular compartment (data not shown), but as part of a complex miRNA network [13, 24] they are undoubtedly important in regulating granulosa cell behaviour and warrant future investigation into their precise functions. Cytoplasm-enriched miRNAs included well known, abundant members like let-7 isomirs and miR-27b . No changes in let-7 family expression associated with atresia in the pig follicle [4, 68] were found between SGC and LGC groups. We also did not observe any pattern in the sequences of mature miRNA between nucleus- and cytoplasm-enriched (data not shown). While a hexanucleotide motif has been shown to be important for miR-29b , this has not been observed by others [5, 34], which agrees with our results presented here. Of note, we did observe strong nuclear enrichment of miR-29b in both SGC and LGC (Fig. 4a; Supplemental Tables 1 and 2). Other lines of evidence have suggested that the subcellular distribution of miRNA is dependent on target abundance, and miRNA are shuttled to the appropriate cellular compartment as needed . Argonaute-loaded miRNAs have been found to target nuclear RNA species, including introns and lncRNAs, which supports the idea that post-transcriptional regulation occurs in the nucleus and cytoplasm [14, 61]. The complete list of determinants for nuclear import of miRNA, however, remain to be elucidated.
In addition to post-transcriptional repression, a more recently appreciated role of nuclear miRNA is RNA activation (RNAa) whereby these small RNA can be recruited to complex nucleic acid scaffolds to recruit and stabilize transcriptional co-factors . Gene ontology analysis of promoter regions targeted by nucleus-enriched miRNA revealed enrichment of direct DNA-binding proteins. Parallel mechanisms involving either small RNA or DNA-binding proteins have been shown to be responsible for chromatin modifications in yeast . While no reports in mammalian cells link these two phenomena, the involvement of small RNA in heterochromatin is thought to be conserved (reviewed in ) and an autoregulatory loop between DNA-binding proteins and miRNA has been observed . The miRNA observed in our study that bind promoters of DNA-binding proteins could serve as redundant or feed-forward mechanisms to regulate transcription in granulosa cells during the final stages of follicle maturation. Several genes were associated with developmental processes, including the ubiquitously expressed upstream binding transcription factor (UBTF), a putative transcriptional target of let-7e that has been shown to regulate transcription to maintain homeostasis during cell growth . A study of human follicles showed higher oocyte expression levels of UBTF in antral compared with secondary follicles , consistent with a role in cell proliferation. A second predicted developmental target of let-7e is reticulon (RTN)4. While no in-depth studies exist for this family of endoplasmic reticulum-associated proteins in the ovary, other reticulon genes and receptor-like proteins in the follicle have been correlated to better embryo development in vitro . In mouse granulosa cells, RTN4 was a direct target of AKT1, itself a critically important regulator of cell growth and organized locomotion . Roles for let-7 family miRNAs other than let-7e have been demonstrated in the pig ovarian follicle [4, 68]. Given their high sequence homology, it is likely that some degree of target redundancy exists between members of this family. The steroid receptor RNA activator (SRA)1 gene was a predicted target of miR-339-5p. Knockdown of SRA1 in ovarian endometriosis increased estrogen receptor (ER)-alpha expression at the expense of ER-beta . Although SRA1 is expressed in the ovary, neither the long noncoding RNA nor the protein it codes for have been studied in granulosa cells. In ovarian cancer cells, miR-339-5p inhibits proliferation , although its role in normal follicle growth has not been established. These putatively targeted genes are also involved with biological processes that include differentiation and signaling in neurogenesis. Given the role of neurotrophins in ovarian development [7, 27, 37], our results suggest miRNA-mediated RNAa in granulosa cells may be another mechanism by which these small RNA regulate follicle growth.
We also profiled the snoRNA and piRNA that were captured by our sequencing data. The most studied function of mammalian piRNAs is the maintenance of genomic integrity in gametes . In spite of this, changes in piRNA populations have been implicated in ovarian cancers and low ovarian follicle reserves [2, 6, 56]. The repetitive nature of piRNA make mapping short reads to the genome difficult, but chromosomal piRNA expression appeared to be consistent between subcellular compartments. We also observed higher levels of piRNA and a larger difference between SGC and LGC in the nucleus. As the primary transposon silencing activity of piRNA takes place in the cytoplasm of gametes, our results could suggest nuclear functions like transcriptional gene silencing predominate in granulosa cells . A similar pattern of nuclear piRNA expression has been observed in cancer cells , although discerning whether these piRNA represent nascent transcripts originating in the nucleus or mature sequencing functioning there requires further study.
A similarly-limited focus has been placed on the role of snoRNA in the ovary. As both intact structures, and as a source of short regulatory RNAs, the role of snoRNA in cells continues to expand . A study looking at the box H/ACA snoRNA U17/SNORA73 found that it inhibits hypoxia-upregulated mitochondrial movement regulator (HUMMR) during cholesterol trafficking . Indeed, the authors find a reciprocal relationship between decreasing SNORA73 and increasing HUMMR during mouse ovary development; SNORA73 inhibition increased levels of ovarian pregnenolone and progesterone . This role in steroid synthesis is supported by our finding that SNORA73 is one of the most significantly cytoplasm-enriched snoRNA in granulosa cells (Fig. 5; accession RF00045). SnoRNAs may yet prove to be play an important role in ovarian function.
Overall, we have provided a survey of the subcellular small RNAome in pig preovulatory follicle granulosa cells. The data used for this analysis is freely available at github.com/derektoms/s3RNA, and we anticipate that it will be a valuable resource for others studying the role of small RNA in the pig ovary. Increasing granularity in small RNA profiling, including from single cells , will continue to increase our understanding of the complex role these molecules play in ovarian granulosa cells.
Granulosa cell isolation
Details of the isolation and culture of granulosa cells has been described previously . Briefly, porcine ovaries were removed from gilts at a local slaughterhouse and returned to the laboratory within 1.5 h in sterile 1x PBS at 22 °C. The ovaries were rinsed at least three times with 1x PBS. Granulosa cells were removed from small-sized follicles (1–3 mm in diameter) and large-sized follicles (3–6 mm in diameter) with a 20-gauge needle fixed to a 20-ml disposable syringe. Each sample represents a biologically distinct pool of granulosa cells from follicles of multiple ovaries that were subsequently fractioned.
Preparation of cytosolic and nuclear fractions
The separation of cytoplasmic and nuclear was performed in accordance with the manufacture’s guidelines (Nuclei EZ Prep Nuclei Isolation Kit, Product Code: NUC-101). Briefly, fresh granulosa cells were transferred into a separate 15-ml centrifuge tube and centrifuged at 500 x g for 5 min at 4 °C. Cells were kept on ice and washed with 10 ml ice cold PBS twice. Nuclei EZ lysis buffer and brief vortexing was used to lyse cells; nuclei were collected by centrifugation at 500 x g for 5 min at 4 °C. The cytoplasmic fraction was selected from the upper 250 ul and the nuclear fraction from the lower 45 ul, discarding an intermediate layer of around 200 ul with the intention of controlling cross-contamination. Cold Nuclei EZ storage buffer was added to nuclear fraction, followed by further vortexing and trituration to help break up clumps of nuclei.
Nuclei for RNA and protein analysis were immediately frozen at − 80 °C in Nuclei EZ storage buffer. For assessing morphological quality and yield, nuclei were placed on a glass slide and microscopically examined. The nuclei were stained with either Hoechst 33342 (10 ng μl− 1), or 0.2% (vol/vol) trypan blue solution in PBS to confirm the purity.
Western blotting analysis and antibodies
The cell nuclear pellets were lysed by using the Immunoprecipitation & Western Blotting RIPA Lysis Buffer (Cat 20–188, Sigma USA). Western blots were carried out using 4-20% Mini-PROTEAN® TGX Stain-Free™ Precast Gels (Cat. 4568095, Bio-Rad, Hercules, CA, USA) loaded with 50 to 100 μg of lysate and subsequently transferred to 0.45 μm Immun-Blot® low fluorescence PVDF Membrane (Cat. 1620260, Bio-Rad) using a Trans-Blot Turbo apparatus (Bio-Rad). Gels were activated by ultraviolet (UV) exposure for 2 min using a Bio-Rad ChemiDoc MP imager and total protein for each loading well were assessed. After protein transfer, membranes were imaged for stain-Free staining and total protein was quantified using Image-Lab 4.1 software (Bio-Rad). After transfer, the gels were stained with Bio-Safe™ Coomassie Stain (Cat. 1610786, BioRad) for 15 to 60 min and quickly de-stained in water to remove nonspecific staining. The membranes were then imaged, and total leftover proteins were measured using Image-Lab 4.1. The transferred membranes were blocked with 5% (w/v) nonfat milk in TBST for 1 h and then were cut into two separate pieces according to the different sizes. The upper blot (bigger than 50 kDa) was incubated with the nuclear envelope marker Anti-Lamin B1 antibody (1:4,000, Abcam, ab16048, Rabbit polyclonal, Cambridge, MA) followed by incubation with secondary antibody (anti-Rabbit IgG horseradish peroxidase, 1:4,000; Cell Signaling Technology, Danvers, MA, USA) for 1 h at room temperature; The lower blots were incubated with the Anti-GAPDH antibody [6C5] (1:4,000, abcam, ab8245, Mouse monoclonal, 1:2,000, MA, USA) for 1 h, following by incubation with secondary antibody (anti-mouse IgG horseradish peroxidase, 1:4,000; Cell Signaling Technology, Danvers, MA, USA). Membranes were again washed with TBST three times for 10 min each, incubated with Clarity chemiluminescence substrate (Bio-Rad), and imaged on the ChemiDoc MP, and then the bands were detected and analyzed with ImageLab 4.1.
RNA isolation and quality control
Total RNA Purification Kit (Norgen Biotek, Cat. 17200; Thorold, ON, Canada) was used to isolate total RNA from both cytoplasm and nuclear samples. Briefly, around 3 × 106 granulosa nucleus or 100 μl cytoplasm lysate were mixed with 350 μL of Buffer RL, then vortexed for up to 60 s to ensure the mixture becomes transparent. 200 μL of 96–100% ethanol was added into to the lysate following by vortexing for 10–30 s. Lysate with ethanol was added onto the column and centrifuge for 1 min at 3,500 x g to allow the entire lysate to pass through the column. RNase-Free DNase I Kit was used to minimize genomic DNA contamination in according with the manufacture’s On-Column DNA Removal Protocol (Norgen, Cat. 25710). After genomic DNA removal, the column was washed and centrifuge twice. 50 μL of elution solution was used to recover the total RNA. Quality and quantity of total RNA were determined using a NanoDrop ND-1000 Spectrophotometer at the Genomics Facility in the University of Guelph’s Advanced Analysis Centre. Using Agilent 2100 Bioanalyzer DNA and RNA chips, we verified that DNA concentration in the nuclear fraction was >100X that of the cytoplasmic fraction and that 28S ribosomal RNA peak fluorescence in the cytoplasmic fraction was >10X that of the nuclear fraction. These two ratios for the same samples imply good fractionation.
Library preparation and sequencing
Small RNA library preparation and sequencing was performed by Personalized Genomics and Innovative Medicine in Toronto, Canada.
Briefly, small RNAs ranging from 18 to 52-nt were gel-purified and ligated to the 39 adaptor (59-pUCGUAUGCCGUCUUCUGCUUGidT-39; p, phosphate; idT, inverted deoxythymidine) and 59 adaptors (59-GUUCAGAGUUCUACAGUCCGACGAUC-39). Ligation products were gel-purified, reverse transcribed, and amplified using Illumina’s sRNA primer set (59-CAAGCAGAAGACG GCATACGA-39; 59-AATGATACGGCGACCACCGA39). Samples were sequenced on an Illumina 1G Genome Analyzer.
Samples were sequenced on the rapid run mode of Illumina HiSeq 2500 sequencing platform. Base calling (BCL) files from sequencer were converted to the FASTQ files using bcl2fastq2-v220.127.116.11 tool. The raw reads generated were single-end and of length 51-nt. Raw reads quality in FASTQ format was assessed using FastQC tool v0.11.0  and identified over-represented sequences (e.g. small RNA adapter and RNA PCR primer sequences). Small RNA sequencing data analysis conducted using miRDeep2 package . Mapper module from miRDeep2 package requires raw reads (FASTQ format) as the input. Mapper module contains the data preprocessing functionalities that involves the conversion of FASTQ to FASTA, RNA to DNA, remove non-ATGCUN/atgcun bases, trimming Illumina small RNA adapter sequences (TGGAATTCTCGGGTGCCAAGG) and RNA primer sequences, filtering low-quality reads with < 18 nt length and collapsing the identical reads to remove redundancy. To quantify the number of reads independently mapped to pig mature and hairpin miRNAs sequences from miRBase v21 [17, 29, 30], Ensembl CDS, other classes of non-coding RNAs [lincRNA, scaRNA, snoRNAs, snRNAs, rRNAs, sRNAs, tRNAs, mitochondrial tRNAs, mitochondrial rRNAs, miscRNAs; piRNAs] and finally to pig genome (http://igenomes.illumina.com.s3-website-us-east-1.amazonaws.com/Sus_scrofa/UCSC/susScr3) using the mapper module of miRDeep2 package with the following parameters (−e -h -i -j -k TGGAATTCTCGGGTGCCAAGG -l 18 -m -p –s –t –v). Followed by alignment of the collapsed reads against the pig reference genome (susScr3) using bowtie to generate mapped output in ARF format and collapsed reads in fasta format. The core miRDeep2 module is used to detect known and novel small RNAs, that involves the following input files, the mapped output (ARF), collapsed reads (FASTA), pig reference genome (FASTA) and the miRBase v21 mature and hairpin reference sequences (FASTA). In-house perl scripts were written to detect and quantify all types of isomiRs, which were treated as reads for a single mature miRNA.
The R package DESeq2 (v 1.22. 2;) was used to detect differential expression in small RNA-Seq data, specifically miRNA and snoRNA. Matrices of read counts (small RNA vs samples) were analyzed by DESeq2 using two different models, one based on two variables (GC size and subcellular localization) and the other on a single “group” variable combining size and localization (e.g. “SGC-nuclear”). Both models accounted for batch effects. Statistically significant RNAs were be filtered using an adjusted p-value < 0.1 based on the negative binomial distribution.
For hierarchical cluster analysis, data was transformed using a regularized log method  before unsupervised clustering was performed with the pheatmap package based on the Euclidean distances between transformed expression values using a complete linkage method (v 1.0.1 2;). For other visualization and ranking of differentially expressed genes, shrinkage of the log fold change estimates (effect size) was performed using the apeglm method .
Prediction of miRNA binding to the promoter
Target genes were predicted using computational prediction software MicroPIR2 (biotech4/2). This software contains a database with over 80 million miRNA predicted targets in the promoter sequences of the human genome. MicroRNA targets were searched by miRNA name, with an average p-value ≤0.05, average conservation score ≥ 0.85, and a maximum number of four unpaired nucleotides. The binding pattern and sequence was observed with consideration towards near perfect matching in the seed region. Ensembl, a genetic database, was then used to determine if the predicted human promoter target sequence is conserved in Sus scrofa by identifying the sequence within the genome of Sus scrofa using the “BLAST” feature. The selected genes were then searched in the Ensembl Sus scrofa database to determine the target sequences’ relative location to the start codon, ATG in order to confirm the target sequence is indeed in the promoter region, which is ~ 1000 bp upstream of ATG and ~ 200 bp downstream of ATG.
Availability of data and materials
Raw sequences have been deposited to https://www.ncbi.nlm.nih.gov/sra/PRJNA604978 and the project is available at https://www.ncbi.nlm.nih.gov/bioproject/PRJNA604978
Processed read count data for miRNA and snoRNA, and R code used for analysis can be accessed at github.com/derektoms/s3RNA
Large preovulatory (4–6 mm diameter) follicle-derived granulosa cells;
Small preovulatory (0.5–3 mm diameter) follicle-derived granulosa cells
Small nucleolar RNA
Ahn HW, Morin RD, Zhao H, Harris RA, Coarfa C, Chen Z-J, et al. MicroRNA transcriptome in the newborn mouse ovaries determined by massive parallel sequencing. Mol Hum Reprod. 2010;7:463–71.
Bachmayr-Heyda A, Auer K, Sukhbaatar N, Aust S, Deycmar S, Reiner AT, et al. Small RNAs and the competing endogenous RNA network in high grade serous ovarian cancer tumor spread. Oncotarget. 2016.
Bagg MA, Nottle MB, Armstrong DT, Grupen CG. Relationship between follicle size and oocyte developmental competence in prepubertal and adult pigs. Reprod Fertil Dev. 2007;19(7):797. https://doi.org/10.1071/RD07018.
Cao R, Wu WJ, Zhou XL, Xiao P, Wang Y, Liu HL. Expression and preliminary functional profiling of the let-7 family during porcine ovary follicle atresia. Mol Cell. 2015;38(4):304–11. https://doi.org/10.14348/molcells.2015.2122.
Chen B, Zhang B, Luo H, Yuan J, Skogerbø G, Chen R. Distinct microRNA subcellular size and expression patterns in human cancer cells. Int J Cell Biol. 2012:672462.
Chen D, Zhang Z, Chen B, Ji D, Hao Y, Zhou P, et al. Altered microRNA and Piwi-interacting RNA profiles in cumulus cells from patients with diminished ovarian reserve. Biol Reprod. 2017;97(1):91–103. https://doi.org/10.1093/biolre/iox062.
Childs AJ, Bayne RAL, Murray AA, Martins Da Silva SJ, Collins CS, Spears N, et al. Differential expression and regulation by activin of the neurotrophins BDNF and NT4 during human and mouse ovarian development. Dev Dyn. 2010;239(4):1211–9. https://doi.org/10.1002/dvdy.22252.
Du X, Li Q, Pan Z, Li Q. Androgen receptor and miRNA-126∗ axis controls follicle-stimulating hormone receptor expression in porcine ovarian granulosa cells. Reproduction. 2016;152(2):161–9. https://doi.org/10.1530/REP-15-0517.
El-Brolosy MA, Kontarakis Z, Rossi A, Kuenne C, Günther S, Fukuda N, et al. Genetic compensation triggered by mutant mRNA degradation. Nature. 2019;7751:193–7.
Elzaiat M, Herman L, Legois B, Léger T, Todeschini AL, Veitia RA. High-throughput exploration of the network dependent on AKT1 in mouse ovarian granulosa cells. Mol Cell Proteomics. 2019;18(7):1307–19. https://doi.org/10.1074/mcp.RA119.001461.
Falaleeva M, Stamm S. Processing of snoRNAs as a new source of regulatory non-coding RNAs: snoRNA fragments form a new class of functional RNAs. Bioessays. 2013;1:46–54.
Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;1:37–52.
Friedman RC, Farh KK-H, Burge CB, Bartel DP. Most mammalian mRNAs are conserved targets of microRNAs. Genome Res. 2009;1:92–105.
Gagnon KT, Li L, Chu Y, Janowski BA, Corey DR. RNAi factors are present and active in human cell nuclei. Cell Rep. 2014;1:211–21.
Gebremedhn S, Salilew-Wondim D, Hoelker M, Rings F, Neuhoff C, Tholen E, et al. MicroRNA-183-96-182 cluster regulates bovine granulosa cell proliferation and cell cycle transition by coordinately targeting FOXO1. Biol Reprod. 2016;6:127.
Grewal SIS. RNAi-dependent formation of heterochromatin and its diverse functions. Curr Opin Genet Dev. 2010;20(2):134–41. https://doi.org/10.1016/j.gde.2010.02.003.
Griffiths-Jones S, Grocock RJ, van Dongen S, Bateman A, Enright AJ. miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res. 2006;34(90001):D140–4. https://doi.org/10.1093/nar/gkj112.
Hagemann-Jensen M, Abdullayev I, Sandberg R, Faridani OR. Small-seq for single-cell small-RNA sequencing. Nat Protoc. 2018;10:2407–24.
Hicks JA, Li L, Matsui M, Chu Y, Volkov O, Johnson KC, et al. Human GW182 Paralogs are the central organizers for RNA-mediated control of transcription. Cell Rep. 2017;7:1543–52.
Holley CL, Li MW, Scruggs BS, Matkovich SJ, Ory DS, Schaffer JE. Cytosolic accumulation of small nucleolar RNAs (snoRNAs) is dynamically regulated by NADPH oxidase. J Biol Chem. 2015;18:11741–8.
Houseley J, Tollervey D. The many pathways of RNA degradation. Cell. 2009;4:763–76.
Hwang H-W, Wentzel EA, Mendell JT. A hexanucleotide element directs microRNA nuclear import. Science (80-). 2007;5808:97–100.
Ioannidis J, Donadeu FX. Changes in circulating microRNA levels can be identified as early as day 8 of pregnancy in cattle. PLoS One. 2017;4:e0174892.
Ivanovska I, Cleary MA. Combinatorial microRNAs: working together to make a difference. Cell Cycle. 2008;20:3137–42.
Jia S. RNAi-independent heterochromatin nucleation by the stress-activated ATF/CREB family proteins. Science (80-). 2004;5679:1971–6.
Jinn S, Brandis KA, Ren A, Chacko A, Dudley-Rucker N, Gale SE, et al. SnoRNA U17 regulates cellular cholesterol trafficking. Cell Metab. 2015;6:855–67.
Kawamura K, Kawamura N, Mulders SM, Gelpke MDS, Hsueh AJW. Ovarian brain-derived neurotrophic factor (BDNF) promotes the development of oocytes into preimplantation embryos. Proc Natl Acad Sci U S A. 2005;102(26):9206–11. https://doi.org/10.1073/pnas.0502442102.
Kolde R. Pheatmap: pretty heatmaps; 2019.
Kozomara A, Griffiths-Jones S. MiRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 2014;42(D1):D68–73. https://doi.org/10.1093/nar/gkt1181.
Kozomara A, Birgaoanu M, Griffiths-Jones S. MiRBase: from microRNA sequences to function. Nucleic Acids Res. 2019.
Li MW, Sletten AC, Lee J, Pyles KD, Matkovich SJ, Ory DS, et al. Nuclear export factor 3 regulates localization of small nucleolar RNAs. J Biol Chem. 2017;49:20228–39.
Li Q, Du X, Liu L, Pan Z, Cao S, Li Q. MiR-126* is a novel functional target of transcription factor SMAD4 in ovarian granulosa cells. Gene. 2019.
Li Y, Xiang Y, Song Y, Wan L, Yu G, Tan L. Dysregulated miR-142, −33b and −423 in granulosa cells target TGFBR1 and SMAD7: a possible role in polycystic ovary syndrome. Mol Hum Reprod. 2019;10:638–46.
Liao J-Y, Ma L-M, Guo Y-H, Zhang Y-C, Zhou H, Shao P, et al. Deep sequencing of human nuclear and cytoplasmic small RNAs reveals an unexpectedly complex subcellular distribution of miRNAs and tRNA 3′ trailers. PLoS One. 2010;5(5):e10563. https://doi.org/10.1371/journal.pone.0010563.
Liao Y, Wang J, Jaehnig EJ, Shi Z, Zhang B. WebGestalt 2019: gene set analysis toolkit with revamped UIs and APIs. Nucleic Acids Res. 2019;47(W1):W199–205. https://doi.org/10.1093/nar/gkz401.
Lin K, Zhan H, Ma J, Xu K, Wu R, Zhou C, et al. Silencing of SRA1 regulates ER expression and attenuates the growth of stromal cells in ovarian endometriosis. Reprod Sci. 2017.
Linher-Melville K, Li J. The roles of glial cell line-derived neurotrophic factor, brain-derived neurotrophic factor and nerve growth factor during the final stage of folliculogenesis: a focus on oocyte maturation. Reproduction. 2013;2:R43–54.
Liu E-L, Zhou Y-X, Li J, Zhang D-H, Liang F. Long-chain non-coding RNA SNHG3 promotes the growth of ovarian cancer cells by targeting miR-339-5p/TRPC3 Axis. Onco Targets Ther. 2020:10959–71.
Liu H, Lei C, He Q, Pan Z, Xiao D, Tao Y. Nuclear functions of mammalian MicroRNAs in gene regulation, immunity and cancer. Mol Cancer. 2018;1.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550. https://doi.org/10.1186/s13059-014-0550-8.
Martinez G, Choudury SG, Slotkin RK. tRNA-derived small RNAs target transposable element transcripts. Nucleic Acids Res. 2017;9:5142–52.
Menon B, Sinden J, Franzo-Romain M, Botta RB, Menon KMJ. Regulation of LH receptor mRNA binding protein by miR-122 in rat ovaries. Endocrinology. 2013;12:4826–34.
van Montfoort APA, Geraedts JPM, Dumoulin JCM, Stassen APM, Evers JLH, Ayoubi TAY. Differential gene expression in cumulus cells as a prognostic indicator of embryo viability: a microarray analysis. Mol Hum Reprod. 2008.
Pan B, Li J. MicroRNA-21 up-regulates metalloprotease by down-regulating TIMP3 during cumulus cell-oocyte complex in vitro maturation. Mol Cell Endocrinol. 2018.
Pan B, Toms D, Shen W, Li J. MicroRNA-378 regulates oocyte maturation via the suppression of aromatase in porcine cumulus cells. Am J Physiol Endocrinol Metab. 2015;6:E525–34.
Pan B, Toms D, Li J. MicroRNA-574 suppresses oocyte maturation via targeting hyaluronan synthase 2 in porcine cumulus cells. Am J Phys Cell Physiol. 2017:ajpcell000652017.
Park CW, Zeng Y, Zhang X, Subramanian S, Steer CJ. Mature microRNAs identified in highly purified nuclei from HCT116 colon cancer cells. RNA Biol. 2010;5:606–14.
Paugh SW, Coss DR, Bao J, Laudermilk LT, Grace CR, Ferreira AM, et al. MicroRNAs form triplexes with double stranded DNA at sequence-specific binding sites; a eukaryotic mechanism via which microRNAs could directly alter gene expression. PLoS Comput Biol. 2016;2:e1004744.
Pitchiaya S, Heinicke LA, Park JI, Cameron EL, Walter NG. Resolving subcellular miRNA trafficking and turnover at single-molecule resolution. Cell Rep. 2017;3:630–42.
Post C, Clark JP, Sytnikova YA, Chirn GW, Lau NC. The capacity of target silencing by Drosophila PIWI and piRNAs. RNA. 2014;20(12):1977–86. https://doi.org/10.1261/rna.046300.114.
Rosen MP, Shen S, Dobson AT, Rinaudo PF, McCulloch CE, Cedars MI. A quantitative assessment of follicle size on oocyte developmental competence. Fertil Steril. 2008;90(3):684–90. https://doi.org/10.1016/j.fertnstert.2007.02.011.
Russell SJ, Stalker L, LaMarre J. PIWIs, piRNAs and Retrotransposons: complex battles during reprogramming in gametes and early embryos. Reprod Domest Anim. 2017;52:28–38. https://doi.org/10.1111/rda.13053.
Sanij E, Diesch J, Lesmana A, Poortinga G, Hein N, Lidgerwood G, et al. A novel role for the pol I transcription factor ubtf in maintaining genome stability through the regulation of highly transcribed pol II genes. Genome Res. 2015;25(2):201–12. https://doi.org/10.1101/gr.176115.114.
Schauer SN, Sontakke SD, Watson ED, Esteves CL, Donadeu FX. Involvement of miRNAs in equine follicle development. Reproduction. 2013;3:273–82.
Sen A, Prizant H, Light A, Biswas A, Hayes E, Lee H-J, et al. Androgens regulate ovarian follicular development by increasing follicle stimulating hormone receptor and microRNA-125b expression. Proc Natl Acad Sci U S A. 2014;8:3008–13.
Singh G, Roy J, Rout P, Mallick B. Genome-wide profiling of the PIWI-interacting RNA-mRNA regulatory networks in epithelial ovarian cancers. PLoS One. 2018.
Sun X-F, Li Y-P, Pan B, Wang Y-F, Li J, Shen W. Molecular regulation of miR-378 on the development of mouse follicle and the maturation of oocyte in vivo. Cell Cycle. 2018:1–13.
Sylvestre Y, De Guire V, Querido E, Mukhopadhyay UK, Bourdeau V, Major F, et al. An E2F/miR-20a autoregulatory feedback loop. J Biol Chem. 2007;4:2135–43.
Toms D, Xu S, Pan B, Wu D, Li J. Progesterone receptor expression in granulosa cells is suppressed by microRNA-378-3p. Mol Cell Endocrinol. 2014.
Wang D, Sun X, Wei Y, Liang H, Yuan M, Jin F, et al. Nuclear miR-122 directly regulates the biogenesis of cell survival oncomiR miR-21 at the posttranscriptional level. Nucleic Acids Res. 2018;4:2012–29.
Weinmann L, Höck J, Ivacevic T, Ohrt T, Mütze J, Schwille P, et al. Importin 8 is a gene silencing factor that targets argonaute proteins to distinct mRNAs. Cell. 2009;136(3):496–507. https://doi.org/10.1016/j.cell.2008.12.023.
Wingett SW, Andrews S. FastQ screen: a tool for multi-genome mapping and quality control. F1000Research. 2018.
Woo I, Christenson LK, Gunewardena S, Ingles SA, Thomas S, Ahmady A, et al. Micro-RNAs involved in cellular proliferation have altered expression profiles in granulosa of young women with diminished ovarian reserve. J Assist Reprod Genet. 2018;35(10):1777–86. https://doi.org/10.1007/s10815-018-1239-9.
Xu S, Linher-Melville K, Yang BB, Wu D, Li J. Micro-RNA378 (miR-378) regulates ovarian estradiol production by targeting aromatase. Endocrinology. 2011;10:3941–51.
Zhang M, Zhang Q, Hu Y, Xu L, Jiang Y, Zhang C, et al. miR-181a increases FoxO1 acetylation and promotes granulosa cell apoptosis via SIRT1 downregulation. Cell Death Dis. 2017;10:e3088.
Zhang Q, Sun H, Jiang Y, Ding L, Wu S, Fang T, et al. MicroRNA-181a suppresses mouse granulosa cell proliferation by targeting activin receptor IIA. PLoS One. 2013;3:e59667.
Zhang Y, Yan Z, Qin Q, Nisenblat V, Chang HM, Yu Y, et al. Transcriptome landscape of human folliculogenesis reveals oocyte and granulosa cell interactions. Mol Cell. 2018;72(6):1021–1034.e4. https://doi.org/10.1016/j.molcel.2018.10.029.
Zhou J, Liu J, Pan Z, Du X, Li X, Ma B, et al. The let-7g microRNA promotes follicular granulosa cell apoptosis by targeting transforming growth factor-β type 1 receptor. Mol Cell Endocrinol. 2015;409:103–12. https://doi.org/10.1016/j.mce.2015.03.012.
Zhu A, Ibrahim JG, Love MI. Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences. Bioinformatics. 2019;35(12):2084–92. https://doi.org/10.1093/bioinformatics/bty895.
The authors wish to thank the workers at Conestoga Farm Fresh pork for providing the ovaries used in this study and Ramanandan Prabhakaran at Personalized Genomics and Innovative Medicine for his help with sequencing and analysis.
This work was funded by a Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant [grant number 400259] to JL. The funding body had no input on the design of the study and collection, analysis, and interpretation of data, nor in writing the manuscript.
All animal tissue was obtained from a commercial slaughterhouse; no ethics approval was necessary.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Toms, D., Pan, B., Bai, Y. et al. Small RNA sequencing reveals distinct nuclear microRNAs in pig granulosa cells during ovarian follicle growth. J Ovarian Res 14, 54 (2021). https://doi.org/10.1186/s13048-021-00802-3
- Granulosa cells
- Reproductive biology
- Next generation sequencing
- Cellular fractionation