Exploring differentially expressed genes of microspore embryogenesis under heat stress in sweet pepper

Stress is considered to be the inducer of microspore embyogenesis (ME), and heat stress is indispensible in the ME of sweet pepper. The aim of the study was to explore differentially expressed genes of microspore embryogenesis under heat stress in sweet pepper. The swollen rate of microspore was significantly affected by heat stress, while no green plant could be acquired without heat pretreatment. Anthers with or without heat stress were used for whole transcriptome analysis by RNA sequencing to provide new insights on how cells adapt to stress. A total of 5031 differentially expressed genes were identified, among which 2657 differentially expressed genes were up-regulated and 2374 differentially expressed genes were down-regulated in the early stage of heat stress. KEGG pathway analysis identified "plant hormone signal transduction" (67; 11.20%), followed by starch and sucrose metabolism (63; 10.54%). RNA-Seq data and quantitative real-time polymerase chain reaction showed that 224 genes related to glutathione metabolism, starch and sucrose metabolism, plant hormone signal transduction and phenylpropanoid biosynthesis were the most likely specific genes in ME under heat stress. This research provides new insights into molecular regulation during the early stage of ME in sweet pepper under heat stress.


INTRODUCTION
Stress is considered as a kind of mutation on microspore embryos: unconstrained microspores form flower powder along the normal binding pathway (Touraev et al., 1997). Among various stress treatments, heat stress is widely used to initiate microspore embyogenesis (ME) in many crops (Asadi et al., 2018;Bhatia et al., 2018;Cimò et al., 2017;Dubas et al., 2014). The anthers of sweet pepper cultured under heat stress can be induced to develop into haploid embryos with complete functions instead of mature pollen (Bárány et al., 2005). Despite extensive advances have been made, compared to cruciferous and cereal species, sweet pepper is still considered recalcitrant to ME and doubled haploid production, which limit the use of this technology in breeding programs (Seguí-Simarro et al., 2011). In the past few years, most molecular biology researchers have focused on the culture system of sweet pepper (Popva et al., 2016;Heidari-Zefreh et al., 2019;Sánchez et al., 2020). We do not know the mechanism of cell fate transformation and the regulation of gene expression that initiates *Corresponding author. E-mail: chengyan820620@163.com.
Author(s) agree that this article remain permanently open access under the terms of the Creative Commons Attribution License 4.0 International License embryogenesis, and this knowledge gap is particularly obvious (Tsuwamoto et al., 2007). Many biochemical and morphological changes are closely related to the forced manipulation and alteration of gene expression patterns in embryos (Soriano et al., 2013).
Today, transcriptome analysis is an important means to study possible mechanisms and identify potential genes. Using this method, the roles of stress to induce the embryogenesis and alternation of gene expression have been examined in several crops (Elhiti et al., 2012;Liu et al., 2016;Zhang et al., 2019). However, it is still unclear for the differentially expressed genes in the early stage of ME under heat stress in sweet pepper.
In our previous work, an efficient ME system for pepper was established (Cheng et al., 2013). The present study aimed at exploring differentially expressed genes of microspore embryogenesis under heat stress in sweet pepper using high-throughput sequencing technology. This discovery will provide new insights into the molecular mechanism of sweet pepper micro-ecosystem.

Plant material and heat treatment
Sweet pepper (Capsicum annumm L.) variety Jinjiao 203, a responsive genotype in ME, grew up in the greenhouse of Shanxi Academy of Agricultural Sciences, China. In February 2017, seeds were sown in the soil, and flower buds were collected from the donor plants. Anther pretreatment was performed as previously decribed by Cheng et al. (2013). Briefly, flower buds were sterilized immediately after collection. They were dissected and cultured on pretreatment medium (10 mM CaCl2, 1 mM MnSO4·7H2O, 1 mM KNO3, 200 M KH2PO4, 1 M KI, 100 nM CuSO4·5H2O, 0.37 M mannitol and 0.5% agar) for embryogenesis induction. Anthers cultured at 25°C were used as control group and anthers cultured at 34°C as heat treatment group. Each group was repeated four times. Anther load of each biological replica was isolated from at least five feed plants and divided into two groups on average. After 7 days of pre-culture, the microspores of 5 anthers were separated in 1 ml sterile water, and the microspore expansion rate (the total expanded microspores and the fraction of the total microspores in 5 microscope fields) of each treatment was analyzed. At the same time, a small amount of anthers (about 100 mg) were collected into test tubes by liquid nitrogen quick freezing method, and then stored at -80°C or RNA isolation. As described by Cheng et al. (2013), the remaining anthers were isolated and cultured for embryogenesis. After 2 months of culture, the number of green plants in every 100 anthers was counted.

RNA isolation and cDNA library construction
In order to separate RNA from pretreated anthers, frozen samples were ground into fine powder in a microcentrifuge tube. Total RNA was isolated by triazole reagent (Invitrogen, Carlsbad, California, USA), which produced about 10 micrograms of total RNA per sample. The nanodrop 2000 spectrophotometer (Thermo) was used to detect the concentration, and the RNA nano6000 detection kit of Agilent Bioanalyzer 2100 system (Agilent Technologies, California, USA) was used to evaluate the integrity. After adjusting their concentration to 10 nm, appropriate total RNA samples were collected in the same volume in 4 repeated pretreated anthers.
A total of 1 microgram of ribonucleic acid was used as input material in each sample to prepare a ribonucleic acid sample. Following the manufacturer's recommendations, the NEBNext UltraTM RNA Library Preparation Kit for Illumina (NEB USA) was used to create a sequencing library and the index code to the attribute sequence of each sample was added. Briefly, mRNA was purified from total RNA using magnetic beads and addition of oligonucleotide poly-T. At high temperature, divalent cations was used to adhere to the buffer of NEBNext synthesis reaction (5×). The first cDNA strand was synthesized with 6 polysilicon random seeds and mulv arz. Then deoxyribonucleic acid polymerase I and ribonuclease H were used to synthesize the second strand of deoxyribonucleic acid. The remaining overhangs were converted to blunt ends by exonuclease/polymerase activity. After adenylation at the 3' end of deoxyribonucleic acid fragment, NEBNext adapter with hairpin loop structure was connected to prepare for hybridization. To select a cDNA fragment with a length of 240 bp, the library fragment was purified by AMPure XP system (Beckman Coulter, Beverly, USA). Then, before polymerase chain reaction, 3 mU·L -1 of user enzyme (National Biological Laboratory, USA) was used to treat the cDNA ligated with the adapter selected in size at 37°C for 15 min, and then treated at 95°C for 5 min. Then polymerase chain reaction was carried out with polymerase, universal primer and index primer. Finally, the PCR products (Ample XP system) were purified on Agilent Bioanalyzer 2100 system, and the library quality was evaluated.

RNA sequencing and transcript analysis
According to the manufacturer's instructions, TruSeq PE Cluster Suite v4-cBot-HS (Illumia) was used to cluster the index coded samples on cBot cluster generation system. After the cluster is created, the preparation order of libraries on the Illumina Hiseq Xten platform is specified, and generates an end reading pair. The raw data in Fastq format (raw reading) was first processed by internal perl script. This step deletes the adapter-containing read, policy-containing read and low-quality read from the original data to obtain clean data (clean read). Q20, Q30, GC-content and sequence iteration level of clean data are also calculated. All downstream analysis is based on high quality and clean data. The adaptor sequences and low quality sequence reads would be removed from the dataset. The original sequence was converted to clean reading after data processing. These clean reads were then mapped to the reference genome sequence database Zunla-1. Only completely consistent or inconsistent readings are further analyzed and annotated according to the reference genome. Tophat2 tool software is used to map reference genomes. The Kyoto Encyclopedia of Genes and Genomes predicts metabolic and cellular pathways. The gene expression level was calculated by reading millions of Millennium bug (FPKM), and the fragments used were mapped to reference sequences. The two groups of DEGs were screened according to gene expression levels using the DESeq R software package (1.10.1). Binary negative sign distribution is used to identify differential expression in digital genetic expression data. More than twice the expression level changes and significantly different expressions (P<0.05) are considered as the differential expressions among different treatments.

Quantitative real-time PCR analysis
To verify the expression of SDR, 20 candidate genes identified by KEGG concentration analysis were randomly selected for real-time PCR quantitative analysis (QR-PCR). Selected gene names and primer information are listed in Table 1. Before assembling RNA sequences, the total RNA (1 μg) of each sample was used as a

Gene id Forward (5'→3') Rewerse (5'→3')
CATCTGAGGCTACTTGGTGTC GTTCCAGTTCATGCTTCCATTG template for synthesizing the first strand cDNA. Gene names and information of selected seeds are shown in Table 1. Before compiling the nucleic acid sequence, the whole DNA [1 of each sample (1 microgram)] was used as a model for synthesizing the first cDNA strand. On the Bio-Rad CFX96 instrument, trans start top green qpcr super mix (AQ131) was used to replicate each biological sample for three times, and quantitative reverse transcription polymerase chain reaction analysis was carried out Using capsicum actin nanny gene as internal control. Quantitative comparison CT method ( ΔΔ CT method) is used to quantify the relative expression of specific genes (Livak and Schmittgen, 2001).

Effects of heat stress on the induction of ME in sweet pepper
In order to check the influence of thermal stress on electromagnetic induction, 7 days 34°C anther pretreatments were processed. In all four biological replicates, the anther volume increased after heat stress compared with the control group ( Figure 1A). Moreover, more microspores in anthers were inflated, while fewer expanded microspores and no green plants were observed in the control ( Figure 1B). The results indicated that both the swollen rate of microspore and green plants production were significantly affected by the heat stress (34°C, 7 days) in anther pretreatment, which means that heat stress could improve the embryogenesis initiation and possibility of microspore to develop into green plants.

Illumina sequence analysis and verification of selected genes by quantitative reverse transcription polymerase chain reaction
The cDNA library of embryogenic microspores of sweet pepper was sequenced on Illumina Hiseq Xten platform, and paired terminal readings were generated. A total of 143,541,722 and 151,420,406 original readings were obtained from two cDNA libraries composed of control group and thermal pressure group (Table 2). 92.91% of the readings were useful after data quality inspection and screening with quality greater than 30 (NQ30), among which 89.88% (128,990,381) of the control group readings and 82.90% (125,760,669) of the heat stress group readings were located in the pepper genome. About 86.39% of the readings were located on genes, and 13.61% were not located on genes, which indicated that most of the readings were located on reference genes. In order to confirm the differentially expressed genes confirmed in sequencing and computational analysis, 20 DEGs concentrated in KEGG pathway were randomly selected for quantitative real-time polymerase chain reaction including 2 signal transduction mechanisms protein (Capana00g002630, Capana01g000883), 1 probable glutathione S-transferase (Capana00g003106), 1 agamous-like MADS-box protein (Capana01g002457), 1 peroxidase of carbohydrate transport and metabolism (Capana00g001129), 1 ent-copalyl diphosphate synthase Compared with the control group, the single asterisk and double star numbers showed statistical differences (T test, P < 0.05 or 0.01, respectively)  (Capana01g000278), 1 protein P21 of thaumatin family (Capana01g000279), 1 hydrophobic seed protein (Capana01g000500), 1 probable lipid transfer protein (Capana01g000731), 1 pollen-specific leucine-rich repeat elongation protein 1 (Capana01g003124), 1 C1-like domain (Capana01g003125), 1 serum kinase protein /suaminotransfer LRR receptor (Capana00g004543), 1 shikimate O-hydroxycinnamoyltransferase (Capana01g004373), and 3 function unknown protein (Capana00g005078, Capana01g001352, Capana01g001414).
The results of quantitative reverse transcription polymerase chain reaction (Figure 2) on the expression of 20 selected genes showed that compared with the control group, the average expression level of 12 genes in the heat stress group was significantly increased. However, the average gene expression level of heat stress group was significantly down-regulated. The expressions of these genes were consistent in RNA RNA-Seq and qRT-PCR data.

KEGG enrichment analysis on DEGs
In addition to 32832 unchanged genes, 2657 upregulated gene expression regions and 2374 downregulated gene expression regions were found in the microenvironment between the control group and the hot group. In order to further identify the possible functional pathways, 598 DEGs were found in 20 KEGG functional pathways by KEGG pathway analysis. Among these pathways, the category of 'plant hormone signal transduction' (67; 11.20%) represented the largest group, followed by 'starch and sucrose metabolism' (63; 10.54%), 'biosynthesis of amino acids' (61; 10.20%), 'Carbon metabolism' (57; 9.53%), 'phenylpropanoid biosynthesis' (51; 8.53%) and 'Plant-pathogen interaction' (50; 8.36%). KEGG enrichment analysis showed that heat stress regulatory genes were enriched in four main pathways (p<0.05), including 'glutathione metabolism', starch and sucrose metabolism', 'plant hormone signal transduction' and 'phenylpropanoid biosynthesis' ( Figure  3). Among the four main pathways, 224 estimated heat stress specific genes have been confirmed, of which 138 genes are up-regulated under heat pressure and 86 genes are down-regulated under heat pressure (Supplementary Table 1).

DISCUSSION
High stress tolerance could be presumed the first prerequisite for successful ME induction, for the reason that ME was induced or at least strongly stimulated by a stress pre-treatment (Zoriniants et al., 2005). One of the most important components of stress resistance is an effective antioxidant system composed of enzymes and low molecular antioxidants, which protects cells from the production of reactive oxygen species (Mittler, 2002). Glutathione was a major antioxidant in all forms of life and an indicator of cellular oxidative stress (Belmonte and Stasolla, 2009). In a reduced form, glutathione was metabolized in multiple ways leading to the biosynthesis of mercapturonate, glutamate, glycine, cysteine and other amino acids (Noctor et al., 2012). It also participates in the regulation of cell cycle, proliferation and programmed cell death as a signal molecule. Our results proved the importance of glutathione metabolism under heat stress in ME of sweet pepper. The key role of glutathione in embryo and meristem development was confirmed by analyzing the phenotype of glutathione-deficient Arabidopsis mutants (Noctor et al., 2012). The published results show that some glutathione response genes have encoded transcription factors and proteins, and are involved in cell division, redox potential (such as thioredoxin, glutamoren), auxin biosynthesis, transport and regulation of transcription reaction (Schnaubelt et al., 2013). It revealed that the initial environment of embryogenesis requires a reduced environment (high GSH/GSH + GSSG), which may promote cell proliferation by enhancing nucleotide synthesis and mitotic activity (Stasolla, 2010). More and more recently published data indicate that ROS accumulation initiates signal transduction leading to microspore reprogramming and embryogenic development (Żur et al., 2019).
In this study, starch and sucrose metabolism were detected as the second major pathways in ME of sweet pepper under heat stress. The metabolism of starch and sucrose fuels all aspects of plant growth and development. The kinetics of starch synthesis is related to the differentiation process and the clear change of cell wall structure and organization, which is characterized by the de-esterification of pectin and the increase of RGII and XG components (Bárány et al., 2005;Satpute et al., 2005). The accumulation of plastid starch occurred in differentiated cells, which showed high levels of deesterified pectin and rich RGII and XG, while the proliferation cells rich in esterified pectin did not show starch deposition (Bárány et al., 2010). Data show that one of the important changes of carbohydrate Daisy network is related to the transformation of embryogenesis and development program and cell fate (Corral-Martínez et al., 2019). The results reported in this paper show that many genes related to cell wall, major carbohydrate and starch metabolism are expressed differently during early embryonic development.
It was found that among the 20 KEGG functional pathways, plant hormone signal transduction accounted for the largest number of DEGs. Plant growth regulators are considered as the core signaling molecules for controlling plant growth and development, which respond to environmental stimuli and initiate signaling pathways (Kohli et al., 2013). Moreover, plant growth regulators interfere with the interaction between plant genotypes and environmental factors, and play a very important role in the micro-ecosystem, controlling the differentiation and development of embryos derived from microspheres and the regeneration of haploid/diploid plants (Divi et al., 2010). In the early stage of ME, numerous genes related to axins (Dubas et al., 2014), cytokinins, abscisic acid (Dubas et al., 2013), gibberellins, brassinosteroids, jasmonic acid, salicylic acid, or ethylene may produce translational regulation to adapt the stress. Hormone homeostasis seems to be one of the most important factors that determine the embryogenesis ability of cells, and a more comprehensive method is needed to understand the mechanism controlling this process, so as to break the barrier of self-resistance (Żur et al., 2015).
Phenylpropanoid biosynthesis was induced by several stresses (Dixon and Paiva, 1995), and became the main pathway in ME of sweet pepper under heat stress. In plants, the phenylpropanoid pathway was responsible for the synthesis of secondary metabolites, including lignin monomers, flavonoids, and coumarins, which play essential roles in determining plant structure, biomass recalcitrance, and stress tolerance (Vogt, 2010). It had been the significantly enriched pathways with specific common DEGs mainly including peroxidase, phenylalnine ammonialyase, and β-glucosidase (Zhang et al., 2019).
The up-regulated DEGs related to phenylpropanoid biosynthesis in this study might have important function to keep microspore from death and promote embryonic microspore dedifferentiation.

Conclusions
Heat stress is indispensible in the ME of sweet pepper. In this study, transcriptome analysis of the anther of sweet pepper in the early stage of ME showed that the DEGs between heat treatment and control were mainly associated with glutathione metabolism, starch and sucrose metabolism, plant hormone signaling and phenylpropyl ester biosynthesis. Our research provides new insights into molecular regulation during the early stage of ME in sweet pepper under heat stress.