Comparative Transcriptome Analysis of the Anthers from the Cytoplasmic Male-Sterile Pepper Line HZ1A and Its Maintainer Line HZ1B

Cytoplasmic male-sterility (CMS) is important for the utilization of crop heterosis and study of the molecular mechanisms involved in CMS could improve breeding programs. In the present study, anthers of the pepper CMS line HZ1A and its maintainer line HZ1B were collected from stages S1, S2, and S3 for transcriptome sequencing. A total of 47.95 million clean reads were obtained, and the reads were assembled into 31,603 unigenes. We obtained 42 (27 up-regulated and 15 down-regulated), 691 (346 up-regulated and 345 down-regulated), and 709 (281 up-regulated and 428 down-regulated) differentially expressed genes (DEGs) in stages S1, S2, and S3, respectively. Through Gene Ontology (GO) analysis, the DEGs were found to be composed of 46 functional groups. Two GO terms involved in photosynthesis, photosynthesis (GO:0015986) and photosystem I (GO:0009522), may be related to CMS. Through Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, oxidative phosphorylation (ko00190) and phenylpropanoid biosynthesis (ko00940) were significantly enriched in the S1 and S2 stages, respectively. Through the analysis of 104 lipid metabolism-related DEGs, four significantly enriched KEGG pathways may help to regulate male sterility during anther development. The mitochondrial genes orf470 and atp6 were identified as candidate genes of male sterility for the CMS line HZ1A. Overall, the results will provide insights into the molecular mechanisms of pepper CMS.


Introduction
Pepper (Capsicum annuum L.) is an important global vegetable crop, and most of the varieties grown are hybrids [1]. Utilization of the maximum heterosis was responsible for approximately 163.8% [2] or 331.11% [3] of fruit pepper yield. The three-line hybrid breeding system (cytoplasmic male sterility (CMS) line, maintainer line, and restorer line) is critical and widely used for the efficient production of pepper hybrids. In this system, CMS plants must be maintained and restored for seed production. Improved knowledge of the molecular mechanisms governing CMS could help improve breeding programs in the future.
CMS is a common phenomenon that is observed in a diverse array of plant species, such as rice, maize, cotton, and soybean [4][5][6]. It occurs via interactions between the mitochondrial and nuclear genes. There are two types of genes associated with CMS: the restorer of male

Cytological Examination
In the morning, anthers of different-sized buds from the CMS line HZ1A and its maintainer line HZ1B were dissected. For identifying the fertility of HZ1A and HZ1B, the pollen grains of fully bloomed flowers were tested using Alexander's stain [49]. The fertile pollen grains were stained magenta red, and the sterile pollen grains were stained blue-green. For observing the corresponding relationship between the development of flower buds and pollen, microspores were stained with hematoxylin solution and observed under an optical microscope.
For scanning electron microscopy (SEM) observations, the specimens were prepared according to Nie et al. [50]. Briefly, anthers with mature pollen of HZ1A or HZ1B were treated as follows: fixed with 4% glutaraldehyde (0.1 M sodium phosphate buffer, pH 7.2), overnight at 4 • C, rinsed with 0.1 M sodium phosphate buffer (pH 7.2), substituted in 1% OsO 4 for 2 h at 4 • C, rinsed again, dehydrated in a graded ethanol series, sputter coated with gold, and then observed under a SEM (S3000N, Hitachi, Krefeld, Germany).

cDNA Preparation and RNA Sequencing
Total RNA was isolated from the anthers of each specimen (S1_HZ1A, S1_HZ1B, S2_HZ1A, S2_HZ1B, S3_HZ1A, and S3_HZ1B) using an RNAprep Pure Plant Plus Kit (DP441, Tiangen Biotech, Beijing, China) according to the manufacturer's instructions.
The quantity and quality of the mRNA were assayed using Nanodrop 2000 (NanoDrop Technologies, Wilmington, DE, USA) and LabChip GX Touch (Perkin Elmer, Boston, MA, USA), respectively. Then, 1 µg total RNA was used to construct cDNA libraries by a PCR-cDNA kit (SQK-PCS109) following the manufacturer's instructions. Reverse transcriptase was used to enrich full-length cDNAs. Thereafter, PCR adapters were added directly to both ends of the first-strand cDNA following 14 cycles of cDNA PCR with the LongAmp Taq PCR Kit (New England Biolabs, Ipswich, MA, USA) according to the manufacturer's instructions. The PCR products were subjected to Oxford Nanopore Technologies (ONT, Oxford, UK) adaptor ligation using T4 DNA ligase (New England Biolabs, Ipswich, MA, USA). Thereafter, the products were purified using Agencourt AMPure XP beads (Beckman Coulter, Indianapolis, IN, USA). Finally, the cDNA libraries were added to FLO-MIN109 flow cells and run on the PromethION sequencing platform (ONT, Oxford, UK).

Bioinformatic Analysis of the Sequencing Data
Initially, raw reads were filtered with a minimum average read quality score of 6 and a minimum read length of 500 bp using the NanoFilt package, and the ribosomal RNA reads were discarded. The Phred scores (Q) were calculated from the formula: Q = −10 × log 10 P, where P refers to the probability of base recognition error. Then, fulllength, non-chimeric (FLNC) transcripts were determined by searching for primers at both ends of reads. Clusters of FLNC transcripts were obtained after mapping to the reference genome with mimimap2 [51]. Clusters of full-length, non-chimeric transcripts and consensus isoforms were obtained. Consensus sequences were mapped to the reference genome using the minimap2. Mapped reads were further collapsed using the cDNA_Cupcake package with min-coverage = 85% and min-identity = 90%. A 5 difference was not considered when the redundant transcripts collapsed [52].

Differential Gene Expression Analysis
Differential expression analysis was performed using the DESeq2 R package (1.6.3) [56]. The resultant p values were then adjusted using Benjamini and Hochberg's approach [57] to control the false discovery rate. Both |log 2 (HZ1A vs. HZ1B)| and the p value were used to determine the threshold for the DEGs. Genes with |log 2 (HZ1A vs. HZ1B)| ≥ 1.0 and p value < 0.01 were defined as DEGs. The TFs were identified by predicting all the DEGs using iTAK (Plant Transcription factor & Protein Kinase Identifier and Classifier; http://itak.feilab.net/cgi-bin/itak/index.cgi, accessed on 30 September 2021) [58] based on the databases PlnTFDB and PlantTFDB.

Functional Enrichment Analysis
Functional annotation of the DEGs was performed based on the GO [59], and the enrichment analysis of the DEGs was implemented using the GOseq R packages based on Wallenius non-central hyper-geometric distribution [60], which can adjust for gene length bias in DEGs. KEGG pathway annotation was carried out using a BLAST search against the KEGG database (http://www.genome.jp/kegg/, accessed on 22 September 2021), and pathway enrichment analysis of the DEGs was performed using the KEGG Orthology Based Annotation System [61]. GO terms and KEGG pathways with a corrected p < 0.05 were defined as significantly enriched.

Quantitative Real-Time PCR and Candidate Genes Amplification
All primers used in the present study were designed using Primer-Blast software (NCBI, Bethesda, MD, USA) and synthesized by Sangon Biotech (Shanghai, China). The pepper β-actin gene [62] was used as an internal control to normalize the expression data. The primer pairs used are listed in Table S1. Quantitative real-time PCR (qRT-PCR) was conducted according to the instructions of 2 × Hieff UNICON ® Universal Blue qPCR SYBR Green Master Mix (YEASEN, Shanghai, China), and analyzed using the Quant Studio 5 real-time PCR system (Thermo Fisher Scientific, Inc., Waltham, MA, USA). The 2 −∆∆Ct method was used for calculating the relative mRNA expression. Three biological replicates and three technical replicates were used for each experiment. The reaction mixture (50 µL) for the PCR amplification contained 25 µL of 2 × Taq Plus Master Mix II (Dye Plus) (Vazyme Biotech Co., Ltd., Nanjing, China), 2 µL each of 10 µmol/L forward and reverse primers, 17 µL of ddH 2 O, and 4 µL of DNA (30 ng). The PCR products were sequenced by Sangon Biotech (Shanghai, China).

Morphological Characteristics of Male Fertility and Sterility
Male-fertility phenotypes of the CMS line HZ1A and its maintainer line HZ1B were observed during their flowering ( Figure 1A,B). HZ1A had the typical features of male-sterility, including no pollen dissemination around its anthers. The pollen grains of the HZ1A plants were shriveled, less in number, and stained blue-green with Alexander's stain. The pollen grains of the HZ1B plants were plump and stained magenta red with Alexander's stain. method was used for calculating the relative mRNA expression. Three biological replicates and three technical replicates were used for each experiment. The reaction mixture (50 µL) for the PCR amplification contained 25 µL of 2 × Taq Plus Master Mix II (Dye Plus) (Vazyme Biotech Co., Ltd., Nanjing, China), 2 µL each of 10 µmol/L forward and reverse primers, 17 µL of ddH2O, and 4 µL of DNA (30 ng). The PCR products were sequenced by Sangon Biotech (Shanghai, China).

Morphological Characteristics of Male Fertility and Sterility
Male-fertility phenotypes of the CMS line HZ1A and its maintainer line HZ1B were observed during their flowering ( Figure 1A,B). HZ1A had the typical features of malesterility, including no pollen dissemination around its anthers. The pollen grains of the HZ1A plants were shriveled, less in number, and stained blue-green with Alexander's stain. The pollen grains of the HZ1B plants were plump and stained magenta red with Alexander's stain.

The Cytological Characteristics of CMS Line HZ1A and Its Maintainer Line HZ1B Microspores during the Five Developmental Events
The length of flower buds showed a significant correlation with different stages of pollen development in pepper. The corresponding relationship in the HZ1B was as follows: about 1.5 mm, the sporogenous tissue stage; about 2.3 mm, the meiotic stage; about 3.5 mm, the tetrad period; about 4.5 mm, the mononuclear stage; and about 6.5 mm, the pollen mature stage (Figure 2A).
There was no differences between the microspores of the CMS line HZ1A and its maintainer line HZ1B during the sporogenous tissue stage and meiotic stages ( Figure  2B,C). In the tetrad stage, few microspores in the HZ1A were similar to that in the HZ1B in terms of their cytological features; however, few microspores (indicated by red arrows) in the HZ1A were anomalous and irregular. It was difficult to find microspores in the  The length of flower buds showed a significant correlation with different stages of pollen development in pepper. The corresponding relationship in the HZ1B was as follows: about 1.5 mm, the sporogenous tissue stage; about 2.3 mm, the meiotic stage; about 3.5 mm, the tetrad period; about 4.5 mm, the mononuclear stage; and about 6.5 mm, the pollen mature stage (Figure 2A).
There was no differences between the microspores of the CMS line HZ1A and its maintainer line HZ1B during the sporogenous tissue stage and meiotic stages ( Figure 2B,C). In the tetrad stage, few microspores in the HZ1A were similar to that in the HZ1B in terms of their cytological features; however, few microspores (indicated by red arrows) in the HZ1A were anomalous and irregular. It was difficult to find microspores in the mononuclear stage of the HZ1A as the number of microspores was less. These results indicate that the male sterility in HZ1A to begin at the tetrad stage of microspore development.

RNA Sequencing and DEG Identification
To clarify the basic molecular mechanism regulating male sterility in the CMS line HZ1A, gene expression levels in the HZ1A and HZ1B anthers during three developmental stages (S1, S2, and S3) were compared using RNA-Seq with the PromethION platform (ONT, Oxford, UK). RNA sequencing of S1_HZ1A, S2_HZ1A, S3_HZ1A, S1_HZ1B, S2_HZ1B, and S3_HZ1B, resulted in a total of 47.95 million clean reads being obtained from the raw reads ( Table 1). The number of clean reads from the six samples ranged from 7.29 to 8.52 million reads. These reads had a mean length of 711-923 bp. The data for the six samples were on an average 84.04-87.91% that of the full-length reads. The mapped rate ratios of the six samples to the pepper reference genome were 76.83-87.45% on an average. This indicated that the aligned data were suitable for subsequent analyses.
Horticulturae 2021, 7, x FOR PEER REVIEW 6 of 18 mononuclear stage of the HZ1A as the number of microspores was less. These results indicate that the male sterility in HZ1A to begin at the tetrad stage of microspore development.

RNA Sequencing and DEG Identification
To clarify the basic molecular mechanism regulating male sterility in the CMS line HZ1A, gene expression levels in the HZ1A and HZ1B anthers during three developmental stages (S1, S2, and S3) were compared using RNA-Seq with the PromethION platform (ONT, Oxford, UK). RNA sequencing of S1_HZ1A, S2_HZ1A, S3_HZ1A, S1_HZ1B, S2_HZ1B, and S3_HZ1B, resulted in a total of 47.95 million clean reads being obtained from the raw reads (Table 1). The number of clean reads from the six samples ranged from 7.29 to 8.52 million reads. These reads had a mean length of 711-923 bp. The data for the six samples were on an average 84.04-87.91% that of the full-length reads. The mapped rate ratios of the six samples to the pepper reference genome were 76.83-87.45% on an average. This indicated that the aligned data were suitable for subsequent analyses. There was a total of 59,468 functionally annotated transcripts. The reads were assembled into 31,603 unigenes, including 35,336 genes that were matched to nuclear reference sequences, 190 genes that were matched to mitochondrial reference genome sequences, 84 genes to chloroplast reference genome sequences, and 4958 new genes. The number of genes ranged from 25,144-25,814 across all six samples ( Figure 3A). The gene expression levels were compared at each developmental stage of the anthers between the CMS line HZ1A and its maintainer line HZ1B, and a total of 1377 genes across all stages were identified as being differentially expressed (Table S2). A total of 42 (27 up-regulated and 15 downregulated), 691 (346 up-regulated and 345 down-regulated), and 709 (281 up-regulated and 428 down-regulated) DEGs were obtained from stages S1, S2, and S3 for the CMS line HZ1A, respectively ( Figure 3C). Only four DEGs were common in all three stages ( Figure 3D). More up-regulated DEGs were detected in the CMS line HZ1A at the S1 stage. In contrast, more down-regulated DEGs were detected in the CMS line HZ1A at the S2 and S3 stages.

Functional Annotation by GO
A total of 1377 DEGs from S1, S2, and S3 were classified into three main categories: 'biological process (BP)', 'cellular component (CC)', and 'molecular function (MF)'. These three categories were composed of 46 functional groups based on GO assignments ( Figure  4). In the 'BP' category, 'metabolic process', 'cellular process' and 'single-organism process' were the major functional groups, including a total of 277, 276, and 223 genes in the three stages, respectively. The top six significant GO terms were different in the three stages. Photosynthesis (GO:0015986) was significantly enriched in the S1 and S3 stages (Table S3). In the 'CC' category, 'membrane', 'membrane part', 'cell part', 'cell' and 'organelle' were the dominant functional types, including a total of 308, 283, 272, 272, and 210 genes in three stages, respectively. The top six significant GO terms in the three stages are shown in Table S4. Photosystem I (GO:0009522), endomembrane system (GO:0012505),

Functional Annotation by GO
A total of 1377 DEGs from S1, S2, and S3 were classified into three main categories: 'biological process (BP)', 'cellular component (CC)', and 'molecular function (MF)'. These three categories were composed of 46 functional groups based on GO assignments (Figure 4). In the 'BP' category, 'metabolic process', 'cellular process' and 'single-organism process' were the major functional groups, including a total of 277, 276, and 223 genes in the three stages, respectively. The top six significant GO terms were different in the three stages. Photosynthesis (GO:0015986) was significantly enriched in the S1 and S3 stages (Table S3). In the 'CC' category, 'membrane', 'membrane part', 'cell part', 'cell' and 'organelle' were the dominant functional types, including a total of 308, 283, 272, 272, and 210 genes in three stages, respectively. The top six significant GO terms in the three stages are shown in Table S4. Photosystem I (GO:0009522), endomembrane system (GO:0012505), and endoplasmic reticulum (GO:0005783) were significantly enriched in the S2 and S3 stages. In the 'MF' category, 'binding' and 'catalytic activity' were the dominant functional types, including a total of 453 and 449 genes in the three stages, respectively. The top six significant GO terms in the three stages are shown in Table S5. Hydrogen ion transmembrane transporter activity (GO:0015078), serine-type endopeptidase activity (GO:0004252), and O-methyltransferase activity (GO:0008171) were significantly enriched in S1, S2, and S3. Interestingly, the photosynthesis (GO:0015986) term in the 'BP' category and photosystem I (GO:0009522) term in the 'CC' category were both involved in photosynthesis. Therefore, the photosynthesis process may be related to CMS.

Pathway Mapping by KEGG
Functional annotation and pathway enrichment analysis of the DEGs was performed using the KEGG database. DEGs were enriched in 13, 100, and 97 (Tables S6-S8) pathways in the three stages, respectively, as observed by the KEGG pathway analysis, respectively. Only two pathways, oxidative phosphorylation (ko00190) ( Figure 5A) and phenylpropanoid biosynthesis (ko00940) ( Figure 5B) were significantly enriched in the S1 and S2 stages, respectively. Although the numbers of DEGs in the S2 and S3 stages were almost the same, there was no significantly enriched KEGG pathway. There were four DEGs in the oxidative phosphorylation pathway and 24 DEGs in the phenylpropanoid biosynthesis pathway (ko00940). In the latter pathway, almost all the DEGs (20) were upregulated in the CMS anthers. This indicated that the up-regulation of phenylpropanoid biosynthesis might cause sterility in the pepper CMS line, HZ1A.

Identification of Differentially Expressed TFs
The TFs of the DEGs were analyzed at the S1, S2, and S3 stages of the CMS line HZ1A and its maintainer line HZ1B. A total of 72 TFs were identified by predicting all the DEGs, and 1 (1 up-regulated), 52 (34 up-regulated and 18 down-regulated), and 21 (15 upregulated and 6 down-regulated) TFs were identified in S1, S2, and S3 stages of the CMS line HZ1A, respectively (Table S9). These TFs were classified into 30 families, including NAC (9), LOB (8), bZIP (5), MADS-M-type (4), MYB (4), and 25 other TFs ( Figure 6). Only three TFs (Capana02g000496, Capana02g002747, and Capana06g001532) were up-regulated in the S2, and S3 stages simultaneously. There were 45 TFs that were classified into 25 families that were up-regulated in the CMS line, whereas 27 TFs that were classified into 12 families that were down-regulated in the CMS line. In the S1 stage, only one TF (Capana02g000758) was up-regulated in the CMS line. Most of the TFs were up-regulated or down-regulated during the S2 and S3 developmental stages. These results suggested that there is a complex transcriptional network was involved in the anther development.

Validation of the DEGs Using Quantitative Real-Time PCR
Lipid metabolism and endoplasmic reticulum-related genes are associated with male sterility [29,37]. To validate the reliability of the RNA-Seq, nine DEGs related to lipid metabolism (Capana00g002971, Capana04g001148, Capana07g002334, Capana08g001266, and Capana12g002508), and endoplasmic reticulum (Capana01g000067, Capana06g000385, Ca-pana06g002873, and Capana09g001928) were randomly selected and used for qRT-PCR analysis. The expression levels of the nine DEGs showed similar trends in both RNA-Seq analysis and qRT-PCR results (Figure 7). These results indicate that the RNA-Seq results were reliable. . Gene Ontology (GO) enrichment analysis results for the differentially expressed genes in anthers of the cytoplas mic male-sterility line HZ1A and its maintainer line HZ1B during S1, S2, and S3 stages. There were 46 functional group in three categories.

Pathway Mapping by KEGG
Functional annotation and pathway enrichment analysis of the DEGs was perfo using the KEGG database. DEGs were enriched in 13, 100, and 97 (Tables S6-S8) path in the three stages, respectively, as observed by the KEGG pathway analysis, respect Only two pathways, oxidative phosphorylation (ko00190) ( Figure 5A) and pheny panoid biosynthesis (ko00940) ( Figure 5B) were significantly enriched in the S1 a stages, respectively. Although the numbers of DEGs in the S2 and S3 stages were a the same, there was no significantly enriched KEGG pathway. There were four DE the oxidative phosphorylation pathway and 24 DEGs in the phenylpropanoid biosy sis pathway (ko00940). In the latter pathway, almost all the DEGs (20) were up-regu in the CMS anthers. This indicated that the up-regulation of phenylpropanoid biosy sis might cause sterility in the pepper CMS line, HZ1A. . Gene Ontology (GO) enrichment analysis results for the differentially expressed genes in anthers of the cytoplasmic male-sterility line HZ1A and its maintainer line HZ1B during S1, S2, and S3 stages. There were 46 functional groups in three categories.

Mitochondrial Genes Involved in Energy Production Are Related to CMS
Mitochondria play an important role in energy production. Mitochondrial deficiencies may lead to failure in acquiring sufficient energy for normal pollen development [62]. In some studies, the sterility-related genes in plant mitochondria, such as orf138, orf79, atp4, and atp6 [8], cause abnormal energy supplies and eventually leading CMS. In the present study, RNA sequencing data were compared to the mitochondrial reference genome and the chloroplast reference genome. In total, 12 mitochondrial DEGs were obtained. Of these DEGs, of which three (orf470, atp6 and cox1) were associated with energy production. Comparison of the gene expression levels from the RNA-Seq analysis suggested that the genes orf470 and cox1 in the CMS line HZ1A showed lower expression than those in the HZ1B during the S1 and S3 stages, respectively; however, the atp6 gene showed higher expression levels in HZ1A than in HZ1B during the S1 and S2 stages. Primers (Table S1) for these three DEGs were designed to amplify in the CMS line HZ1A and its maintainer line HZ1B. The results showed that orf470 was amplified only in the HZ1B; however, the other two DEGs were effectively amplified ( Figure 8A). According to the sequencing results of the PCR products, there was no difference in the sequences of cox1 between the HZ1A and HZ1B; however, several nucleotide differences in the sequence of atp6 were observed between the HZ1A and HZ1B, which lead to changes in the amino acid sequences of corresponding proteins ( Figure 8B). These results indicate that orf470 and atp6 may be related to male sterility and were thus considered the candidate genes for male sterility in the CMS line, HZ1A. This will require further verification through functional characterization and other methods in the future.

Identification of Differentially Expressed TFs
The TFs of the DEGs were analyzed at the S1, S2, and S3 stages of the CMS line and its maintainer line HZ1B. A total of 72 TFs were identified by predicting all the and 1 (1 up-regulated), 52 (34 up-regulated and 18 down-regulated), and 21 (15 up lated and 6 down-regulated) TFs were identified in S1, S2, and S3 stages of the CM HZ1A, respectively (Table S9). These TFs were classified into 30 families, including (9), LOB (8), bZIP (5), MADS-M-type (4), MYB (4), and 25 other TFs ( Figure 6). Only TFs (Capana02g000496, Capana02g002747, and Capana06g001532) were up-regula the S2, and S3 stages simultaneously. There were 45 TFs that were classified into 25 lies that were up-regulated in the CMS line, whereas 27 TFs that were classified i families that were down-regulated in the CMS line. In the S1 stage, only o (Capana02g000758) was up-regulated in the CMS line. Most of the TFs were up-reg or down-regulated during the S2 and S3 developmental stages. These results sug that there is a complex transcriptional network was involved in the anther develop

Validation of the DEGs Using Quantitative Real-Time PCR
Lipid metabolism and endoplasmic reticulum-related genes are associated with male sterility [29,37]. To validate the reliability of the RNA-Seq, nine DEGs related to lipid metabolism (Capana00g002971, Capana04g001148, Capana07g002334, Capana08g001266, and Capana12g002508), and endoplasmic reticulum (Capana01g000067, Capana06g000385, Capana06g002873, and Capana09g001928) were randomly selected and used for qRT-PCR analysis. The expression levels of the nine DEGs showed similar trends in both RNA-Seq analysis and qRT-PCR results (Figure 7). These results indicate that the RNA-Seq results were reliable.

Mitochondrial Genes Involved in Energy Production Are Related to CMS
Mitochondria play an important role in energy production. Mitochondrial defi cies may lead to failure in acquiring sufficient energy for normal pollen development In some studies, the sterility-related genes in plant mitochondria, such as orf138, o atp4, and atp6 [8], cause abnormal energy supplies and eventually leading CMS. In

The Anther Was the Best Organ for Studying Male Sterility in Pepper
In plants, the phenomenon of male fertility is reflected in the anthers and poll grains, and the development of pollen grains is generally related to the anthers. The a thers are consequently the best organs in which to study male sterility. To study ge expression related to male fertility in pepper, directly study of the anther and poll grains would be the best practice.
Peppers have small anthers, especially in their early stages of development, whi makes it difficult to obtain the anther samples separately for RNA extraction across all t developmental stages. Therefore, the pepper floral buds have been used in previous stu ies for transcriptome analysis [62,63]. In wheat and rice, RNA stabilization solutions ha been used to collect anthers and pollen grains, and qualified RNA has been extracted [6 Therefore, the RNA stabilization solution (RNALater) was used to collect the pepper a thers and extract the RNA successfully in the present study. On the other hand, accordin to the results of the expressions of genes in pepper floral buds and anthers using qRT-PC (not published data), the expressions of some genes were different in the flower buds an anthers. This result was the same as in soybean [65]. Consequently, the anther was co sidered the best organ to study male sterility.

Male Sterility Occurs alongside Gene Expression Level Changes in the Middle and Late Stages of Anther Development
In each of the anther development stages assessed (S1, S2, and S3), the number detected genes was more than 25,000, and there was little difference between the sample This shows that a large number of genes are involved in the anther development of pepp

The Anther Was the Best Organ for Studying Male Sterility in Pepper
In plants, the phenomenon of male fertility is reflected in the anthers and pollen grains, and the development of pollen grains is generally related to the anthers. The anthers are consequently the best organs in which to study male sterility. To study gene expression related to male fertility in pepper, directly study of the anther and pollen grains would be the best practice.
Peppers have small anthers, especially in their early stages of development, which makes it difficult to obtain the anther samples separately for RNA extraction across all the developmental stages. Therefore, the pepper floral buds have been used in previous studies for transcriptome analysis [62,63]. In wheat and rice, RNA stabilization solutions have been used to collect anthers and pollen grains, and qualified RNA has been extracted [64]. Therefore, the RNA stabilization solution (RNALater) was used to collect the pepper anthers and extract the RNA successfully in the present study. On the other hand, according to the results of the expressions of genes in pepper floral buds and anthers using qRT-PCR (not published data), the expressions of some genes were different in the flower buds and anthers. This result was the same as in soybean [65]. Consequently, the anther was considered the best organ to study male sterility.

Male Sterility Occurs alongside Gene Expression Level Changes in the Middle and Late Stages of Anther Development
In each of the anther development stages assessed (S1, S2, and S3), the number of detected genes was more than 25,000, and there was little difference between the samples. This shows that a large number of genes are involved in the anther development of pepper at different developmental stages, however, the number of DEGs was significantly different at different developmental stages. In the S1 stage, only 42 DEGs were identified. While in the S2 and S3 stages, the number of DEGs (691 and 709), respectively, increased rapidly, which was similar to previous results for eggplant anthers [66] and pepper buds [62]. This shows that the change in gene expression is not obviously caused by male sterility in the early stage of the anther development, and that male sterility mainly affects the gene expression level in the S2 and S3 stages. These results were in-line with the male sterility beginning at the tetrad stage (S2) of the CMS line HZ1A microspore development. This suggests that male sterility causes changes in the gene expression levels in the middle and late stages of anther development in the CMS line HZ1A.

TFs Related to the CMS in Pepper Line HZ1A
TFs are associated with fertility in some plant species. GAMYB encodes a TF associated with the gibberellin pathway and leads to abnormal meiosis through regulating another development in rice [32,33]. CSA encodes an R2R3 MYB TF that leads to male sterility through reduced levels of carbohydrates in rice [67]. A barley PHD TF, MALE STER-TILITY1 (MS1) is expressed in the anther tapetum and plays a critical role during pollen development [68]. The barley HvMS1 gene encodes a PHD-finger TF that is expressed in the anther tapetum, which is essential for pollen development and causes complete male sterility when overexpressed in barley [68]. The TF ZmDREB1.7 partially restores male fertility of CMS-S maize by activating the expression of the mitochondria-encoded CMS gene orf355 [69]. The SlMS10 gene (Solyc02g079810) encodes the bHLH TF that regulates meiosis and cell death of the tapetum during microsporogenesis in tomato [70].
In our study, the highest number of TFs was observed in the S2 stage, and only one TF was identified in the S1 stage. These results were in-line with the male sterility in the CMS line HZ1A beginning at the tetrad stage (S2) of microspore development. The TFs in the S2 stage are most likely related to male sterility. We analyzed the TFs that belong to the same family of the above described TFs in HZ1A and HZ1B. Only bHLH (two TFs) and MYB (one TF) families were identified in the S2 stage. These three TFs (Capana01g000416, Capana01g000418, and Capana02g000991) might regulate fertility in HZ1A.

Genes Involved in Lipid Metabolism Are Related to Male Sterility
Male sterility is a phenomenon caused by a complex process that occurs widely in flowering plants. Lipid transport and metabolism are related to male fertility development in the anther. In the present study, according to the functional annotation of DEGs, 104 lipid metabolism-related genes were found to be differentially expressed, and 1 (1 down-regulated), 68 (46 up-and 22 down-regulated), and 36 (18 up-and 18 downregulated) genes were identified in stages S1, S2, and S3 of the CMS line HZ1A, respectively. In the S2 stage, the number of lipid metabolism-related genes was higher than that in the S1 and S3 stages. These results indicated that S2 is an important stage in the lipid metabolism pathway.
The lipid metabolism-related genes were analyzed using KEGG pathway enrichment. The significantly enriched terms were fatty acid metabolism (ko01212, 10 genes), biosynthesis of unsaturated fatty acids (ko01040, 5 genes), sphingolipid metabolism (ko00600, 5 genes), and glycerophospholipid metabolism (ko00564, 6 genes) ( Figure S1). To date, many lipid metabolism-related male sterile genes have been identified and characterized, and the transcriptional regulation pathways, lipid metabolism, and control of male fertility have been well investigated in plants [29,71]. Therefore, these four pathways in our study may play important roles in regulating male sterility in the CMS line HZ1A during anther development.
There were 21 DEGs in the four KEGG pathways. Of them, 10 DEGs were downregulated in the CMS line, HZ1A and these DEGs may be related to male sterility. Among these 10 DEGs, Capana02g003412 was down-regulated and showed the lowest expression (log 2 (fold change) = −2.88) in the CMS line HZ1A. According to the KEGG pathway annotation, Capana02g003412 was a glycosphingolipid biosynthesis-related gene. Glycosphingolipids are the major constituents of the outer leaflet of the plasma membrane in eukaryotic cells [72]. On the other hand, we searched the homologous genes using the amino acid sequence of Capana02g003412 through NCBI-Blast (https://blast.ncbi.nlm.nih.gov/Blast.cgi, accessed on 30 October 2021). One of the homologous genes, Arabidopsis thaliana gene AT5G20710, was found to be related to pollen development [73]. These results suggest that Capana02g003412 may be related to pollen development and male sterility.

Conclusions
In the present study, a comparative transcriptome analysis of the anthers was conducted between the CMS line HZ1A and its maintainer line HZ1B during three developmental stages. RNA-Seq analysis identified that 42, 691, and 709 DEGs were obtained in stages S1, S2, and S3, respectively. GO and KEGG analysis of the DEGs revealed that two GO terms and one KEGG pathway may be related to CMS. The mitochondrial genes orf470 and atp6 were identified as candidate genes for CMS. Therefore, our study elucidated the pathways and key genes related to the CMS of pepper, and our results provided insights into the molecular mechanisms of CMS, laying the foundation for further research in pepper.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/horticulturae7120580/s1, Figure S1: Scatterplot of the KEGG pathway enrichment results for the lipid metabolism-related DEGs between HZ1A and HZ1B anthers, Table S1: Primers of qRT-PCR and mitochondrial genes in the present study, Table S2: DEGs in HZ1B vs. HZ1A anthers during three stages, Table S3: The top six significant GO terms in 'biological process' during three stages, Table S4: The top six significant GO terms in 'cellular component' during three stages, Table S5: The top six significant GO terms in 'molecular function' during three stages, Table S6: The pathways of the DEGs at S1 stage based on KEGG metabolic pathways, Table S7: The pathways of the DEGs at S2 stage based on KEGG metabolic pathways, Table S8: The pathways of the DEGs at S3 stage based on KEGG metabolic pathways, Table S9: The transcription factors from the DEGs in the S1, S2, and S3 development stages.