Complementary Transcriptomic and Proteomic Analysis Reveals a Complex Network Regulating Pollen Abortion in GMS (msc-1) Pepper (Capsicum annuum L.)

Pepper (Capsicum annuum L.) is a globally important horticultural crop. Use of the genic male-sterile (GMS) line enables efficient commercial hybrid pepper seed production. However, the mechanisms of pepper GMS functioning remain unclear. In this study, we used proteomic and transcriptomic analysis to identify proteins and genes related to genic male sterility. A total of 764 differentially expressed proteins (DEPs) and 1069 differentially expressed genes (DEGs) were identified in the proteomic and transcriptomic level respectively, and 52 genes (hereafter “cor-DEGs-DEPs” genes) were detected at both levels. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis identified 13 DEPs and 14 DEGs involved in tapetum and pollen development. Among the 13 DEPs identified, eight were involved in pollen exine formation, and they were all up-regulated in the fertile line 16C1369B. For the 14 DEGs identified, ABORTED MICROSPORES (AMS) and DEFECTIVE IN TAPETAL DEVELOPMENT AND FUNCTION1 (TDF1) were involved in tapetum development, and both are possibly regulated by Msc-1. All of these genes were detected and confirmed by qRT-PCR. The presence of these genes suggests their possible role in tapetum and pollen exine formation in GMS pepper. Most key genes and transcription factors involved in these processes were down-regulated in the sterile line 16C1369A. This study provides a better understanding of GMS (msc-1) molecular functioning in pepper.


Introduction
The use of male sterile lines to produce hybrid seeds is one of the most important methods in hybrid breeding; it is used, for instance, in crops such as rice, maize, Chinese cabbage, sunflower and pepper [1,2] The genic male sterile (msc-1) (GMS) pepper is a spontaneous male-sterile mutation found in Shenjiao (C. annuum) in China [3]. To date, this GMS line has been used successfully to create several varieties of Shenjiao pepper; these are No. 3, 4, and 5 etc. widely grown in Northern China.
The defects of stamens provide the biological basis for male sterility in plants; these include defects in both anther and pollen development. The anther consists of four distinct cell layers: epidermis, endothecium, middle layer, and tapetum. Microspores are produced in the locule surrounded by Liu et al. [40] analyzed the anther transcriptomes in pepper CMS lines and its fertility restoration; they found that anther transcriptome analysis may be useful for identifying potential candidate genes associated with the formation or abortion of pollen. It has been reported that the pollen abortion is associated with upregulation of the genes that code for methyltransferase during meiosis in CMS lines [41]. RNA sequencing has been used to identify many of the genes potentially involved in pollen development in fertile pepper [42].
Although we have previously identified a strong candidate gene for msc-1 [29], we focus here on proteins, because they are the final products of genes and are more directly related to cellular metabolism and biological processes. Prior to now, the mechanisms by which these genes and proteins regulate male sterility in the GMS (msc-1) pepper have not been explained. In this study, we use proteomic and transcriptomic analysis to identify the DEPs and DEGs in GMS (msc-1) pepper and its male fertile line at the key stage of pollen abortion. We discuss the relationship between these DEPs and DEGs and male sterility in GMS pepper, examine the possible biological functions of these DEPs and DEGs, and explore their potential effects on anther development and male sterility.

Cytological Observation of Pepper Anthers at Different Developmental Stages
To determine the cause of anther sterility in 16C1369A, pepper anthers of the sterile line 16C1369A and the fertile line 16C1369B were selected at four different pollen developmental stages and paraffin section analysis was performed. At the pollen mother cell stage, the sterile 16C1369A and fertile line 16C1369B did not differ significantly in their cytological structure ( Figure 1A,D), whereas pollen abortion was observed in the sterile line 16C1369A in the tetrad stage ( Figure 1B,E); in this stage, the sterile line 16C1369A tapetal layer cells were over-vacuolized and premature death occurred, whereas in the fertile line 16C1369B, the tapetal layer cells were degraded sufficiently for microspore development. In the sterile line, the 16C1369A tapetal layer cells subsequently tightly surrounded and squeezed the microspores, inhibiting callosum degradation and microspore release from the tetrad at the uninucleate stage ( Figure 1F). After this, the microspores became dissociated and died during the later stages of pollen development. Finally, the collapsed tapetal cells and crushed microspores ultimately formed a dense belt ( Figure 1H). Microspores were released normally from the tetrad in the fertile line 16C1369B, then became rounded, and developed a three-grooved germination pore; the tapetal layer cells dissolved, or partially remained on the endothecium ( Figure 1C,G). Therefore, the key stage of pollen abortion occurred at the tetrad stage, due to premature tapetum death in the GMS pepper.

Overview of Quantitative Proteomics Analysis
In total, 502,216 spectra were generated using label-free analysis in the sterile line 16C1369A and fertile line 16C1369B. A total of 150 852 spectra matched to known spectra, and 29 056 peptides, 21 108 unique peptides, and 4909 protein species were identified (FDR ≤ 0.01) ( Table S1). The correlation coefficients (R 2 ) between the biological replicates were 0.9845 (16C1369A-1 and 16C1369A-2), 0.9795 (16C1369A-1 and 16C1369A-3), 0.9734 (16C1369A-2 and 16C1369A-3), 0.9873 (16C1369B-1 and 16C1369B-2), 0.9895 (16C1369B-1 and 16C1369B-3), and 0.9885 (16C1369B-2 and 16C1369B-3), respectively. The distribution of molecular weights of the identified proteins, number of peptides, length of peptides and the coverage of the proteins by the identified peptides were provided in Figure S1. Based on the threshold for screening DEPs (|log 2 FC| ≥ 1 and p ≤ 0.05), 674 DEPs were identified and quantified; 111 were up-regulated and 71 were down-regulated in the fertile line 16C1369B. Additionally, 276 and 216 protein species accumulated specifically in the sterile line 16C1369A and the fertile line 16C1369B, respectively (Figure 2A; Table S2). The results indicate that these DEPs may form a complex regulatory network for pollen abortion in pepper.

Comparison of Transcriptome and Proteome Data
To evaluate the relationship between transcript levels and protein abundance, we compared all the transcripts and their corresponding quantified proteins in the sterile line 16C1369A and the fertile line 16C1369B at the tetrad stage. Among all the transcripts and the quantified proteins, there were 2478 genes that were found in both the transcriptome and proteome ( Figure 2C). Of the 764 DEPs and 1069 DEGs identified, 52 cor-DEGs-DEPs were regulated at both the transcriptional and translational levels ( Figure 2D).
We conducted a correlation analysis between the transcriptomic and quantitative proteomic data. The expression levels of all the transcripts and their corresponding quantified proteins in sterile line 16C1369A and fertile line 16C1369B showed limited correlation (r = 0.415 and 0.448; Figure 3A,B). However, less correlation was observed between all the transcripts and their corresponding quantified proteins (r = 0.242, Figure 3C). For the 764 DEPs and 1069 DEGs identified, the correlation of the 52 cor-DEGs-DEPs genes was 0.32 ( Figure 3D), and a higher correlation was observed between proteins and their corresponding mRNAs with the same expression trend (r = 0.436; Figure 3E) or opposite expression trend (r = 0.773; Figure 3F).

Comparison of Transcriptome and Proteome Data
To evaluate the relationship between transcript levels and protein abundance, we compared all the transcripts and their corresponding quantified proteins in the sterile line 16C1369A and the fertile line 16C1369B at the tetrad stage. Among all the transcripts and the quantified proteins, there were 2478 genes that were found in both the transcriptome and proteome ( Figure 2C). Of the 764 DEPs and 1069 DEGs identified, 52 cor-DEGs-DEPs were regulated at both the transcriptional and translational levels ( Figure 2D).
We conducted a correlation analysis between the transcriptomic and quantitative proteomic data. The expression levels of all the transcripts and their corresponding quantified proteins in sterile line 16C1369A and fertile line 16C1369B showed limited correlation (r = 0.415 and 0.448; Figure 3A,B). However, less correlation was observed between all the transcripts and their corresponding quantified proteins (r = 0.242, Figure 3C). For the 764 DEPs and 1069 DEGs identified, the correlation of the 52 cor-DEGs-DEPs genes was 0.32 ( Figure 3D), and a higher correlation was observed between proteins and their corresponding mRNAs with the same expression trend (r = 0.436; Figure 3E) or opposite expression trend (r = 0.773; Figure 3F).

Cluster Analysis of Expression Patterns in the cor-DEGs-DEPs Genes
Among the 52 cor-DEGs-DEPs genes, 40 genes showed the same expression trend, whereas 12 genes showed the opposite expression trend at the proteomic and transcriptomic levels ( Figure 4; Table S8). The correlated proteins and genes mostly showed up-regulated expression trends. We suggest that some of these genes might play important roles in pollen development.

Cluster Analysis of Expression Patterns in the cor-DEGs-DEPs Genes
Among the 52 cor-DEGs-DEPs genes, 40 genes showed the same expression trend, whereas 12 genes showed the opposite expression trend at the proteomic and transcriptomic levels ( Figure 4; Table S8). The correlated proteins and genes mostly showed up-regulated expression trends. We suggest that some of these genes might play important roles in pollen development.

GO and KEGG Analysis of the cor-DEGs-DEPs Genes
Gene ontology (GO) analysis indicated that 45 of the 52 cor-DEGs-DEPs genes had GO annotations. The results covered a wide range of cellular components, molecular functions, and biological processes, including 30 important functional groups ( Figure 5A). The largest subcategory found in the "biological process" was "external encapsulating structure organization", followed by "pollen exine formation", "pollen wall assembly", "pollen development", and others. For the "cellular component" category the three largest subcategories were "extracellular region", "external encapsulating structure", and "cell wall". For the "molecular function" category, "carbon-nitrogen lyase activity", "amine-lyase activity", "strictosidine synthase activity", and "fatty acid synthase activity" were the subcategories with the most annotations.
The KEGG pathways of the cor-DEGs-DEPs genes were identified using a p-value ≤ 0.05 as the cut-off, and the 52 cor-DEGs-DEPs genes were mapped to 47 KEGG pathways (Table S9). The results showed that three KEGG pathways were significantly enriched at both mRNA and protein levels; these included the following pathways: "metabolic pathways" (map01100), "Biosynthesis of secondary metabolites" (map01110), and "Phenylpropanoid biosynthesis" (map00940) ( Figure 5B). It has been reported that the "Phenylpropanoid biosynthesis" pathways are associated with pollen exine formation [43][44][45]. Therefore, based on the analysis of GO and the KEGG pathways enrichment, the difference between the sterile line 16C1369A and fertile line 16C1369B was mainly related to "pollen exine formation", "pollen wall assembly", "pollen development", and "Phenylpropanoid biosynthesis". biological processes, including 30 important functional groups ( Figure 5A). The largest subcategory found in the "biological process" was "external encapsulating structure organization", followed by "pollen exine formation", "pollen wall assembly", "pollen development", and others. For the "cellular component" category the three largest subcategories were "extracellular region", "external encapsulating structure", and "cell wall". For the "molecular function" category, "carbon-nitrogen lyase activity", "amine-lyase activity", "strictosidine synthase activity", and "fatty acid synthase activity" were the subcategories with the most annotations.
The KEGG pathways of the cor-DEGs-DEPs genes were identified using a p-value ≤ 0.05 as the cut-off, and the 52 cor-DEGs-DEPs genes were mapped to 47 KEGG pathways (Table S9). The results showed that three KEGG pathways were significantly enriched at both mRNA and protein levels; these included the following pathways: "metabolic pathways" (map01100), "Biosynthesis of secondary metabolites" (map01110), and "Phenylpropanoid biosynthesis" (map00940) ( Figure 5B). It has been reported that the "Phenylpropanoid biosynthesis" pathways are associated with pollen exine formation [43][44][45]. Therefore, based on the analysis of GO and the KEGG pathways enrichment, the difference between the sterile line 16C1369A and fertile line 16C1369B was mainly related to "pollen exine formation", "pollen wall assembly", "pollen development", and "Phenylpropanoid biosynthesis".

DEPs and DEGs Related to Anther and Pollen Development
Pollen development is a complex process that involves many events, especially tapetum development and pollen exine formation. In this study, based on the GO and KEGG analysis, we identified 13 DEPs and 14 DEGs involved the pollen development; most of these, including five cor-DEGs-DEPs genes (TA29, TA1, FAS1, LAP3, and PG2), were up-regulated in the fertile line 16C1369B. Among these genes, we propose that AMS and TDF1 are involved in tapetum development. The genes ACOS5, CYP704B1, PKSA, and PKSB, TKPR1, TKPR2, MS2, and ABCG26 were involved in the pollen exine formation (Table 1).

DEPs and DEGs Related to Anther and Pollen Development
Pollen development is a complex process that involves many events, especially tapetum development and pollen exine formation. In this study, based on the GO and KEGG analysis, we identified 13 DEPs and 14 DEGs involved the pollen development; most of these, including five cor-DEGs-DEPs genes (TA29, TA1, FAS1, LAP3, and PG2), were up-regulated in the fertile line 16C1369B. Among these genes, we propose that AMS and TDF1 are involved in tapetum development. The genes ACOS5, CYP704B1, PKSA, and PKSB, TKPR1, TKPR2, MS2, and ABCG26 were involved in the pollen exine formation (Table 1).

Validation of Gene Expression Levels
qRT-PCR was performed to validate the RNA-seq results and the correspondence between protein levels and mRNA abundance. We selected 14 up-or down-regulated DEGs and 13 DEPs (including five cor-DEGs-DEPs genes) involved in tapetum and pollen development for qRT-PCR analysis. Twelve DEGs were highly expressed and two DEGs (ARO1, AGL62) were down-regulated in the fertile line 16C1369B. This result is consistent with the RNA-Seq results ( Figure 6A). The genes involved in coding for the 13 DEPs were all expressed significantly more in the fertile line 16C1369B than in the sterile line 16C1369A, with the same trend in the label-free data ( Figure 6B). These results indicate that the transcriptome and proteome data were reliable.

Validation of Gene Expression Levels
qRT-PCR was performed to validate the RNA-seq results and the correspondence between protein levels and mRNA abundance. We selected 14 up-or down-regulated DEGs and 13 DEPs (including five cor-DEGs-DEPs genes) involved in tapetum and pollen development for qRT-PCR analysis. Twelve DEGs were highly expressed and two DEGs (ARO1, AGL62) were down-regulated in the fertile line 16C1369B. This result is consistent with the RNA-Seq results ( Figure 6A). The genes involved in coding for the 13 DEPs were all expressed significantly more in the fertile line 16C1369B than in the sterile line 16C1369A, with the same trend in the label-free data ( Figure 6B). These results indicate that the transcriptome and proteome data were reliable.

Discussion
Although genic male sterility provides great potential in terms of promoting pepper heterosis, no prior studies have characterized the msc-1 changes in protein and gene expression in pepper anthers at the proteomic and the transcriptomic levels. As such, in order to gain a better understanding of the mechanisms of GMS regulation in pepper, we have conducted proteomic and transcriptomic analysis comparing sterile and fertile lines, in the tetrad stages of anther development.
In this study, the sterile line 16C1369A tapetal layer cells were over-vacuolized and premature death occurred ( Figure 1E), and finally the collapsed tapetal cells and crushed microspores formed a dense belt ( Figure 1H). In the sterile line, we observed that the tapetal layer cells were over-vacuolized and premature death occurred; this finding suggests abnormal processing of programmed cell death (PCD) in the sterile line, causing the tapetum to have defects in sporopollenin synthesis and nutrient supply, which affect later anther development. These results are consistent with the findings of Guo et al. [39] and Fang et al. [38] and suggest that the abnormal development of the tapetal cells and microspores leads to pollen and anther abortion. In addition, we have previously identified a strong candidate gene for msc-1 [29]; Msc-1 is a homolog of AtDYT1, which is a bHLH transcription factor involved in early tapetal development [20].
The anther structure consists of gametophytes surrounded by four distinct cell layers: outer epidermis, endothecium, middle cell layer, and the innermost tapetum [46]. These layers, especially the tapetum is known to act as a supplier of metabolites, nutrients, and sporopollenin precursors [44]. Tapetum degeneration is induced through the process of PCD during the late developmental stage of the anther. Furthermore, the normal process of PCD in the tapetum is essential for normal exine pattern formation [47,48].
Exine consists primarily of sporopollenin, a tough and chemically inert biopolymer, which ensures the production of functional pollens grains [8]. Establishment of exine ornamentation by sporopollenin polymerization needs three key steps: synthesis, secretion, and translocation of sporopollenin precursors. Synthesis of a majority of the sporopollenin precursors is conducted in the

Discussion
Although genic male sterility provides great potential in terms of promoting pepper heterosis, no prior studies have characterized the msc-1 changes in protein and gene expression in pepper anthers at the proteomic and the transcriptomic levels. As such, in order to gain a better understanding of the mechanisms of GMS regulation in pepper, we have conducted proteomic and transcriptomic analysis comparing sterile and fertile lines, in the tetrad stages of anther development.
In this study, the sterile line 16C1369A tapetal layer cells were over-vacuolized and premature death occurred ( Figure 1E), and finally the collapsed tapetal cells and crushed microspores formed a dense belt ( Figure 1H). In the sterile line, we observed that the tapetal layer cells were over-vacuolized and premature death occurred; this finding suggests abnormal processing of programmed cell death (PCD) in the sterile line, causing the tapetum to have defects in sporopollenin synthesis and nutrient supply, which affect later anther development. These results are consistent with the findings of Guo et al. [39] and Fang et al. [38] and suggest that the abnormal development of the tapetal cells and microspores leads to pollen and anther abortion. In addition, we have previously identified a strong candidate gene for msc-1 [29]; Msc-1 is a homolog of AtDYT1, which is a bHLH transcription factor involved in early tapetal development [20].
The anther structure consists of gametophytes surrounded by four distinct cell layers: outer epidermis, endothecium, middle cell layer, and the innermost tapetum [46]. These layers, especially the tapetum is known to act as a supplier of metabolites, nutrients, and sporopollenin precursors [44]. Tapetum degeneration is induced through the process of PCD during the late developmental stage of the anther. Furthermore, the normal process of PCD in the tapetum is essential for normal exine pattern formation [47,48].
Exine consists primarily of sporopollenin, a tough and chemically inert biopolymer, which ensures the production of functional pollens grains [8]. Establishment of exine ornamentation by sporopollenin polymerization needs three key steps: synthesis, secretion, and translocation of sporopollenin precursors. Synthesis of a majority of the sporopollenin precursors is conducted in the tapetum [9]. Although the mechanisms of sporopollenin biosynthesis are not well understood, it appears that genes in pathways which produce fatty acids, phenylpropanoids, as well as some anther-specific genes, are important for exine synthesis [49][50][51].
In our proteomic data, many key proteins have been identified for pollen exine formation. These include ACOS5 (Capana02g003302), CYP704B1 (Capana01g002203), PKSA (Capana01g003460) and PKSB (Capana08g002676), TKPR1 (Capana05g000665), TKPR2 (Capana01g002831), MS2 (Capana03g003125) and ABCG26 (Capana07g002406), and these DEPs were all up-regulated in 16C1369B (Table 1) which was confirmed by the qRT-PCR analysis ( Figure 6). In Arabidopsis, the proteins that we identified as being involved in pollen exine formation are reported to participate in the synthesis of sporopollenin precursors [9]. Among these, ACOS5 catalyzes mid-/long-chain fatty acids into fatty acyl-CoA in plastids; this is then hydroxylated by members of the cytochrome P450 superfamily, especially CYP703A2 and CYP704B1. The hydroxylated products are then catalyzed by PKSA and PKSB into triketide and tetraketide α-pyrones, which are the substrates of TKPR1and TKPR2, respectively [52][53][54]. MS2, a fatty acyl reductase, encodes a fatty acid reductase and catalyzes the palmitoyl acyl carrier protein into a fatty alcohol [55]. The resultant sporopollenin precursors are transported to the anther locule by a member of ATP-binding cassette transporter superfamily (ABCG26) and by lipid transfer proteins [56][57][58]. So due to the down-regulated of these proteins in the male sterile line 16C1369A, it may inhibit the synthesis and translocation of the sporopollenin precursors.
In our transcriptomic data, two DEGs, AMS (Capana08g000254) and TDF1 (Capana04g001901), were identified and were up-regulated in 16C1369B (Table 1; Figure 6). In Arabidopsis, both of these DEGs are essential for tapetum development [17,21,59]. In our previous study, we used genome resequencing to identify a strong candidate gene for msc-1 [29]; there was a 7-bp deletion in the 3rd exon, which leads to a premature stop codon, and may cause a loss of function mutation in 16C1369A. In addition to the fact thatMsc-1 is a homology of DYT1, it was up-regulated in 16C1369B. In Arabidopsis, DYT1 is the most upstream regulator in tapetum development, and directly binds the promoter region of TDF1; it thereby affects many downstream genes related to tapetum development and pollen wall formation [60]. TDF1 directly regulates AMS via an AACCT cis-element [61,62]. AMS can bind to the promoters of CYP703A2, PKSB, TKPR1, and CYP704B1 in vivo [59]. MS188 interacts with AMS to regulate the expression of CYP703A2 [63]. MS188 directly regulates the expression of PKSA, PKSB, MS2, and CYP703A2; MS188 and other MYB family members play redundant roles in activating the expression of CYP704B1, ACOS5, and TKPR1 [54]. Therefore, in the sterile line 16C1369A, AMS (Capana08g000254) and TDF1 (Capana04g001901) may be regulated by Msc-1, thus affecting tapetum development.
By comparing the transcriptome and proteome, we identified several cor-DEGs-DEPs involved in tapetum development and pollen exine formation. And all of these were significantly up-regulated in 16C1369B.
In this study, Msc-1 was not listed as a DEP or a DEG in the proteome and the transcriptome, however Msc-1 was detected both in the sterile line 16C1369A and fertile line 16C1369B with the log 2 (16C1369B/16C1369A) was 0.558 in the transcriptome. And in our qRT-PCR analysis, Msc-1 was highly expressed in the pollen mother cell stage in the fertile line 16C1369B and was decreased significantly with the anther development, besides that the expression level of Msc-1 was much higher in the fertile line 16C1369B than in the sterile line 16C1369A in the four anther development stage ( Figure 7A). So, based on our interpretations of our findings, we propose a model to explain the mechanisms of GMS in pepper ( Figure 7B). Due to the down-regulation of Msc-1 (DYT1) in the 16C1369A [29], down-regulation of TDF1 and AMS may occur. This may result in abnormal PCD in the tapetum, severe inhibition of sporopollenin biosynthesis and transport, and defects in the formation of exine. Ultimately, this leads to pollen abortion in GMS pepper.

Plant Materials
Genic male-sterile line 16C1369A and its near isogenic male-fertile line 16C1369B were grown in plastic-covered tunnels at the China Agricultural University, Beijing, China, in 2016. Flowers at different developmental stages were collected for paraffin section analysis, using Safranin O-Fast Green Staining. For proteomic analysis, RNA sequencing, and gene expression analysis, anthers were collected from flower buds at the key stage of pollen abortion, frozen in liquid nitrogen immediately, and then stored at −80 °C for later use. For proteomicand RNA-seq analysis, three biological replicates of sterile and fertile anthers (A1-A3; B1-B3) were collected from 150 to 200 flower buds at the tetrad stage, which were selected randomly from 30 pepper plants. For real-time PCR analysis, three biological and three technical replicates were conducted.

Protein Isolation, Enzymolysis, and Label-Free Analysis
Anther protein extraction and label-free quantitative mass spectrometry analysis were performed according to Guo et al. [39]. For protein quantitation, we required at least two unique peptides in three biological replicates to be detected; the relative quantitative protein ratios between the two groups were calculated by comparing the average abundance values obtained in the three biological replicates. For analysis of differentially expressed proteins (DEPs), a 2-fold cutoff was used determine distinguish between up-regulated and down-regulated proteins. Additionally, proteins detected in only one of the groups (the sterile line 16C1369A or fertile line 16C1369B) at least two replicates were considered to be presence/absence proteins. To analyze differences between presence/absence protein species, the fold-change calculated was based on the mean value of the lowest protein intensity detected in each group.

RNA Isolation, Library Preparation, and Sequencing
Total RNA was extracted from anthers using TRIzol (Invitrogen, Waltham, MA, USA) according to the manufacturer's instructions. Sequencing libraries were constructed using NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) following

Plant Materials
Genic male-sterile line 16C1369A and its near isogenic male-fertile line 16C1369B were grown in plastic-covered tunnels at the China Agricultural University, Beijing, China, in 2016. Flowers at different developmental stages were collected for paraffin section analysis, using Safranin O-Fast Green Staining. For proteomic analysis, RNA sequencing, and gene expression analysis, anthers were collected from flower buds at the key stage of pollen abortion, frozen in liquid nitrogen immediately, and then stored at −80 • C for later use. For proteomicand RNA-seq analysis, three biological replicates of sterile and fertile anthers (A1-A3; B1-B3) were collected from 150 to 200 flower buds at the tetrad stage, which were selected randomly from 30 pepper plants. For real-time PCR analysis, three biological and three technical replicates were conducted.

Protein Isolation, Enzymolysis, and Label-Free Analysis
Anther protein extraction and label-free quantitative mass spectrometry analysis were performed according to Guo et al. [39]. For protein quantitation, we required at least two unique peptides in three biological replicates to be detected; the relative quantitative protein ratios between the two groups were calculated by comparing the average abundance values obtained in the three biological replicates. For analysis of differentially expressed proteins (DEPs), a 2-fold cutoff was used determine distinguish between up-regulated and down-regulated proteins. Additionally, proteins detected in only one of the groups (the sterile line 16C1369A or fertile line 16C1369B) at least two replicates were considered to be presence/absence proteins. To analyze differences between presence/absence protein species, the fold-change calculated was based on the mean value of the lowest protein intensity detected in each group.

RNA Isolation, Library Preparation, and Sequencing
Total RNA was extracted from anthers using TRIzol (Invitrogen, Waltham, MA, USA) according to the manufacturer's instructions. Sequencing libraries were constructed using NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) following manufacturer's recommendations.
The purity and quality of the libraries was assessed using an Agilent 2100 Bioanalyzer and Qubit 2.0. Then, the six libraries were sequenced using the Illumina Hiseq2500/X platform and 125/150 bp paired-end reads were generated. Clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N, and low-quality reads from raw data. At the same time, the Q20, Q30, and GC content of the clean data was calculated. All the downstream analyses were based on the clean data using high quality reads. TopHat2 was used to map clean reads to reference gene and reference genome [64]. Cufflinks were used to discover novel genes [65]. Cuffquant and Cuffnorm (v 2.2.1) were used to calculate Fragments Per Kilobase of transcript per Million mapped reads (FPKMs) of genes in each sample. Gene FPKMs were computed by summing the FPKMs of transcripts in each gene group. We used the DESeq2 R package [66] to analyze the differentially expressed genes (DEGs) in the sterile line 16C1369A and the fertile line 16C1369B.

Bioinformatics Analysis
Gene ontology (GO) enrichment analysis of DEGs and DEPs was implemented using the GOseq R package [67] in which gene length bias was corrected. Pathway analysis of DEGs and DEPs was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/ kegg). A p-value ≤ 0.05 was used as the threshold to determine the significant enrichments of GO and KEGG pathways.
The RNA-Seq raw sequence data are deposited in the Short Read Archive (SRA) of National Centre for Biotechnology (NCBI) and are available under the accession: PRJNA524267. And the mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifier PXD012762. The reviewer account: Username: reviewer22878@ebi.ac.uk; Password: U964wSDm.

Quantitative Real-Time PCR (qRT-PCR) Analysis
Total RNA of anthers from the sterile line 16C1369A and the fertile line 16C1369B was exacted using the SV Total RNA Isolation System Kit (Promega, Madison, WI, USA), following manufacturer's instructions. We performed reverse transcription using 1 µg of total RNA using a PrimeScript™ RT Kit (Takara, Kusatsu, Shiga, Japan). Real-time PCR was carried out with using GoTaq ® qPCR Master Mix (Promega), following manufacturer's instructions, on an ABI 7500 real-time PCR system. The thermocycling conditions were as follows: 95 • C for 2 min, 40 cycles of 95 • C for 15 s, and 60 • C for 1 min. Relative expression values were calculated using the 2 −∆∆Ct method, and the UBI-3(AY486137.1) as the reference gene. The primers used for Real-time PCR were designed using Primer 5.0 (http://www.premierbiosoft.com/primerdesign/) and are listed in Table S10.

Conclusions
In this study, we identified numerous proteins and genes associated with pollen development and genic male sterility, using proteomic and transcriptomic analysis. In total, 764 DEPs and 1069 DEGs were identified at the proteomic and transcriptomic levels, respectively. Fifty-two cor-DEGs-DEPs were detected at both the proteomic and transcriptomic level. Bioinformatics analysis suggested that these DEPs and DEGs are involved in pollen exine formation, pollen wall assembly, pollen development, and phenylpropanoid biosynthesis. Based on our results, we propose a gene regulation network model to explain the mechanisms of genic male sterility in pepper (Figure 7). This study will improve our understanding of the genes associated with pollen development in the GMS pepper and contribute to the improvement of pepper hybrid breeding.