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).