# Transcriptome Analysis of Bovine Ovarian Follicles at Predeviation and Onset of Deviation Stages of a Follicular Wave

Authors: Pengfei Li; Jinzhu Meng; Wenzhong Liu; George W. Smith; Jianbo Yao; Lihua Lyu
Journal: International Journal of Genomics
 (2016)

---

## Abstract

For two libraries (PDF1 and ODF1) using Illumina sequencing 44,082,301 and 43,708,132 clean reads were obtained, respectively. After being mapped to the bovine RefSeq database, 15,533 genes were identified to be expressed in both types of follicles (cut-off RPKM > 0.5), of which 719 were highly expressed in bovine follicles (cut-off RPKM > 100). Furthermore, 83 genes were identified as being differentially expressed in ODF1 versus PDF1, where 42 genes were upregulated and 41 genes were downregulated. KEGG pathway analysis revealed two upregulated genes in ODF1 versus PDF1, CYP11A1, and CYP19A1, which are important genes in the steroid hormone biosynthesis pathway. This study represents the first investigation of transcriptome of bovine follicles at predeviation and onset of deviation stages and provides a foundation for future investigation of the regulatory mechanisms involved in follicular development in cattle.


---

## Body

## 1. Introduction

The ovarian follicle is an essential component of the reproductive process. It plays an important role in controlling the estrous cycle, determining estrous behaviour, ensuring oocyte competency and subsequent embryo survival rate, and determining both postovulation corpus luteum function and progesterone synthesis [1]. In a number of species, follicular growth is characterized by a wave-like pattern, with two or three waves occurring during the normal course of estrous cycles in cattle [2]. During each wave of follicular development, a cohort of antral follicles are induced to begin accelerated growth [3]. After a period of concurrent growth, a species specific number of follicles will then be selected to become dominant, while the remaining follicles will be lost through a process known as atresia. Diameter deviation is defined as the divergence in growth rates between the two largest follicles in a follicular wave [3]. The onset of diameter deviation occurs when the largest follicle reaches 8.5 mm in dairy cattle and marks initiation of divergence in growth rate and estradiol producing capacity between the F1 or largest (future dominant follicle) and F2 or second largest (future subordinate) growing follicles culminating the process of dominant follicle selection. While the exact mechanisms of dominant follicle selection are not completely understood, there have been many studies on the hormones and factors involved in follicular development. Antral follicles are dependent upon FSH for growth and each follicular wave is preceded by a transient rise in FSH concentrations [4]. Many growth factors linked to regulation of follicular development, such as inhibins, activins, and insulin-like growth factors and their binding proteins, have been identified in follicular fluid of individual bovine follicles [5]. These molecules can regulate follicular cell survival, proliferation, or death. Recent studies have attempted to understand the molecular regulation of follicular development in cattle [6, 7]. However, the molecular mechanisms governing the wave-like pattern of follicular development are incompletely described, particularly at the onset of diameter deviation which is the first morphological indication of follicular dominance.Traditionally, gene expression studies in the field of follicular development focus on the study of expression of candidate genes of interest. With the development of next-generation sequencing technologies, transcriptome profiling has become a powerful approach for identification of genes globally expressed in various tissues including ovarian follicles [8]. In the present study, we performed RNA-Seq of granulosa cell RNA from bovine ovarian follicles at predeviation (PD) and onset of deviation (OD) stages of a follicular wave in cattle to catalog the transcriptome and identify potential differentially expressed genes associated with these key stages of a follicular wave. This study provides a comprehensive sequence resource for future studies on follicular development in cattle.


## 2. Materials and Methods

### 2.1. Materials

All materials were obtained from Sigma-Aldrich (St. Louis, MO) unless otherwise stated.


### 2.2. Animal Model and Sample Collection

All animal procedures were approved by the Institutional Animal Care and Use Committee at Michigan State University. Estrus was synchronized in nonlactating Holstein dairy cows with two injections of prostaglandin F2α (PGF2α; Prostamate; IVX Animal Health, St. Joseph, MO) administered 14 days apart, and follicular growth was observed and recorded by daily ultrasonography.Ovaries were removed from cows at the following stages of the first follicular wave: predeviation (PD; approximately Day 3 after estrus; 1.5 days after emergence [emergence is the first scan where a new follicle at least 4 mm is detected by ovarian ultrasonography]) and onset of deviation (OD; first scan where growth of the F1 [largest; future dominant] follicle to >8.5 mm was detected by ovarian ultrasonography). The F1 follicles were isolated from the PD and OD groups. Granulosa cells were isolated from the two types of follicles (PDF1 and ODF1), lysed, and stored at −80°C immediately.


### 2.3. RNA Isolation

Total RNA was isolated from the lysed granulosa cells using the RNeasy mini kit (Qiagen) and DNase treated on column according to the manufacturer’s protocol. The RNA integrity was evaluated by Agilient Bioanalyser and the RNA concentration was measured using a Nanodrop-1000 spectrophotometer. RNA samples with a RNA integrity number greater than 8 were selected for deep sequencing.


### 2.4. Library Preparation and Illumina Sequencing

RNA sequencing was performed by the WM Keck Center for Comparative and Functional Genomics at the University of Illinois at Urbana-Champaign. RNA samples from four follicles were pooled within each group (ODF1 or PDF1). RNA-Seq libraries were prepared with a TruSeq RNA Sample Preparation kit (Illumina) according to the manufacturer’s instructions. The cDNA libraries were sequenced on one lane for 100 cycles using Illumina HiSeq™ 2000 by a TruSeq SBS kit v5 (Illumina) and analyzed with pipeline version 1.8.


### 2.5. Identification of Differentially Expressed Genes and Pathway Analysis

The CLC Genomics Workbench (CLC bio, Aarhus, Denmark) was used to map the sequence reads to the bovine RefSeq database. The reads per kilobase per million reads (RPKM) values were calculated as the normalized transcript expression values [9]. A Z-test [10] was used to identify differentially expressed genes between ODF1 and PDF1 (FDR corrected p value <
0.05, RPKM cut-off > 0.5, and RPKM fold change >1.5) using the CLC genomics workbench. DAVID software (https://david.ncifcrf.gov/gene2gene.jsp) was used to perform GO annotations and KEGG pathway analysis for highly expressed (RPKM > 100) and differentially expressed genes.


## 2.1. Materials

All materials were obtained from Sigma-Aldrich (St. Louis, MO) unless otherwise stated.


## 2.2. Animal Model and Sample Collection

All animal procedures were approved by the Institutional Animal Care and Use Committee at Michigan State University. Estrus was synchronized in nonlactating Holstein dairy cows with two injections of prostaglandin F2α (PGF2α; Prostamate; IVX Animal Health, St. Joseph, MO) administered 14 days apart, and follicular growth was observed and recorded by daily ultrasonography.Ovaries were removed from cows at the following stages of the first follicular wave: predeviation (PD; approximately Day 3 after estrus; 1.5 days after emergence [emergence is the first scan where a new follicle at least 4 mm is detected by ovarian ultrasonography]) and onset of deviation (OD; first scan where growth of the F1 [largest; future dominant] follicle to >8.5 mm was detected by ovarian ultrasonography). The F1 follicles were isolated from the PD and OD groups. Granulosa cells were isolated from the two types of follicles (PDF1 and ODF1), lysed, and stored at −80°C immediately.


## 2.3. RNA Isolation

Total RNA was isolated from the lysed granulosa cells using the RNeasy mini kit (Qiagen) and DNase treated on column according to the manufacturer’s protocol. The RNA integrity was evaluated by Agilient Bioanalyser and the RNA concentration was measured using a Nanodrop-1000 spectrophotometer. RNA samples with a RNA integrity number greater than 8 were selected for deep sequencing.


## 2.4. Library Preparation and Illumina Sequencing

RNA sequencing was performed by the WM Keck Center for Comparative and Functional Genomics at the University of Illinois at Urbana-Champaign. RNA samples from four follicles were pooled within each group (ODF1 or PDF1). RNA-Seq libraries were prepared with a TruSeq RNA Sample Preparation kit (Illumina) according to the manufacturer’s instructions. The cDNA libraries were sequenced on one lane for 100 cycles using Illumina HiSeq™ 2000 by a TruSeq SBS kit v5 (Illumina) and analyzed with pipeline version 1.8.


## 2.5. Identification of Differentially Expressed Genes and Pathway Analysis

The CLC Genomics Workbench (CLC bio, Aarhus, Denmark) was used to map the sequence reads to the bovine RefSeq database. The reads per kilobase per million reads (RPKM) values were calculated as the normalized transcript expression values [9]. A Z-test [10] was used to identify differentially expressed genes between ODF1 and PDF1 (FDR corrected p value <
0.05, RPKM cut-off > 0.5, and RPKM fold change >1.5) using the CLC genomics workbench. DAVID software (https://david.ncifcrf.gov/gene2gene.jsp) was used to perform GO annotations and KEGG pathway analysis for highly expressed (RPKM > 100) and differentially expressed genes.


## 3. Results

### 3.1. Illumina Sequencing

To identify differentially expressed genes involved in bovine follicular development, Illumina sequencing was used on two libraries constructed from RNA isolated from ODF1 and PDF1 follicles. After filtering, a total of 44,082,301 and 43,708,132 clean reads were obtained from PDF1 and ODF1 libraries, respectively. The clean reads were mapped to the bovine RefSeq database (containing 35,325 annotated transcripts). Using a cut-off value of RPKM > 0.5, a total of 15,533 genes were identified in both types of follicles (Additional file 1: Table S1 in Supplementary Material available online athttp://dx.doi.org/10.1155/2016/3472748), among which 719 are considered to be highly expressed (RPKM cut-off > 100) in bovine follicles (Additional file 2: Table S2).


### 3.2. GO Functional Classification and KEGG Pathway Analysis of Highly Expressed Genes

The top 30 highly expressed genes in granulosa cells of bovine follicles at the predeviation and onset of deviation stages are shown in Table1. Many of them are known to be important for follicular growth and development, such as Serpin peptidase inhibitor clade E member 2 (SERPINE2), Inhibin alpha (INHA), Inhibin beta A (INHBA), and Follistatin (FST). GO functional classification of these highly expressed genes was performed using DAVID software. All 719 highly expressed genes can be assigned into 22 groups under three categories (biological process, 39%; cellular component, 44%; and molecular function, 17%) based on their putative functions (Figure 1). Many of the highly expressed genes are involved in metabolic process, multicellular organismal process, and binding. KEGG pathway analysis showed that the highly expressed genes are involved in 12 major pathways (Figure 2), of which the most significantly enriched genes are involved in ribosome pathway.Table 1
Top 30 highly expressed genes in ODF1 and PDF1 follicles.
Gene name
PDF1-RPKM
ODF1-RPKM
Glutathione S-transferase alpha 3
12445.32
17611.44
Serpin peptidase inhibitor clade E member 2
6816.847
11075.44
Inhibin alpha
6341.615
9658.832
Inhibin beta A
5824.462
9440.244
Serglycin
8808.27
7315.526
Follistatin
3435.171
4998.088
Cytochrome c oxidase subunit I-like
5017.207
3707.555
Unknown
7026.759
3115.129
Cytochrome c oxidase subunit III-like
3508.943
2440.522
Cytochrome c oxidase subunit I-like
3252.12
2372.119
Unknown
5741.943
2262.251
Milk fat globule-EGF factor 8 protein
2439.639
2077.699
Vimentin
2032.584
2050.732
Lysosomal protein transmembrane 4 beta
1539.453
2025.655
Gap junction protein alpha 1
1408.231
1975.614
Cytochrome P450 family 19 subfamily A polypeptide 1
656.8769
1944.887
Eukaryotic translation elongation factor 1 alpha 1
1864.128
1924.407
Low density lipoprotein receptor-related protein 8 apolipoprotein E receptor
1061.542
1897.644
Enolase 1
1005.447
1707.612
Heat shock protein 8
1615.28
1700.711
Ribosomal protein L18a
1559.355
1680.399
Glyceraldehydes 3 phosphate dehydrogenase
1342.963
1670.066
Ribosomal protein S27a
1576.131
1589.545
Ribosomal protein
1627.751
1566.558
ST3 beta-galactoside alpha-2,3-sialyltransferase 4
1096.954
1542.973
Ribosomal protein L4
1443.707
1502.241
Ribosomal protein S8
1563.222
1490.917
Tribbles homolog 2
1411.494
1461.935
Ribosomal protein S3A
1409.527
1450.511
Cytochrome P450, family 11, subfamily A, polypeptide 1
908.9737
1425.869Figure 1
Highly expressed genes GO analysis.Figure 2
Highly expressed genes KEGG pathways analysis.


### 3.3. Differentially Expressed Genes in ODF1 versus PDF1

Using RPKM cut-off > 0.5 and fold change cut-off > 1.5 at FDR correctedp value <
0.05, a total of 83 differentially expressed genes were identified, with 41 downregulated genes (Table 2) and 42 upregulated genes (Table 3) in ODF1 versus PDF1. To understand the functions of these differentially expressed genes, GO analysis was performed. The upregulated genes were categorized into 14 functional groups under 3 major GO classifications: biological process (35%), cellular component (30%), and molecular function (35%) (Table 4). Many of the differentially expressed genes are known to play a role in ovarian follicular development (Table 5). For example, serine protease 23 (PRSS23) is expressed in granulosa cells and may play a crucial role in follicular atresia, whereas serine protease 35 (PRSS35) is also expressed in granulosa cells and may be involved in ovulation and CL formation and regression. KEGG pathway analysis of the upregulated genes demonstrated that two important genes (CYP11A1 and CYP19A1) in the steroid hormone biosynthesis pathway are upregulated in ODF1 versus PDF1.Table 2
List of downregulated genes in ODF1 versus PDF1 and their functions.
Gene symbol
GenBank number
PDF1 RPKM
ODF1 RPKM
Fold change
FDR correctedp value
Gene product functions
ACTR1A
NM_001193248.1
17.31
0.43
−39.45
2.63 × 10−2
Vesicle motility
LOC787803
XM_002700116.1
45.66
2.15
−20.79
3.41 × 10−7
Unknown
PPP1R14A
XM_002694966.1
68.29
3.99
−16.76
7.75 × 10−11
Protein phosphatase inhibitor
OLA1
NM_001046045.1
22.71
1.74
−12.79
1.10 × 10−2
Hydrolase activity and GTP binding
QRFPR
NM_001192681.1
19.15
1.55
−12.07
4.40 × 10−2
Modulate adenylate cyclase
RMRP
NR_036646.1
149.79
13.57
−10.82
Lncrna class
RN5-8S1
NR_036643.1
9330.16
856.09
−10.68
Unknown
LOC100335749
XR_083021.1
35.26
4.58
−7.54
8.45 × 10−4
Senescence-associated protein-like
C11H2orf40
NM_001038113.1
41.02
5.84
−6.88
2.22 × 10−4
Esophageal cancer
BOLA
NM_001040532.1
74.93
11.19
−6.56
1.01 × 10−8
Transcription
ANGPT2
NM_001098855.1
39.04
7.91
−4.83
3.56 × 10−3
Angiogenic signal
VNN1
NM_001024556.2
33.51
7.26
−4.52
1.94 × 10−2
Amidohydrolase
KRT2
XM_001254015.1
40.64
9.41
−4.23
6.20 × 10−3
Keratinocyte activation
IHH
NM_001076870.2
42.84
9.98
−4.21
3.93 × 10−3
Smoothened
RN18S1
NR_036642.1
4431.16
1040.31
−4.17
Unknown
BOLA
NM_001038518.1
101.2
26.48
−3.74
6.16 × 10−8
Transcription
LOC100335409
XM_002705970.1
1012.71
275.19
−3.61
Unknown
LOC100140002
XR_084188.1
38.53
10.81
−3.49
3.63 × 10−2
Envelope glycoprotein-like
ADM
NM_173888.3
203.61
61
−3.27
Adrenomedullin
4-Sep
NM_001034651.1
69.77
22.91
−2.98
9.40 × 10−4
Cytokinesis, platelet secretion
ITPR1
NM_174841.2
61.79
22.18
−2.73
9.78 × 10−3
Intracellular channel
RN28S1
NR_036644.1
537.61
205.32
−2.56
Unknown
LOC100336997
XM_002685421.1
5741.94
2262.25
−2.49
Unknown
PRSS35
NM_001035457.3
134.05
55.07
−2.38
1.21 × 10−5
Ovulation, CL formation and regression
LOC100337434
XR_083937.1
7026.76
3115.13
−2.21
Unknown
LOC100140226
XM_001787664.2
1183.02
532.17
−2.18
Zinc finger protein 347-like
LOC100337402
XR_083935.1
908.3
413.36
−2.15
Unknown
LOC511901
XM_589328.5
102.29
49.83
−2.01
1.42 × 10−2
H1 histone
LOC100137883
XM_002706880.1
98.48
49.64
−1.94
3.47 × 10−2
Thymosin beta-4-like
CDH2
NM_001166492.1
103.68
53.3
−1.91
3.44 × 10−2
Neuronal recognition
APOA1
NM_174242.3
126.56
66.22
−1.87
1.07 × 10−2
Activates spermatozoa motility
PAPSS2
NM_001076075.1
119.75
65.23
−1.8
3.68 × 10−2
Skeletogenesis
LOC100299201
XR_084007.1
327.62
178.8
−1.79
1.28 × 10−7
Ribosomal protein
GSTA5
NM_001099016.1
138.14
80.02
−1.69
4.93 × 10−2
Glutathione transferase
AKR1B1
NM_001012519.1
182.74
108.97
−1.64
1.33 × 10−2
Electron carrier activity
SLCO1A2
NM_174654.2
192.47
119.34
−1.58
2.65 × 10−2
Mediates transport
HERC1
NM_001103282.1
222.29
138.46
−1.57
1.01 × 10−2
Membrane trafficking
CWC25
NM_001105359.1
433.81
274.72
−1.55
7.12 × 10−6
Alternatively spliced transcripts
ACOT11
NM_001103275.1
693.81
440.12
−1.54
6.15 × 10−10
Acyl-Coa thioesterase activity
LOC615589
NM_001098467.1
245.33
157.07
−1.53
1.10 × 10−2
Keratin-like protein
C12H13orf18
NM_001102041.1
269.27
173.55
−1.52
6.34 × 10−3
UnknownTable 3
List of upregulated genes in ODF1 versus PDF1 and their functions.
Gene symbol
GenBank number
PDF1 RPKM
ODF1 RPKM
Fold change
FDR correctedp value
Gene product functions
GAPDH
XM_001252511.3
1.91
59.65
31.91
9.83 × 10−11
Microtubule and NAD binding
BOLA-N
NM_001105651.1
14.28
174.88
12.50
Unknown
PPP1R14A
NM_001193070.1
6.06
65.75
11.08
8.98 × 10−10
Smooth muscle contraction
LOC505676
NM_001193296.1
14.04
111.19
8.09
Unknown
LOC100337308
XM_002684003.1
3.72
23.27
6.39
4.40 × 10−2
Unknown
MT1A
NM_001040492.2
26.50
150.75
5.81
Bind heavy metals
LOC100125916
NM_001105487.1
14.80
83.88
5.78
1.51 × 10−9
Unknown
TNFAIP6
NM_001007813.1
11.17
46.61
4.26
1.06 × 10−3
Cell-cell and cell-matrix interactions
BEX2
NM_001077034.1
93.97
345.00
3.75
Mitochondrial apoptosis
GPR85
NM_001075150.2
10.06
35.54
3.61
4.05 × 10−2
G-protein coupled receptor
PPM1K
NM_001046474.1
34.75
108.32
3.18
2.52 × 10−7
Cellular survival and development
CYP19A1
NM_174305.1
656.88
1944.89
3.02
Estrogen biosynthesis
MT1E
NM_001114857.1
54.11
152.16
2.87
2.61 × 10−9
Bind heavy metals
ETNK2
XM_002693881.1
22.35
61.96
2.83
4.65 × 10−3
Ethanolamine phosphorylation
CHST11
NM_001192668.1
57.39
154.72
2.75
6.77 × 10−9
Biosynthesis chondroitin sulfate
MT2A
NM_001075140.1
49.37
130.49
2.70
4.16 × 10−7
Bind heavy metals
PRSS23
NM_001080306.1
58.01
151.79
2.67
2.74 × 10−8
Follicular atresia
TXNIP
NM_001101875.2
98.64
231.47
2.40
6.64 × 10−11
Oxidative stress mediator
NPR3
NM_174127.2
47.55
106.35
2.28
5.93 × 10−4
Natriuretic peptide hormone receptor
GREB1
NM_001205631.1
35.17
78.20
2.27
1.30 × 10−2
Estrogen-stimulated cell proliferation
EIF4EBP1
NM_001077893.1
138.11
294.07
2.17
1.30 × 10−11
Mediates protein translation regulation
PIK3R1
NM_174575.1
113.42
210.95
1.90
1.21 × 10−5
Insulin actions metabolic
LRP8
NM_001097565.1
1061.54
1897.64
1.82
Sperm maturation
SCD5
NM_001076945.1
124.08
210.97
1.74
4.25 × 10−4
Energy metabolism
ENO1
NM_174049.2
1005.45
1707.61
1.73
Tumor suppressor
LDHA
NM_174099.2
221.85
366.73
1.69
3.41 × 10−7
Affiliated with lncrna
SERPINE2
NM_174669.2
6816.85
11075.44
1.66
Serine protease inhibitor
INHBA
NM_174363.2
5824.46
9440.24
1.65
Regulate gonadal stromal cell proliferation
TMEM176B
NM_001099145.1
106.45
170.39
1.63
2.16 × 10−2
Dendritic cells maturation
OAT
NM_001034240.1
395.42
630.73
1.63
1.71 × 10−11
Ornithine aminotransferase
CYP11A1
NM_176644.2
908.97
1425.87
1.60
Cholesterol to pregnenolone
OBSL1
XM_002685586.1
221.33
338.40
1.56
1.16 × 10−4
Regulate ubiquitin ligase complex
ARFGAP3
NM_001075974.1
273.92
418.22
1.56
5.26 × 10−6
Gtpase-activating protein
ITGB5
NM_174679.2
153.24
233.56
1.56
7.68 × 10−3
Fibronectin receptor
INHA
NM_174094.3
6341.61
9658.83
1.55
Gonadal hormone secretion
OPTN
NM_001034602.1
160.13
243.40
1.55
5.92 × 10−3
Affect cell death
PTGR1
NM_001035281.1
200.07
301.66
1.54
9.40 × 10−4
Inactivation of the chemotactic factor
STBD1
XM_002688357.1
143.24
215.66
1.54
2.15 × 10−2
Bind to carbohydrates
LOC532189
XR_083049.1
340.18
510.18
1.53
5.69 × 10−7
Carboxypeptidase
TMEM20
NM_001076470.1
191.10
286.48
1.53
2.10 × 10−3
Solute carrier
ECE1
NM_181009.2
461.24
691.33
1.53
8.98 × 10−10
Converts big endothelin-1 to endothelin-1
GNG10
NM_001114512.1
617.73
918.26
1.52
7.59 × 10−13
Signal transducerTable 4
GO analysis of genes upregulated in ODF1 versus PDF1.
Category
Term
Genes
Biological process
Hemoglobin biosynthetic process
INHBA, INHA
Biological process
Antigen processing and presentation
LOC505676, LOC100125916, BOLA-N
Biological process
Insulin receptor signaling pathway
EIF4EBP1, PIK3R1
Biological process
Antigen processing and presentation of peptide antigen
LOC100125916, BOLA-N
Biological process
Response to insulin stimulus
EIF4EBP1, PIK3R1
Biological process
Cellular hormone metabolic process
ECE1, CYP11A1
Biological process
Cellular response to hormone stimulus
EIF4EBP1, PIK3R1
Biological process
Regulation of myeloid cell differentiation
INHBA, PIK3R1
Cellular component
Inhibin complex
INHBA, INHA
Cellular component
MHC protein complex
LOC505676, LOC100125916, BOLA-N
Cellular component
Plasma membrane
LOC505676, ARFGAP3, ECE1, GNG10, ITGB5, INHA, LOC100125916, BOLA-N, ENO1
Molecular function
Metal ion binding
ARFGAP3, PTGR1, ECE1, MT1A, CYP11A1, PPM1K, LRP8, SCD5, CYP19A1, ENO1
Molecular function
Peptidase activity
ECE1, SERPINE2, PRSS23, ENO1
Molecular function
Iron ion binding
CYP11A1, SCD5, CYP19A1Table 5
Differentially expressed genes between ODF1 and PDF1 with known or predicted roles in follicular development.
Gene symbol
Fold change
FDR correctedp value
Up- or downregulation in ODF1/PDF1
Gene product functions
PPM1K
3.18
2.52 × 10−7
Upregulation
Cellular survival and development
BEX2
3.75
Upregulation
Mitochondrial apoptosis
CYP19A1
3.02
Upregulation
Estrogen biosynthesis
PRSS23
2.67
2.74 × 10−8
Upregulation
Follicular atresia
GREB1
2.27
1.30 × 10−2
Upregulation
Estrogen-stimulated cell proliferation
SERPINE2
1.66
Upregulation
Serine protease inhibitor
INHBA
1.65
Upregulation
Regulate gonadal stromal cell proliferation
CYP11A1
1.60
Upregulation
Cholesterol to pregnenolone
INHA
1.55
Upregulation
Gonadal hormone secretion
TNFAIP6
4.26
1.06 × 10−3
Upregulation
Cell-cell and cell-matrix interactions
OPTN
1.55
5.92 × 10−3
Upregulation
Affect cell death
PRSS35
−2.38
1.21 × 10−5
Downregulation
Ovulation and CL formation and regression
APOA1
−1.87
1.07 × 10−2
Downregulation
Activates spermatozoa motility
GSTA5
−1.69
4.93 × 10−2
Downregulation
Glutathione transferase


## 3.1. Illumina Sequencing

To identify differentially expressed genes involved in bovine follicular development, Illumina sequencing was used on two libraries constructed from RNA isolated from ODF1 and PDF1 follicles. After filtering, a total of 44,082,301 and 43,708,132 clean reads were obtained from PDF1 and ODF1 libraries, respectively. The clean reads were mapped to the bovine RefSeq database (containing 35,325 annotated transcripts). Using a cut-off value of RPKM > 0.5, a total of 15,533 genes were identified in both types of follicles (Additional file 1: Table S1 in Supplementary Material available online athttp://dx.doi.org/10.1155/2016/3472748), among which 719 are considered to be highly expressed (RPKM cut-off > 100) in bovine follicles (Additional file 2: Table S2).


## 3.2. GO Functional Classification and KEGG Pathway Analysis of Highly Expressed Genes

The top 30 highly expressed genes in granulosa cells of bovine follicles at the predeviation and onset of deviation stages are shown in Table1. Many of them are known to be important for follicular growth and development, such as Serpin peptidase inhibitor clade E member 2 (SERPINE2), Inhibin alpha (INHA), Inhibin beta A (INHBA), and Follistatin (FST). GO functional classification of these highly expressed genes was performed using DAVID software. All 719 highly expressed genes can be assigned into 22 groups under three categories (biological process, 39%; cellular component, 44%; and molecular function, 17%) based on their putative functions (Figure 1). Many of the highly expressed genes are involved in metabolic process, multicellular organismal process, and binding. KEGG pathway analysis showed that the highly expressed genes are involved in 12 major pathways (Figure 2), of which the most significantly enriched genes are involved in ribosome pathway.Table 1
Top 30 highly expressed genes in ODF1 and PDF1 follicles.
Gene name
PDF1-RPKM
ODF1-RPKM
Glutathione S-transferase alpha 3
12445.32
17611.44
Serpin peptidase inhibitor clade E member 2
6816.847
11075.44
Inhibin alpha
6341.615
9658.832
Inhibin beta A
5824.462
9440.244
Serglycin
8808.27
7315.526
Follistatin
3435.171
4998.088
Cytochrome c oxidase subunit I-like
5017.207
3707.555
Unknown
7026.759
3115.129
Cytochrome c oxidase subunit III-like
3508.943
2440.522
Cytochrome c oxidase subunit I-like
3252.12
2372.119
Unknown
5741.943
2262.251
Milk fat globule-EGF factor 8 protein
2439.639
2077.699
Vimentin
2032.584
2050.732
Lysosomal protein transmembrane 4 beta
1539.453
2025.655
Gap junction protein alpha 1
1408.231
1975.614
Cytochrome P450 family 19 subfamily A polypeptide 1
656.8769
1944.887
Eukaryotic translation elongation factor 1 alpha 1
1864.128
1924.407
Low density lipoprotein receptor-related protein 8 apolipoprotein E receptor
1061.542
1897.644
Enolase 1
1005.447
1707.612
Heat shock protein 8
1615.28
1700.711
Ribosomal protein L18a
1559.355
1680.399
Glyceraldehydes 3 phosphate dehydrogenase
1342.963
1670.066
Ribosomal protein S27a
1576.131
1589.545
Ribosomal protein
1627.751
1566.558
ST3 beta-galactoside alpha-2,3-sialyltransferase 4
1096.954
1542.973
Ribosomal protein L4
1443.707
1502.241
Ribosomal protein S8
1563.222
1490.917
Tribbles homolog 2
1411.494
1461.935
Ribosomal protein S3A
1409.527
1450.511
Cytochrome P450, family 11, subfamily A, polypeptide 1
908.9737
1425.869Figure 1
Highly expressed genes GO analysis.Figure 2
Highly expressed genes KEGG pathways analysis.


## 3.3. Differentially Expressed Genes in ODF1 versus PDF1

Using RPKM cut-off > 0.5 and fold change cut-off > 1.5 at FDR correctedp value <
0.05, a total of 83 differentially expressed genes were identified, with 41 downregulated genes (Table 2) and 42 upregulated genes (Table 3) in ODF1 versus PDF1. To understand the functions of these differentially expressed genes, GO analysis was performed. The upregulated genes were categorized into 14 functional groups under 3 major GO classifications: biological process (35%), cellular component (30%), and molecular function (35%) (Table 4). Many of the differentially expressed genes are known to play a role in ovarian follicular development (Table 5). For example, serine protease 23 (PRSS23) is expressed in granulosa cells and may play a crucial role in follicular atresia, whereas serine protease 35 (PRSS35) is also expressed in granulosa cells and may be involved in ovulation and CL formation and regression. KEGG pathway analysis of the upregulated genes demonstrated that two important genes (CYP11A1 and CYP19A1) in the steroid hormone biosynthesis pathway are upregulated in ODF1 versus PDF1.Table 2
List of downregulated genes in ODF1 versus PDF1 and their functions.
Gene symbol
GenBank number
PDF1 RPKM
ODF1 RPKM
Fold change
FDR correctedp value
Gene product functions
ACTR1A
NM_001193248.1
17.31
0.43
−39.45
2.63 × 10−2
Vesicle motility
LOC787803
XM_002700116.1
45.66
2.15
−20.79
3.41 × 10−7
Unknown
PPP1R14A
XM_002694966.1
68.29
3.99
−16.76
7.75 × 10−11
Protein phosphatase inhibitor
OLA1
NM_001046045.1
22.71
1.74
−12.79
1.10 × 10−2
Hydrolase activity and GTP binding
QRFPR
NM_001192681.1
19.15
1.55
−12.07
4.40 × 10−2
Modulate adenylate cyclase
RMRP
NR_036646.1
149.79
13.57
−10.82
Lncrna class
RN5-8S1
NR_036643.1
9330.16
856.09
−10.68
Unknown
LOC100335749
XR_083021.1
35.26
4.58
−7.54
8.45 × 10−4
Senescence-associated protein-like
C11H2orf40
NM_001038113.1
41.02
5.84
−6.88
2.22 × 10−4
Esophageal cancer
BOLA
NM_001040532.1
74.93
11.19
−6.56
1.01 × 10−8
Transcription
ANGPT2
NM_001098855.1
39.04
7.91
−4.83
3.56 × 10−3
Angiogenic signal
VNN1
NM_001024556.2
33.51
7.26
−4.52
1.94 × 10−2
Amidohydrolase
KRT2
XM_001254015.1
40.64
9.41
−4.23
6.20 × 10−3
Keratinocyte activation
IHH
NM_001076870.2
42.84
9.98
−4.21
3.93 × 10−3
Smoothened
RN18S1
NR_036642.1
4431.16
1040.31
−4.17
Unknown
BOLA
NM_001038518.1
101.2
26.48
−3.74
6.16 × 10−8
Transcription
LOC100335409
XM_002705970.1
1012.71
275.19
−3.61
Unknown
LOC100140002
XR_084188.1
38.53
10.81
−3.49
3.63 × 10−2
Envelope glycoprotein-like
ADM
NM_173888.3
203.61
61
−3.27
Adrenomedullin
4-Sep
NM_001034651.1
69.77
22.91
−2.98
9.40 × 10−4
Cytokinesis, platelet secretion
ITPR1
NM_174841.2
61.79
22.18
−2.73
9.78 × 10−3
Intracellular channel
RN28S1
NR_036644.1
537.61
205.32
−2.56
Unknown
LOC100336997
XM_002685421.1
5741.94
2262.25
−2.49
Unknown
PRSS35
NM_001035457.3
134.05
55.07
−2.38
1.21 × 10−5
Ovulation, CL formation and regression
LOC100337434
XR_083937.1
7026.76
3115.13
−2.21
Unknown
LOC100140226
XM_001787664.2
1183.02
532.17
−2.18
Zinc finger protein 347-like
LOC100337402
XR_083935.1
908.3
413.36
−2.15
Unknown
LOC511901
XM_589328.5
102.29
49.83
−2.01
1.42 × 10−2
H1 histone
LOC100137883
XM_002706880.1
98.48
49.64
−1.94
3.47 × 10−2
Thymosin beta-4-like
CDH2
NM_001166492.1
103.68
53.3
−1.91
3.44 × 10−2
Neuronal recognition
APOA1
NM_174242.3
126.56
66.22
−1.87
1.07 × 10−2
Activates spermatozoa motility
PAPSS2
NM_001076075.1
119.75
65.23
−1.8
3.68 × 10−2
Skeletogenesis
LOC100299201
XR_084007.1
327.62
178.8
−1.79
1.28 × 10−7
Ribosomal protein
GSTA5
NM_001099016.1
138.14
80.02
−1.69
4.93 × 10−2
Glutathione transferase
AKR1B1
NM_001012519.1
182.74
108.97
−1.64
1.33 × 10−2
Electron carrier activity
SLCO1A2
NM_174654.2
192.47
119.34
−1.58
2.65 × 10−2
Mediates transport
HERC1
NM_001103282.1
222.29
138.46
−1.57
1.01 × 10−2
Membrane trafficking
CWC25
NM_001105359.1
433.81
274.72
−1.55
7.12 × 10−6
Alternatively spliced transcripts
ACOT11
NM_001103275.1
693.81
440.12
−1.54
6.15 × 10−10
Acyl-Coa thioesterase activity
LOC615589
NM_001098467.1
245.33
157.07
−1.53
1.10 × 10−2
Keratin-like protein
C12H13orf18
NM_001102041.1
269.27
173.55
−1.52
6.34 × 10−3
UnknownTable 3
List of upregulated genes in ODF1 versus PDF1 and their functions.
Gene symbol
GenBank number
PDF1 RPKM
ODF1 RPKM
Fold change
FDR correctedp value
Gene product functions
GAPDH
XM_001252511.3
1.91
59.65
31.91
9.83 × 10−11
Microtubule and NAD binding
BOLA-N
NM_001105651.1
14.28
174.88
12.50
Unknown
PPP1R14A
NM_001193070.1
6.06
65.75
11.08
8.98 × 10−10
Smooth muscle contraction
LOC505676
NM_001193296.1
14.04
111.19
8.09
Unknown
LOC100337308
XM_002684003.1
3.72
23.27
6.39
4.40 × 10−2
Unknown
MT1A
NM_001040492.2
26.50
150.75
5.81
Bind heavy metals
LOC100125916
NM_001105487.1
14.80
83.88
5.78
1.51 × 10−9
Unknown
TNFAIP6
NM_001007813.1
11.17
46.61
4.26
1.06 × 10−3
Cell-cell and cell-matrix interactions
BEX2
NM_001077034.1
93.97
345.00
3.75
Mitochondrial apoptosis
GPR85
NM_001075150.2
10.06
35.54
3.61
4.05 × 10−2
G-protein coupled receptor
PPM1K
NM_001046474.1
34.75
108.32
3.18
2.52 × 10−7
Cellular survival and development
CYP19A1
NM_174305.1
656.88
1944.89
3.02
Estrogen biosynthesis
MT1E
NM_001114857.1
54.11
152.16
2.87
2.61 × 10−9
Bind heavy metals
ETNK2
XM_002693881.1
22.35
61.96
2.83
4.65 × 10−3
Ethanolamine phosphorylation
CHST11
NM_001192668.1
57.39
154.72
2.75
6.77 × 10−9
Biosynthesis chondroitin sulfate
MT2A
NM_001075140.1
49.37
130.49
2.70
4.16 × 10−7
Bind heavy metals
PRSS23
NM_001080306.1
58.01
151.79
2.67
2.74 × 10−8
Follicular atresia
TXNIP
NM_001101875.2
98.64
231.47
2.40
6.64 × 10−11
Oxidative stress mediator
NPR3
NM_174127.2
47.55
106.35
2.28
5.93 × 10−4
Natriuretic peptide hormone receptor
GREB1
NM_001205631.1
35.17
78.20
2.27
1.30 × 10−2
Estrogen-stimulated cell proliferation
EIF4EBP1
NM_001077893.1
138.11
294.07
2.17
1.30 × 10−11
Mediates protein translation regulation
PIK3R1
NM_174575.1
113.42
210.95
1.90
1.21 × 10−5
Insulin actions metabolic
LRP8
NM_001097565.1
1061.54
1897.64
1.82
Sperm maturation
SCD5
NM_001076945.1
124.08
210.97
1.74
4.25 × 10−4
Energy metabolism
ENO1
NM_174049.2
1005.45
1707.61
1.73
Tumor suppressor
LDHA
NM_174099.2
221.85
366.73
1.69
3.41 × 10−7
Affiliated with lncrna
SERPINE2
NM_174669.2
6816.85
11075.44
1.66
Serine protease inhibitor
INHBA
NM_174363.2
5824.46
9440.24
1.65
Regulate gonadal stromal cell proliferation
TMEM176B
NM_001099145.1
106.45
170.39
1.63
2.16 × 10−2
Dendritic cells maturation
OAT
NM_001034240.1
395.42
630.73
1.63
1.71 × 10−11
Ornithine aminotransferase
CYP11A1
NM_176644.2
908.97
1425.87
1.60
Cholesterol to pregnenolone
OBSL1
XM_002685586.1
221.33
338.40
1.56
1.16 × 10−4
Regulate ubiquitin ligase complex
ARFGAP3
NM_001075974.1
273.92
418.22
1.56
5.26 × 10−6
Gtpase-activating protein
ITGB5
NM_174679.2
153.24
233.56
1.56
7.68 × 10−3
Fibronectin receptor
INHA
NM_174094.3
6341.61
9658.83
1.55
Gonadal hormone secretion
OPTN
NM_001034602.1
160.13
243.40
1.55
5.92 × 10−3
Affect cell death
PTGR1
NM_001035281.1
200.07
301.66
1.54
9.40 × 10−4
Inactivation of the chemotactic factor
STBD1
XM_002688357.1
143.24
215.66
1.54
2.15 × 10−2
Bind to carbohydrates
LOC532189
XR_083049.1
340.18
510.18
1.53
5.69 × 10−7
Carboxypeptidase
TMEM20
NM_001076470.1
191.10
286.48
1.53
2.10 × 10−3
Solute carrier
ECE1
NM_181009.2
461.24
691.33
1.53
8.98 × 10−10
Converts big endothelin-1 to endothelin-1
GNG10
NM_001114512.1
617.73
918.26
1.52
7.59 × 10−13
Signal transducerTable 4
GO analysis of genes upregulated in ODF1 versus PDF1.
Category
Term
Genes
Biological process
Hemoglobin biosynthetic process
INHBA, INHA
Biological process
Antigen processing and presentation
LOC505676, LOC100125916, BOLA-N
Biological process
Insulin receptor signaling pathway
EIF4EBP1, PIK3R1
Biological process
Antigen processing and presentation of peptide antigen
LOC100125916, BOLA-N
Biological process
Response to insulin stimulus
EIF4EBP1, PIK3R1
Biological process
Cellular hormone metabolic process
ECE1, CYP11A1
Biological process
Cellular response to hormone stimulus
EIF4EBP1, PIK3R1
Biological process
Regulation of myeloid cell differentiation
INHBA, PIK3R1
Cellular component
Inhibin complex
INHBA, INHA
Cellular component
MHC protein complex
LOC505676, LOC100125916, BOLA-N
Cellular component
Plasma membrane
LOC505676, ARFGAP3, ECE1, GNG10, ITGB5, INHA, LOC100125916, BOLA-N, ENO1
Molecular function
Metal ion binding
ARFGAP3, PTGR1, ECE1, MT1A, CYP11A1, PPM1K, LRP8, SCD5, CYP19A1, ENO1
Molecular function
Peptidase activity
ECE1, SERPINE2, PRSS23, ENO1
Molecular function
Iron ion binding
CYP11A1, SCD5, CYP19A1Table 5
Differentially expressed genes between ODF1 and PDF1 with known or predicted roles in follicular development.
Gene symbol
Fold change
FDR correctedp value
Up- or downregulation in ODF1/PDF1
Gene product functions
PPM1K
3.18
2.52 × 10−7
Upregulation
Cellular survival and development
BEX2
3.75
Upregulation
Mitochondrial apoptosis
CYP19A1
3.02
Upregulation
Estrogen biosynthesis
PRSS23
2.67
2.74 × 10−8
Upregulation
Follicular atresia
GREB1
2.27
1.30 × 10−2
Upregulation
Estrogen-stimulated cell proliferation
SERPINE2
1.66
Upregulation
Serine protease inhibitor
INHBA
1.65
Upregulation
Regulate gonadal stromal cell proliferation
CYP11A1
1.60
Upregulation
Cholesterol to pregnenolone
INHA
1.55
Upregulation
Gonadal hormone secretion
TNFAIP6
4.26
1.06 × 10−3
Upregulation
Cell-cell and cell-matrix interactions
OPTN
1.55
5.92 × 10−3
Upregulation
Affect cell death
PRSS35
−2.38
1.21 × 10−5
Downregulation
Ovulation and CL formation and regression
APOA1
−1.87
1.07 × 10−2
Downregulation
Activates spermatozoa motility
GSTA5
−1.69
4.93 × 10−2
Downregulation
Glutathione transferase


## 4. Discussion

Follicular growth occurs in a characteristic wave-like pattern in monotocous species such as cattle [3, 5, 11]. A transient increase in FSH triggers initiation of each follicular wave [5, 11, 12]. Emergence is defined as the first day a new follicle >4 mm in diameter is detected and is the first chronological event marking a new follicular wave that is detectable by ultrasonography. After emergence, follicles in the cohort initially grow at a similar rate (common growth phase) prior to deviation [3]. However, the molecular mechanisms involved regulating the onset of deviation are not well understood, in order to characterize the differences in gene expression that associated with follicular development in different follicles sized in diameter, which the previous studies examined using microarray technology [13–17]. To further investigate the bovine granulosa cell transcriptome and molecular alterations associated with onset of deviation, we examined the transcriptome at specific stages of the estrous cycle.Illumina sequencing technology was used to determine gene expression levels in ODF1 and PDF1 follicles. A total of 15,533 genes were identified in both types of follicles and 83 of them were identified as differentially expressed between ODF1 and PDF1. Our study provided novel information on the bovine granulosa cell transcriptome and identified specific transcripts highly expressed in granulosa cells of bovine follicles prior to and at onset of deviation, including transcripts encoding for several housekeeping genes (e.g., ribosomal proteins L18a, S27a, and L4) and genes with well-established roles in regulation of ovarian function (e.g., INHBA, INHBB, and FST). Of particular interest was SERPINE2, which is abundantly expressed in granulosa cells of follicles collected at both the predeviation and onset of deviation stages of a follicular wave, illustrating its potential importance in bovine ovarian follicular development. Estradiol and SERPINE2 secretion are positively correlated, but estradiol treatment cannot alter the expression of SERPINE2. FSH and growth factors can directly regulate the expression and secretion of SERPINE2 in granulosa cells, and SERPINE2 is an antiapoptotic factor, which may regulate atresia in bovine follicles [18]. Eleven SERPINE genes are expressed in bovine follicles, but only SERPINE2, SERPINE5, and SERPINE6 are expressed in the granulosa cells [19].KEGG analysis revealed upregulated genes associated with onset of deviation (CYP11A1 and CYP19A1) involved in the steroid hormone biosynthesis pathways that play an essential role during follicular development. Proteins encoded by CYP11A1 and CYP19A1 genes are members of the cytochrome P450 superfamily, which are monooxygenases that catalyze many reactions involved in steroidogenesis. Previous studies suggested that CYP19A1 was regulated by multiple pathways, including estrogen receptors and cAMP/protein kinase A which are activated by FSHR in granulosa cells, and these regulatory mechanisms are likely critical for acquisition of follicular dominance in cattle [20]. Our transcriptome sequencing data is consistent with these results. The increase in transcript abundance for CYP19A1 in ODF1 versus PDF1 follicles is consistent with the increase in estradiol producing capacity associated with diameter deviation [3].It is acknowledged that study design was not optimal due to limited biological replication because single pooled samples (n
=
4 per group) were used in Illumina sequencing analysis. Despite such limitations, results have significantly enhanced understanding of bovine follicle transcriptome composition and potential differences in gene expression associated with follicular development that are foundational to further study in the future; several interesting candidates were revealed for future investigation, particularly genes linked to regulation of cell proliferation and survival. For example, results of present studies suggest that PPM1K, a Mn2+/Mg2+-dependent protein phosphatase of PPM family, is potentially upregulated with the onset of deviation in the granulosa cell layer of bovine follicles. This protein is critical for cell survival and embryonic development and can regulate the mitochondrial membrane permeability transition pore opening [21]. Potential upregulation of granulosa cell BEX2 and GREB1 transcript abundance was also noted in association with onset of dominance and may be associated with enhanced granulosa cell survival. BEX2 can downregulate apoptosis and activate the JNK (Jun NH2-terminal kinase) pathway, and these effects can be abolished by administration of a JNK specific inhibitor [22]. GREB1 is an estrogen receptor and coactivator linked to cell proliferation and GREB1 expression is estrogen dependent. It is possible that increased expression of PPM1K, BEX2, and GREB1 may be associated with granulosa cell proliferation and survival during the onset of deviation [23]. TNFAIP6/TSG6 is tumor necrosis factor and alpha-induced protein 6; it is suggested that TSG-6 plays a role in cell-cell or cell-cell matrix interactions during inflammation and tumorigenesis. High LH/hCGR gene expression intensity was associated with TNFAIP6/TSG6 gene expression which has a pivotal importance in the mucification of the COC during the preovulatory period [24]. It is suggested that the expression levels of TNFAIP6/TSG6 were nearly 280-fold in granulosa cells of large follicles than that of small follicles [14]. In our study, TNFAIP6/TSG6 was also differentially expressed in ODF1 and PDF1 with a 4.26-fold change. Altogether, these characteristics suggest that TNFAIP6/TSG6 plays a crucial role in accelerating follicle growth during follicular waves in cattle.


## 5. Conclusions

The present study characterized the granulosa cell transcriptome of bovine follicles at specific stages of follicular development and identified 83 differentially expressed genes potentially associated with onset of deviation, many of which are linked to regulation of follicular development. The study provides a foundation for future studies to investigate regulation of granulosa cell expressed genes and the regulatory mechanisms controlling antral follicle development during follicular waves in cattle.



---

*Source: 3472748-2016-03-21.xml*