Transcriptomic and proteomic analyses between small yellow follicles and the smallest hierarchical follicles reveal a role of VLDLR in chicken follicle selection

DOI: https://doi.org/10.21203/rs.2.17247/v2

Abstract

Background Follicle selection in chicken refers to the process of selecting one from a group of small yellow follicles (SY, 6-8mm in diameter) to enter the 12-15 mm hierarchical follicles (usually F6 follicles), which is a an important process affecting laying performance in the poultry industry. Although transcripromic analysis on chicken ovarian follicles was reported, integrated analysis on chicken follicles around selection by using both transcripromic and proteomic approaches was still lacking. In this study, we compared the proteomes and transcriptomes of SY and F6 follicles of laying hens and found some genes involved in chicken follicle selection.

Results Transcriptomic analysis revealed 855 differentially expressed genes (DEGs) between SY follicles and F6 follicles of laying hens, among which 202 were upregulated and 653 were downregulated. Proteomic analysis revealed 259 differentially expressed proteins (DEPs), including 175 upregulated and 84 downregulated proteins. Among the identified DEGs and DEPs, the expression changes of seven genes including VLDLR1,WIF1, NGFR, AMH, BMP15, GDF6 and MMP13 , and nine proteins including VLDLR, VTG1, VTG3, PSCA, APOB, APOV1, F10, ZP2 and ZP3L2 were validated. In addition VLDLR expression was significantly down-regulated in F6 follicles compared with SY follicles, was signifcantly higher in the GCs than in the TCs and was stimulated by FSH in GCs of both hierarchical and prehierarchical follicles.

Conclusions By comparing the proteomes and transcriptomes of SY follicles and F6 of laying hens, we identified some differentially expressed proteins/genes that might play certain roles in chicken follicle selection. These data may contribute to identification of the functional genes and proteins involved in chicken follicular development and selection.

Background

The ovary is a dynamic organ and a pivotal component of the reproductive system in hens. In the abdomen of laying hens, ovarian follicles of various sizes exist, including small white follicles less than 3 mm in diameter, large white follicles 3-5 mm in diameter, small yellow (SY) follicles 6-8 mm in diameter, large yellow follicles 9-12 mm in diameter and five to six hierarchical follicles with increased sizes, i.e. F6 to F1 [1]. In a reproductively active domestic hen, the hierarchical follicles are recruited daily from a pool of 6-8 mm prehierarchical small yellow follicles, which is called follicle selection [2]. The hierarchical follicles then continue to develop rapidly from the F6 follicle to F1 follicle until ovulation [3]. In the process of follicle selection, the SY follicle being selected absorbs egg yolk at accelerating rate, its granulosa cells rapidly proliferate and differentiate to produce high levels of progesterone [4, 5]. The absorption of the egg yolk by the SY follicles is mediated by receptors, and vitellin synthesized by the liver enters the oocytes via the very low-density lipoprotein receptor (VLDLR) [6, 7].

Transcriptome analysis of chicken prehierarchical follicles revealed changes in the transcripts involved in steroidogenesis, paracrine signaling and transcription during the early stage of follicular growth and development [8]. Comparison of chicken small white, F1 and post-ovulatory follicles transcriptomes identified differentially expressed genes involved in the adhesion, apoptosis and steroid biosynthesis pathways [9]. Candidate genes including ANXA2, Wnt4 and transforming growth factors play some roles in chicken follicle selection [10-15]. However, the dynamics of the transcriptome during chicken follicular selection are unclear. In addition, the mRNA abundance may not accurately predict the quantity of the corresponding functional proteins, while the proteomic approach can provide a systemic overview of protein levels [16]; therefore, has certain advantages over mRNA expression profiling [17].

Proteomic analyses on the ovarian function [18], maturation of oocytes [19], and early embryonic development [20-22] in mammals and in human ovarian diseases, such as polycystic ovarian syndrome (PCOS) and cancer [23] were reported, while in chicken, 2889 proteins were identified in the white yolk and ovarian stroma of small white follicles in Bovan's white laying hen [24]. However, the temporal changes in the proteome during chicken follicular selection are unknown. In this study, we compared the proteomes and transcriptomes of 6-8 mm SY follicles and the smallest hierarchical follicles 12-15 mm in diameter in laying hens and found some differentially expressed proteins/genes (DEPs/DEGs) that might play certain roles in chicken follicle selection.

Results

Transcriptomic analysis

RNA-seq was used to compare the transcriptomes of three SY follicles and three smallest heirachical follicles, which are referred to here as S1, S2, S3, and F1, F2 and F3, respectively. High-throughout RNA-seq generated 61.66 Gb clean data for six samples of chicken follicles, and 91.07%-93.42% reads could be mapped to the chicken genome. All six samples had at least 93.55% reads equal to or exceeding Q30 (Table 1).

A total of 855 DEGs, including 202 upregulated and 653 downregulated genes, were identified between the SY follicles (S) and F6 follicles (F) at the significant creteria of (|log2 (FoldChange)| >1 and padj <0.05) (Fig. 1a). A hierarchical clustered map of DEGs was then constructed and is shown in Fig 1b. Detailed analysis of the top 10 up-/down-regulated DEGs is shown in Table 2. The entire list of DEGs is shown in Additional file 1: Table S1.

DEGs were then assessed by GO and KEGG pathway analyses. Go functional analysis revealed that most of the DEGs were assigned to circulatory system process, cell differentiation and transition metal ion binding (Fig.1c). KEGG pathway analysis of DEGs showed that the most enriched pathways are TGF-β signaling pathway, tyrosine metabolism and cytokine-cytokine receptor interaction (Fig. 1d).

To validate the RNA-seq data, seven DEGs (Table 3) including very low density lipoprotein receptor 1 (VLDLR1), nerve growth factor receptor (NGFR), WNT inhibitory factor 1 (WIF1), anti-Mullerian hormone (AMH), bone morphogenetic protein 15 (BMP15), growth differentiation factor 6 (GDF6) and matrix metallopeptidase 13 (MMP13) were chosen and quantified by qRT-PCR. The result showed that the mRNA levels of these genes were similar to the sequencing data, suggesting that the RNA-seq result was reliable (Fig. 2)

 

Table 1 Summary of RNA-seq metrics for chicken follicles.

Sample

Total reads

Mapped reads, %

Uniq mapped reads, %

GC content, %

% ≥ Q30

S1

62,822,554

57,509,858(91.54%)

56,245,873(89.53%)

50.83

93.75

S2

65,638,382

59,775,235(91.07%)

58,415,734(89.0%)

51.07

93.55

S3

67,279,424

61,497,217(91.41%)

60,160,054(89.42%)

50.39

93.55

F1

66,848,634

61,330,120(91.74%)

59,870,822(89.56%)

50.93

93.71

F2

67,659,778

62,284,292(92.06%)

60,929361(90.05%)

50.65

94.01

F3

80,778,754

75,462,807(93.42%)

73,776,873(91.33%)

50.06

94.24

 

Fig. 1 Gene profile analysis of chicken follicles between SY follicles (S) and F6 follicles (F). a Volcano plot of all the genes detected in the six chicken follicle samples. Green spots represent down-regulation, red spots represent up-regulation. b Hierarchical clustering analysis for DEGs between F and S. c GO enrichment of DEGs from F and S transcriptomes. d KEGG signal pathway enrichment analysis of DEGs.

 

Table 2 The top 10 up- and down-regulated genes of chicken F6 vs SY follicles.

Gene name

Gene ID

log2FoldChange

P-value

padj

Regulated

KRT75L2

431299

5.6448359195

0.00089581216788

0.010968512776

Up

SBK2

420929

5.4012005063

0.0028916599926

0.026749987217

Up

TYRP1

395913

4.0884568588

3.25E-06

0.00011942416133

Up

INHA

424197

3.5955900397

0.0004632196608

0.0065938632923

Up

ACTC1

423298

3.4047368046

1.87E-05

0.0005136082463

Up

SPTSSB

425008

3.1879775126

0.00086973054149

0.010723329659

Up

DMRT2

100858556

3.0375474936

4.20E-06

0.00014842983586

Up

EGR4

422950

2.6567614753

0.00082143038769

0.01027699966

Up

MFSD2B

428656

2.5659125526

1.55E-06

6.55E-05

Up

GJD2

395273

2.4535441052

0.00290562839884

0.0268098383411422

Up

SPIRE1L

418362

-7.1514399248

4.01E-07

2.06E-05

Down

SOX3

374019

-6.9225703884

2.13E-06

8.52E-05

Down

MLPH

424019

-6.3375597357

4.03E-06

0.00014300529036

Down

POU4F3

395521

-6.319423575

7.90E-05

0.0016173871827

Down

LHX3

373940

-6.2949552193

0.00012211061251

0.0022934186943

Down

HAUS3L

101751348

-6.1283949694

1.72E-05

0.00048016708266

Down

GCNT3

427492

-6.0294054072

3.27E-16

1.36E-13

Down

GNOT2

396117

-5.9594390066

3.33E-05

0.00082793115753

Down

CDH15

107054331

-5.8252534515

1.13E-06

5.03E-05

Down

EIF4E1B

107054521

-5.7904110653

5.48E-05

0.0012171943551

Down

 

Table 3. Seven DEGs selected for qRT-PCR validation in chicken follicles

Gene ID

Gene name

Gene description

Log2Foldchange

padj

Regulation

396154

VLDLR1

very low density lipoprotein receptor transcript variant X1

-1.93821454837682

2.52E-10

Down

425805

NGFR

nerve growth factor receptor

-1.99046510645846

1.07E-12

Down

417831

WIF1

WNT inhibitory factor 1

-3.07976009002215

2.58E-20

Down

395887

AMH

anti-Mullerian hormone

-3.49764002338442

8.03E-28

Down

428697

BMP15

bone morphogenetic protein 15

-2.82700723287993

8.58E-07

Down

10705276

GDF6

growth differentiation factor 6

2.17362412753088

0.004248664

Up

395683

MMP13

matrix metallopeptidase 13

2.05982002438789

5.10E-12

Up

 

Fig. 2 The mRNA expression levels of genes examined by qRT-PCR. All data are presented as the mean ± SEM. *, P<0.05

 

Proteomics analysis

Proteins from three SY follicles and three smallest heirachical follicles, which are used for above transcriptome analysis, i.e. S1, S2, S3, and F1, F2 and F3 follicles, were used for TMT labeling and HPLC fractionation followed by LC-MS/MS analysis. The first step was to validate MS data. The distribution of mass error was close to zero, and most of the absolute values were less than 5 ppm, meaning that the mass accuracy of the MS data is compliant with the requirements (Fig. 3a). The length of most peptides was between eight and 16 amino acids, in agreement with the general characteristics of tryptic peptides (Fig. 3b). In this study, a total of 5883 proteins were identified in the samples and 5236 proteins were quantified. According to relative levels, the quantified proteins were divided into two categories: a quantitative ratio over 1.5 was considered up-regulation, and a quantitative ratio less than 1/1.5 was considered down-regulation (P < 0.05) (Fig. 3c). In the F6 follicles, the levels of 175 and 84 proteins significantly increased and decreased, respectively, compared with SY follicles. Detailed analysis of the top 10 up-/down-regulated differentially expressed proteins (DEPs) is shown in Table 4. The entire list of DEPs is shown in Additional file 2: Table S2.

The pathways of DEPs were constructed using the KEGG software. Several important pathways are enriched in the F6 follicles compared with the SY follicles (Fig. 3d). The DEPs were mainly enriched in the ribosome, neuroactive ligand-receptor interaction pathway and cytokine-cytokine receptor interaction.

 

Fig. 3 TMT analysis of the differentially expressed proteins (DEPs) data in chicken F6 and SY follicles. a Mass error distribution of all identified peptides. b The length distribution of the majority of the peptides. c Volcano plots of -log10 (P value) versus log2(expression level) in the F6 vs SY follicles. d KEGG signal pathway enrichment analysis of DEPs.

 

Table 4. Top 10 up- and down regulated proteins in chicken F6 vs SY follicles.

Protein accession code

Protein

description

Gene name

F/S ratio

F/S P value

Regulated

P02659

Apovitellenin-1

---

10.52430044

0.018057763

Up

A0A1D5NZ61

Kinesin-like protein KIF20B

---

9.951417004

0.019122892

Up

A0A1D5PYJ4

Transcriptional repressor p66-alpha

GATAD2A

8.772228989

0.019458417

Up

P41366

Vitelline membrane outer layer protein 1

VMO1

7.163371488

0.016976433

Up

F1NV02

Apolipoprotein B

APOB

6.922300706

0.026918903

Up

R4GM71

Phosphatidylcholine-sterol acyltransferase

LCAT

6.531100478

0.026742079

Up

A0A1L1RYU0

Prostate stem cell antigen

PSCA

6.297115385

0.007551582

Up

A0A1D5NVU2

Keratin, type II cytoskeletal 75

KRT75

6.183636364

0.045556059

Up

P25155

Coagulation factor X

F10

6.062727273

0.021845706

Up

A0M8U1

Suppressor of tumorigenicity 7 protein homolog

ST7

5.905829596

0.02867129

Up

A0A1L1RJJ7

Wnt inhibitory factor 1

WIF1

0.327939317

0.012827638

Down

A0A1D5P589

Tudor and KH domain-containing protein

TDRKH

0.331994981

0.001609664

Down

A0A1D5P0E3

Epithelial cell adhesion molecule

EPCAM

0.367965368

0.00222955

Down

A0A1D5PS81

Protein LSM14 homolog B

LSM14B

0.43363064

0.02217556

Down

F1NWH5

Aquaporin-3

AQP3

0.447158789

0.02917623

Down

F1NC54

SH3 domain and tetratricopeptide repeat-containing protein 1

SH3TC1

0.462576038

0.03313208

Down

E1C8L9

Vacuolar protein sorting-associated protein 29

VPS29L

0.475690608

0.003217184

Down

C7ACT2

Unknown

LOC422926

0.475784992

0.000664487

Down

F1N9X0

Folate receptor alpha

FOLR1

0.482795699

0.008705237

Down

F1P337

Sorting nexin

SNX5

0.49547821

0.00675724

Down

 

Nine DEPs were randomly selected for parallel reaction monitoring (PRM) analysis to verify the accuracy of proteome analysis by LC-MS/MS, including apovitellenin-1 (APO1), apolipoprotein B (APOB), prostate stem cell antigen (PSCA), coagulation factor X (F10), vitellogenin-1 (VTG1) and vitellogenin-3 (VTG3) which were significantly increased in F6 follicles, and zona pellucida sperm-binding protein 2 (ZP2), zona pellucida sperm-binding protein 3 (ZP3L2) and very low-density lipoprotein receptor (VLDLR) which were significantly decreased in F6 follicles (Table 5). The PRM results (Fig. 4) showed that the relative abundance of the peptides from above selected nine individual proteins is consistent with the proteome data.

 

Table 5. Nine proteins in chicken follicle proteome data selected for parallel reaction monitoring analysis

Protein accession code

Protein

description

Gene name

Molecular mass (kDa)

F/S ratio

F/S P value

Regulated

A0A1D5P9N5

Vitellogenin-1

VTG1

209.88

3.858272907

0.029564918

Up

A0A1L1RYU0

Prostate stem cell antigen

PSCA

13.282

6.297115385

0.007551582

Up

F1NV02

Apolipoprotein B

APOB

523.35

6.922300706

0.026918903

Up

P02659

Apovitellenin-1

APOV1

11.966

10.52430044

0.018057763

Up

P25155

Coagulation factor X

F10

53.141

6.062727273

0.021845706

Up

Q91025

Vitellogenin-3

VTG3

38.15

4.500719424

0.03965785

Up

E1BUH5

Zona pellucida sperm-binding protein 3

ZP3L2

51.115

0.519704433

0.007108315

Down

F1NNU1

Zona pellucida sperm-binding protein 2

ZP2

77.033

0.534041942

0.040004135

Down

P98165

Very low-density lipoprotein receptor

VLDLR

94.904

0.551820728

0.019373296

Down

 

Fig. 4 The histogram of nine significantly abundant proteins in F6 follicles (F) vs SY follicles (S) using PRM (P < 0.05).

 

Transcriptome and proteome association analysis

The association analysis of the protepme and transcriptome data of F6 and SY follicles revealed a weak relationship between protein and mRNA expression with a Pearson’s correlation coefficient of 0.23 (Fig. 5). To further understand the relationship between transcripts and proteins, we compared the intersection between differentially expressed genes (DEGs) and differentially expressed proteins(DEPs) (Fig. 6). Most genes were significantly expressed at mRNA level but not at protein level. At both the protein and mRNA levels, 14 and 26 genes were revealed as significantly up- and down- regulated in SY follicles as compared with that in F6 follicles, respectively. In addition, the expression of two genes is inconsistent at changes in mRNA levels and protein levels. Table 6 shows specific regulation information of some genes at mRNA and protein levels. The specific comparative analysis results are shown in Additional file 3: Table S3.

 

Fig.5 Association analysis of proteomic and transcriptomic differences between F6 and SY follicles in laying hens.

 

Fig. 6 Venn diagram of differentially expressed genes/proteins from TMT and DEGs analyses.

 

Table 6 Differentially expressed genes at both mRNA and protein levels between F and S follicles in chicken

Gene

Transcript

Protein

Log2FoldChange

FDR

Regulation type

F/S ratio

P-value

Regulation type

 

 

 

 

 

 

 

VLDLR

-1.93821454837682

0.000000000251

Down

 

0.552

0.0194

Down

 

NGFR

-1.99046510645846

1.07E-12

Down

 

0.566

0.0264

Down

 

WIF1

-3.079759991

2.56E-20

Down

 

0.328

0.0128

Down

 

KPNA7

-5.083462853

1.76E-18

Down

 

0.581

0.0362

Down

 

ACTC1

3.404736805

0.000513608

UP

1.673

0.00158

UP

LCAT

2.363270714

2.85E-19

UP

6.531

0.0267

UP

ZPD

2.302937808

0.007012137

UP

2.727

0.0061

UP

ACTA1

1.543125393

0.000000000407

UP

1.9

23

0.00953

UP

AGT

-3.192564345

0.003115348

Down

 

2.478

0.00904

UP

A2ML1

-1.15166374

0.000215557

Down

 

1.929

0.0304

UP

 

Dynamics and regulation of VLDLR mRNA by FSH

Both transcriptomic and proteomic analyses indicated that VLDLR expression was significantly down-regulated in F6 follicles compared with SY follicles, therefore, we further analyzed the expression of VLDLR mRNA in chicken tissues and found that chcken VLDLR is predominantly expressed in the ovary (Fig. 7a), and its expression level in the prehierarchical follicles was signifcantly higher than that in the hierarchical follicles (P < 0.01) (Fig. 7b). In both the hierarchical and prehierarchical follicles, the VLDLR were signifcantly higher in the granulosa cells (GCs) than in the thecal cells (TCs) (P < 0.05) (Fig. 7c). Follicle-stimulating hormone treatment increased the expression of VLDLR in the GCs of both hierarchical and prehierarchical follicles (P < 0.01) (Fig. 7d, 7e).

 

Fig. 7 Expression dynamics of VLDLR mRNA in chicken and effect of follicle-stimulating hormone (FSH) treatment on the VLDLR mRNA levels in the granulosa cells (GCs) of ovarian follicles. a Expression of VLDLR mRNA in different chicken tissues. b Expression levels of VLDLR in 1–2 mm follicles, 6–8 mm follicles (small yellow follicles), the fifth largest follicles (F5), the third largest follicles (F3), the largest follicles (F1), and the newly postovulatory follicles (POFs). c Expression of chicken VLDLR in the granulosa cells (GCs) and thecal cells (TCs) of hierarchical follicles and the GCs (pre-GCs) and theca cells (pre-TCs) of prehierarchical follicles. d Effect of FSH on VLDLR in the GCs of prehierarchical follicles. e Effect of FSH on VLDLR in the GCs of hierarchical follicles. All data are presented as the mean ± SEM. (abcP < 0.05; ABCP < 0.01).

Discussion

Follicular development is important for egg production and laying performance in the poultry industry. In hens, follicular growth is divided into three stages: slow growth stage from small white follicles to small yellow follicles, selection stage from small yellow follicles to 12-15 mm hierarchical follicles (usually F6 follicles) and rapid growth from the F6 follicles to ovulatory follicles. The morphology and function of these follicles change greatly during follicle selection. A comprehensive understanding of the dynamic changes in the gene and protein expression profiles during follicle selection is essential for elucidating variations in egg production. In this study, by a combinational analysis on changes in follicle transcriptome and proteome, we identified some genes and proteins that may play important roles in follicle selection.

Through bioinformatics analysis and qRT-PCR validation on the DEGs identified by RNA-seq, hundreds of DEGs were revealed to be involved in follicular development, some of which have been shown to be involved in follicular development and selection (Fig. 2). Among these DEGs, anti-mullerian hormone (AMH) is a member of the TGF-β superfamily and is mainly expressed in granulosa cells of 1-5 mm follicles in the early stage of follicular development [1]. MMP13, also known as collagenase-3, initiates the breakdown of thefibrillar collagens that form a key structural element of membranes. Our previous study also revealed that MMP13 is mainly expressed in chicken F5 and POF1 follicle and polymorphisms in the promoter region are associated with age at the first laying trait in laying hens [25]. Bone morphogenetic proteins 15 (BMP15) may promote follicle selection and effects on granulosa cells (GCs) proliferation and steroidogenesis in hens [13,26]. In this study, we also found that, compared with SY follicles, the expression of AMH and BMP15 were significantly decreased in F6 follicles, while that of MMP13 was significantly increased (Table 3). We also detected significant changes in the expression of tumor necrosis factor receptor superfamily member 16 (NGFR), Wnt inhibitory factor 1 (WIF1) and growth differentiation factor 6 (GDF6), the function of which in chicken follicle selection requires further investigation.

By using the LC-MS/MS method, we compared the global proteomes of SY follicles and F6 follicles and identified 259 DEPs that may play certain roles in chicken follicle selection. The levels of apovitellenin-1 (APOV1), apolipoprotein B (APOB), prostate stem cell antigen (PSCA), coagulation factor X (F10), vitellogenin-1 (VTG1) and vitellogenin-3 (VTG3) were significantly increased, and those of zona pellucida sperm-binding protein 2 (ZP2), zona pellucida sperm-binding protein 3 (ZP3L2) and very low-density lipoprotein receptor (VLDLR) were significantly decreased in F6 follicles compared with SY follicles. Apovitellenin-1 is one of five apoproteins that forms low-density lipoproteins (LDL) and is expressed in the egg yolk and vitelline membrane [27, 28]. Apolipoprotein B is the major apolipoprotein of very low-density lipoproteins (VLDL) and LDL and is a ligand for LDL receptors. The increased expression of apovitellenin-1 and apolipoprotein B in 12-15 mm chicken hierarchical follicles is consistent with higher yolk incorporation after follicle selection. Prostate stem cell antigen is reported to be involved in some functions of stem cells, such as self-renewal and proliferation [29]. Coagulation factor X is a vitamin K-dependent serine protease that plays a key role in the common coagulation pathway [30, 31]. Vitellogenins are the yolk precursor proteins produced by the liver and are very important for growth of the chicken ovarian follicles. Vitellogenins circulate in the bloodstream until a follicle enters a stage of vitellogenesis that triggers endocytosis of vitellogenins and transport into the yolk [7]. Vitellogenins 1, 2, and 3 were found in the chicken egg yolk plasma and granules [28]. The increased expression of vitellogenins 1 and 3 in F6 follicles from this study is predictable because yolk deposition is active in F6 follicles. The zona pellucida (ZP) is a specialized extracellular matrix that surrounds the oocyte and early embryo. It is composed of three or four glycoproteins including zona pellucida sperm-binding proteins 1-4 (ZP1-4) with various functions in oogenesis and fertilization [32]. zona pellucida sperm-binding protein 3 can promote oocyte maturation in pigs [33]. In this study, we found that the expression levels of zona pellucida sperm-binding proteins 2 and 3 are grdually decreased.

Very low-density lipoprotein receptor (VLDLR) plays a pivotal role in chicken reproduction, and without VLDLR, oocytes are unable to enter the rapid growth stage of follicle development [34]. During the development of the small white follicle, the VLDLR migrates to the follicular wall, enabling endocytosis of vitellogenin into the yolk, followed by follicular differentiation [7]. In chicken granulosa cells, the expressed variant of VLDLR differs from the variant expressed in the oocyte which contains an O-linked sugar domain (VLDLR 1) [36, 37]. The lower level of VLDLR in granulosa cells is suggested to allow more VLDLR to pass through the intercellular gaps to reach the oocytes rather than receptor mediated endocytosis into the granulosa cells [35]. In this study, by RNA-seq and qRT-PCR, we found that VLDLR1 mRNA is significantly down-regulated, while proteomic analysis and PRM validation both indicated that VLDLR protein mRNA and protein is significantly down-regulated. Further analysis on VLDLR expression indicated that it is mainly expressed in chicken ovary, in SW and SY follicles and in the GCs of prehierachical follicles (Fig. 7). Moreover, its expression in the GCs of both prehierachical and hierachical follicles was stimulated by FSH (Fig. 7d, 7e). These data collectively suggest an important role of VLDLR in chicken follicle selection. Similarly, in geese, the expression of VLDLR mRNA is decreased concomitant to an increase in follicular diameter [38].

In this study, by association analysis between the transcriptome and the proteome, we found that the expression changes of most genes are inconsistent at the mRNA and protein levels, which is consistent with other studies showing that most genes can exhibit differential expression at the transcript and protein level [39-43], suggesting that it is necessary to analyze the function of genes at both levels.

Conclusions

By comparison between the transcriptomes and proteomes of SY follicles and F6 follicles in laying hens, this study revealed 855 DEGs and 259 DEPs. The possible functional significance of some DEGs and DEPs including VLDLR was further discussed. These data may contribute to identification of the functional genes and proteins involved in chicken follicular development and selection.

Methods

Tissue collection

Fifteen randomly sampled Hy-line brown hens were randomly divided into three biological groups, five hens in each group. These hens have been laying regularly for at least one month (28 weeks old, with mean body weight of 2.1± 0.12 kg) that were housed under standard conditions with free access to food and water and vaccination was performed according to recommendations from Hy-line International. From each hen, small yellow follicles 6-8 mm in diameter and hierarchical follicles 12-15 mm in diameter were separately collected, egg yolk was carefully squeezed out with tweezers, washed with phosphate buffered saline (Hyclone), immediately frozen in liquid nitrogen and used for transcriptomic and proteomic analyses. All sampled hens were killed by cervical dislocation immediately after oviposition. The Institutional Animal Care and Use Ethics Committee of Shandong Agricultural University reviewed and approved all procedures described in this study. This study was performed according to the Guidelines for Experimental Animals of the Ministry of Science and Technology of China. Three biological replicates for transcriptomic and proteomic analyses were prepared for total RNA and proteins.

Transcriptome analysis

Total RNA was isolated with TRIzol reagent(Invitrogen, UK) in accordance with the manufacturer’s instructions. Then RNA purity and integrity were checked using the NanoPhotometer® spectrophotometer(IMPLEN, CA, USA) and the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA), respectively. Subsequently, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer(5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase(RNase H-). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3’ends of DNA fragments, NEBNext Adaptor with hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 250~300 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). The library quality was assessed on the Agilent

Bioanalyzer 2100 system(Agilent, USA). the library preparations were sequenced on an Illumina Novaseq platform and 150 bp paired-end reads were generated, and clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. Index of the reference genome was built using Hisat2 v2.0.5 and paired-end clean reads were aligned to the reference genome using Hisat2 v2.0.5. Gene expression level was quantified with featureCounts v1.5.0-p3 and FPKM.

Genes with an adjusted P-value <0.05 and |log2FoldChange| > 1 found by DESeq2 were assigned as differentially expressed. Gene Ontology (GO) enrichment analysis  and Kyoto Encyclopedia of Genes and Genomes database (KEGG) pathway enrichment analysis of differentially expressed genes were implemented by the clusterProfiler R package,and the P-value less than 0.05 were considered significantly enriched by differential expressed genes.

Sample processing and Liquid Chromaography Coupled with Tandem Mass Spectrometry (LC-MS/MS)

Proteins in the tissues were extracted with lysis buffer containing 8 M urea (Sigma) and 1% protease inhibitor cocktail (Calbiochem) and the protein concentration was determined with a BCA kit (Beyotime) according to the manufacturer’s instructions. Protein enzymolysis was performed using trypsin (Promega). After trypsin digestion, peptides were desalted on a Strata X C18 SPE column (Phenomenex) and vacuum-dried. Peptides were reconstituted in 0.5 M TEAB and processed according to the manufacturer’s protocol for the TMT kit (Thermo). The tryptic peptides were fractionated by high pH reverse-phase HPLC using an Agilent 300 Extend C 18 column (5 μm particles size, 4.6 mm ID, 250 mm length). The tryptic peptides were dissolved in 0.1% formic acid (Fluka) (solvent A) and directly loaded onto an in-house packed reversed-phase analytical column (150 mm length, 75 μm i.d.). The gradient included the following steps: from 6% to 23% solvent B (0.1% formic acid in 98% acetonitrile) over 26 min, from 23% to 35% solvent B over 8 min and from 35% to 80% over 3 min; then, the column was eluted with 80% solvent B for the last 3 min; the flow rate was constant at 400 nL/min; an EASY-nLC 1000 UPLC system was used. The peptides were subjected to an NSI source followed by tandem mass spectrometry (MS/MS) in a Q ExactiveTM Plus spectrometer (Thermo) coupled online to the UPLC system. The electrospray voltage was 2.0 kV. The m/z scan range was 350 to 1800 for the full scan, and the intact peptides were detected using an Orbitrap at a resolution of 70,000. Peptides were then selected for MS/MS using NCE set at 28, and the fragments were detected using an Orbitrap at a resolution of 17,500. A data-dependent procedure alternated between an MS scan followed by 20 MS/MS scans with 15.0 s dynamic exclusion. Automatic gain control (AGC) was set at 5E4. Fixed first mass was set as 100 m/z.

The resulting MS/MS data were processed using the Maxquant search engine (v.1.5.2.8). Tandem mass spectra were searched against the Gallus gallus database (http://www.uniprot.org/proteomes/UP000000539) concatenated with the reverse decoy database. Trypsin/P was specified as the cleavage enzyme, and up to 2 missing cleavages were allowed. Mass tolerance for the precursor ions was set to 20 ppm in the first search and to 5 ppm in the main search, and the mass tolerance for the fragment ions was set to 0.02 Da. Carbamidomethyl was specified as a fixed modification of Cys, and oxidation of Met was specified as a variable modification. FDR was < 1%, and the minimum score for the peptides was > 40.

The GO proteome was derived from the UniProt-GOA database (www.http://www.ebi.ac.uk/GOA/). The Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used to identify enriched pathways. Proteins with a threshold of P<0.05 and fold change of > 1.5 or <1/1.5 were identified as differentially expressed proteins (DEPs) between SY and F6 follicles.

Parallel reaction monitoring

Protein extraction and trypsin digestion were performed as described above. Then, the tryptic peptides were dissolved in 0.1% formic acid (solvent A) and directly loaded onto an in-house packed reversed-phase analytical column (150 mm length, 75 μm i.d.). The gradient included the following steps: from 6% to 23% solvent B (0.1% formic acid in 98% acetonitrile) over 38 min, from 23% to 35% solvent B over 14 min, from 35% to 80% solvent B over 4 min, and 80% solvent B for 4 min; constant flow rate of 400 nL/min was maintained using an EASY-nLC 1000 UPLC system. The peptides were subjected to an NSI source followed by tandem mass spectrometry (MS/MS) using a Q ExactiveTM Plus spectrometer (Thermo) connected to the UPLC system. The electrospray voltage was 2.0 kV. The m/z scan range was 350 to 1000 for the full scan, and intact peptides were detected using an Orbitrap at a resolution of 35,000. Then, the peptides were selected for MS/MS using NCE set at 27, and the fragments were detected using an Orbitrap at a resolution of 17,500. A data-independent procedure alternated between an MS scan followed by 20 MS/MS scans. Automatic gain control (AGC) was set at 3E6 for the full MS and at 1E5 for MS/MS. The maximal IT was set at 20 ms for the full MS and switched to automatic mode for MS/MS. The isolation window for MS/MS was set to 2.0 m/z.

The resulting MS/MS data were processed using Skyline (v.3.6). Peptide settings: the enzyme was set as trypsin [KR/P]; the maximal missed cleavage was set as 2; the peptide length was set as 8-25; variable modifications were set as carbamidomethylation of Cys and oxidation of Met; and the maximal variable modification was set as 3. Transition settings: precursor charges were set as 2, 3; ion charges were set as 1, 2; and ion types were set as b, y, p. The product ions were derived from ion 3 to the last ion, and the ion match tolerance was set to 0.02 Da.

Real-time quantitative PCR

The total RNA was extracted from the S and F ovarian follicles that were used for proteome analysis using TRIzol reagent (Invitrogen). Synthesis of the cDNA was performed using 1 μg of the RNA by a PrimeScript RT reagent kit with a gDNA Eraser (TaKaRa, Dalian, China) according to the manufacturer’s protocol. Real-time quantitative PCR was performed using a SYBR Premix Ex TaqTM II kit (TaKaRa, Dalian, China) on a Light Cycler 480 real-time PCR system (Roche) as follows: 95°C for 30 s, 40 cycles of denaturation at 95°C for 10 s, and annealing and extension at 58°C for 20 s. The melting curves were obtained, and quantitative analysis of the data was performed using the (2−ΔΔCT) relative quantification method [44]. The mRNAs of VLDLR, VLDLR1WIF1, NGFR, AMH, BMP15, GDF6 and MMP13 were quantitatively analyzed by this method. Quantification was done by standardizing the reactions versus β-actin. All primers are listed in Additional file 4:  Table S4.

Statistical analysis

Differences in mRNA expression between fillicles and follicular cells were evaluated by one-way ANOVA followed by Duncan’s multiple range test (P < 0.05) using the General Linear Model procedure of SAS (version 9.2). In one experiment, each treatment was repeated four times, and at least three independent experiments were performed. All data were presented as the mean ± SEM (n = 4). The GO and KEGG analysis were performed by Fisher’s t-test. P<0.05 was considered significantly different.

Abbreviations

BP: Biological process; CC: Cellular component; DEGs: Differentially expressed genes; DEPs: Differentially expressed proteins; GO: Gene ontology; KEGG: Kyoto encyclopedia of genens and genomes; GCs: granulosa cells; TCs: thecal cells;

Declarations

Acknowledgments

The authors thank PTM-Biolabs Co., Ltd. (Hangzhou, China) for mass spectrometry analysis and help in preparing the manuscript. The authors would also like to thank Novogene Bioinformatics Technology Co., Ltd. (Beijing, china) for their bioinformatics supports.

Funding

This research was financially supported by grants from the National Natural Science Foundation of China (NSFC 31672414, 31272435, 31772588), the Funds of Shandong“Double Tops” Program (SYL2017YSTD12) and Shandong Agricultural Breed Project (2019LZGC019). The funding body had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Availability of data and materials

For transcriptome sequencing, the raw reads have been submission to the Sequence Read Archive (SRA, https://trace.ncbi.nlm.nih.gov/Traces/sra/) repository, https://trace.ncbi.nlm.nih.gov/Traces/sra_sub/sub.cgi?acc=SRP236909. The MS proteomics data have been deposited to the ProteomeXchange Consortium( http://www.proteomexchange.org/ ) via the PRIDE partner repository with the dataset identifiers PXD011470.

Authors’ contributions

QC, LK and YJ designed and drafted the manuscript. QC, YW and ZL carried out animal care, prepared samples and performed the experiments. QC and XG performed the data processing and biological information analysis. YJ, LK and YS conceived the study and the experimental design and helped draft the manuscript. All authors read and approved the final manuscript.

Ethics approval and consent to participate

The animal experiments were carried out in accordance with the protocols of the ‘Guidelines for Experimental Animals’ of the Ministry of Science and Technology (Beijing, China) and all efforts were made to minimize suffering. The animal experiments were approved by the Institutional Animal Care and Use Ethics Committee of Shandong Agricultural University (Permit Number: 2,007,005). We have obtained the informed consents of a written form from Mr. Xiaoyan Yang who was the owner of the Hy-line brown hens.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References

  1. Johnson AL, Woods DC. Dynamics of avian ovarian follicle development: Cellular mechanisms of granulosa cell differentiation. Gen Comp Endocrinol. 2009;163:12-17.
  2. Onagbesan O, Bruggeman V, Decuypere E. Intra-ovarian growth factors regulating ovarian function in avian species: a review. Anim Reprod Sci. 2009;111:121-140.
  3. Johnson AL. Ovarian follicle selection and granulosa cell differentiation. Poult Sci. 2015;94:781.
  4. Fortune JE, Rivera GM, Evans AC, Turzillo AM. Differentiation of dominant versus subordinate follicles in cattle. Biol Reprod. 2001;65:648-654.
  5. Johnson AL, Woods DC. Ovarian dynamics and follicle development. In: Jamieson, B.G.M. (Ed.), Reproductive Biology and Phylogeny of Aves. 2007;Science Publishers, Inc., pp: 243–277 (Chapter 6).
  6. Stifani S, Barber DL, Nimpf J, Schneider W. A single chicken oocyte plasma membrane protein mediates uptake of very low density lipoprotein and vitellogenin. Proc Natl Acad Sci. 1990;87:1955-1959.
  7. Schneider W. Receptor-mediated mechanisms in ovarian follicle and oocyte development. Gen Comp Endocrinol. 2009;163:18–23.
  8. Diaz FJ, Anthony K, Halfhill AN. Early avian follicular development is characterized by changes in transcripts involved in steroidogenesis, paracrine signaling and transcription. Mol Reprod Dev. 2011;78:212–223.
  9. Zhu G, Yong M, Zhou W, Jiang Y. Dynamic changes in the follicular transcriptome and promoter DNA methylation pattern of steroidogenic genes in chicken follicles throughout the ovulation cycle. Plos One. 2015;10:e0146028.
  10. Wang Y, Chen Q, Liu Z, Guo X, Du Y, Yuan Z, et al. Transcriptome analysis on single small yellow follicles reveals that WNT4 is involved in chicken follicle selection. Front Endocrinol. 2017;8:317.
  11. Ocón-Grove OM, Poole DH, Johnson AL. Bone morphogenetic protein 6 promotes FSH receptor and anti-müllerian hormone mRNA expression in granulosa cells from hen prehierarchal follicles. Reproduction. 2012;14:3825.
  12. Kim D, Ocón-Grove O, Johnson AL. Bone morphogenetic protein 4 supports the initial differentiation of hen (Gallus gallus) granulosa cells. Biol Reprod. 2013;88: 161.
  13. Stephens CS, Johnson PA. Bone morphogenetic protein 15 may promote follicle selection in the hen. Gen Comp Endocrinol. 2016;235:170–176.
  14. Knight PG, Glister C. TGF-beta superfamily members and ovarian follicle development. Reproduction. 2016;132:191–206.
  15. Zhu G, Chen X, Mao Y, Kang L, Ma X, Jiang Y. Characterization of annexin A2 in chicken follicle development: Evidence for its involvement in angiogenesis. Anim Reprod Sci. 2015;161:104-11.
  16. Ewen K, Baker M, Wilhelm D, Aitken RJ, Koopman P. Global survey of protein expression during gonadal sex determination in mice. Mol Cell Proteomics. 2009;8:2624.
  17. Ma M, Guo X, Wang F, Zhao C, Liu Z, Shi Z, et al. Protein expression profile of the mouse metaphase-ii oocyte. J Proteome Res. 2008;7:4821-4830.
  18. Wang L, Zhu YF, Guo XJ, et al. A two-dimensional electrophoresis reference map of human ovary. J Mol Med. 2005;83:812-821.
  19. Pfeiffer MJ, Siatkowski M, Paudel Y, Balbach ST, Baeumer N, Crosetto N, et al. Proteomic analysis of mouse oocytes reveals 28 candidate factors of the "reprogrammome". J Proteome Res. 2011;0:2140.
  20. Massicotte L, Coenen K, Mourot M, Sirard MA. Maternal housekeeping proteins translated during bovine oocyte maturation and early embryo development. Proteomics. 2006;6:3811-3820.
  21. Memili E, Peddinti D, Shack LA, Nanduri B, Mccarthy F, Sagirkaya H, et al. Bovine germinal vesicle oocyte and cumulus cell proteomics. Reproduction. 2007;133: 1107-1120.
  22. Peddinti D, Memili E, Burgess SC. Proteomics-based systems biology modeling of bovine germinal vesicle stage oocyte and cumulus cell interaction. Plos One. 2010;5:e11240.
  23. Upadhyay R D, Balasinor NH, Kumar AV, Sachdeva G, Parte P, Dumasia K. Proteomics in reproductive biology: beacon for unraveling the molecular complexities. Biochim Biophys Acta. 2013;1834:8-15.
  24. Nepomuceno AI, Muddiman DC, Petitte JN. Global proteomic analysis of functional compartments in immature avian follicles using laser microdissection coupled to LC-MS/MS. J Proteome Res. 2015;14:3912-3923.
  25. Yuan Z, Chen Y, Chen Q, Guo M, Kang L, Zhu G, et al. Characterization of chicken MMP13 expression and genetic effect on egg production traits of its promoter polymorphisms. G3&#58; Genes|Genomes|Genetics. 2016;6:1305-1312.
  26. Elis S, Dupont J, Couty I, Persani L, Govoroun M, Blesbois E, et al. Expression and biological effects of bone morphogenetic protein-15 in the hen ovary. J Endocrinol. 2007;194:485-497.
  27. Mann K. Proteomic analysis of the chicken egg vitelline membrane. Proteomics. 2008;8:2322–2332.
  28. Mann K, Mann M. The chicken egg yolk plasma and granule proteomes. Proteomics. 2010;8:178-191.
  29. Gu Z, Thomas G, Yamashiro J, Shintaku IP, Dorey F, Raitano A, Witte ON, et al. Prostate stem cell antigen (PSCA) expression increases with high gleason score, advanced stage and bone metastasis in prostate cancer. Oncogene. 2000;19:1288-1296.
  30. Larson P, Camire RM, Wong D, Fasano NC, Monroe DM, Tracy PB, et al. Structure/function analyses of recombinant variants of human factor xa: factor xa incorporation into prothrombinase on the thrombin-activated platelet surface is not mimicked by synthetic phospholipid vesicles. Biochemistry. 1998;37:5029-5038.
  31. Rudolph AE, Mullane MP, Porche-Sorbet R, Daust HA, Miletich JP. The role of the factor x activation peptide: a deletion mutagenesis approach. Thromb Haemostasis. 2002;88:756-762.
  32. Wassarman PM, Jovine L, Qi H, Williams Z, Darie C, Litscher ES. Recent aspects of mammalian fertilization research. Mol Cell Endocrinol. 2005;234:95-103.
  33. Sakaguchi Y, Uzuhashi R, Iwata H, Monji Y, Kuwayama T. Changes in the sperm-zona pellucida binding properties during porcine oocyte maturation. Journal of Mammalian Ova Research. 2010;27:130-135.
  34. Nimpf J, Radosavljevic MJ, Schneider WJ. Oocytes from the mutant restricted ovulator hen lack receptor for very low density lipoprotein. J Biol Chem. 1989;264:1393-1398.
  35. Lin X, Ma Y, Qian T, Yao J, Mi Y, Zhang C.Basic fibroblast growth factor promotes prehierarchical follicle growth and yolk deposition in the chicken. Theriogenology. 2019;139:90-97.
  36. Hayashi K, Nimpf J, Schneider WJ. Chicken oocytes and fibroblasts express different apolipoprotein-B-specific receptors. J Biol Chem. 1989;264:3131-3139.
  37. Eresheim C, Leeb C, Buchegger P, Nimpf J. Signaling by the extracellular matrix protein reelin promotes granulosa cell proliferation in the chicken follicle. J Biol Chem. 2014;289:10182-10191.
  38. Hu S, Liu H, Pan Z, Xia L, Dong X, Li L, et al. Molecular cloning, expression profile and transcriptional modulation of two splice variants of very low density lipoprotein receptor during ovarian follicle development in geese (anser cygnoide). Anim Reprod Sci. 2014;149:281-296.
  39. Song J, Wang C. Transcriptomic and proteomic analyses of genetic factors influencing adductor muscle coloration in QN Orange scallops. BMC Genomics. 2019;20:363.
  40. Abdulghani M, Song G, Kaur H, Walley JW, Tuteja G. (2019). Comparative analysis of the transcriptome and proteome during mouse placental development. J Proteome Res. 2019.
  41. Zhai Y, Zhang Z, Gao H, Chen H, Sun M, Zhang W, et al. (2017). Hormone signaling regulates nymphal diapause in laodelphax striatellus (hemiptera: delphacidae). Sci Rep. 2017;7:13370.
  42. Ghazalpour A, Bennett B, Petyuk VA, Orozco L, Hagopian R, Mungrue IN, et al. (2011). Comparative analysis of proteome and transcriptome variation in mouse. PLoS Genet, 2011;7:e1001393.
  43. Schwanhäusser B, Busse D, Li N, Dittmar G, Schuchhardt J, Wolf J, et al. (2013). Global quantification of mammalian gene expression control. Nature, 2013;495: 126-127.
  44. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 2001;25:402-408.