Comparative Transcriptome Analysis Reveals Different Low-Nitrogen-Responsive Genes in Pepper Cultivars

The molecular mechanisms underlying the variation in N-use efficiency (NUE) in pepper (Capsicum annuum L.) genotypes are poorly understood. In this work, two genotypes (750-1, low-N tolerant; ZCFB, low-N sensitive) with contrasting low-N tolerance were selected from 100 pepper cultivars on the basis of their relative leaf areas, shoot dry weights, root dry weights, and plant dry weights at the seedling stage. Subsequently, using RNA-Seq, the transcriptome of these two pepper genotypes under N starvation for 28 days was analyzed. We detected 2621/2470 and 3936/4218 different expressed genes (DEGs) in the leaves/roots of 750-1 and ZCFB, respectively. The changes in the expression of basic N metabolism genes were similar between 750-1 and ZCFB. However, different DEGs not directly involved in N metabolism were identified between the 750-1 and ZCFB cultivars. In 750-1, 110 unique DEGs were detected in the leaves, of which 103 were down-regulated, including genes associated with protein metabolism, photosynthesis, secondary metabolism, cell wall metabolism, stress response, and disease resistance. In ZCFB, 142 unique DEGs were detected in the roots, of which 117 were up-regulated, resulting in enhancement of processes such as protein degradation, secondary metabolites synthesis, lipid metabolism, endocytosis, the tricarboxylic acid cycle (TCA), transcriptional regulation, stress response, and disease resistance. Our results not only facilitate an understanding of the different regulatory process in low-N-tolerant and low-N-sensitive pepper cultivars, but also provide abundant candidate genes for improving the low-N tolerance of pepper cultivars.


Introduction
Nitrogen (N) is one of the most important macronutrients for plants, and it can be absorbed and assimilated by the roots in various forms, including nitrate, ammonium, and amino acids. Nitrate and ammonium are the most common forms used by plants, with nitrate being the dominant form [1]. N metabolism can be divided into three processes: uptake, assimilation, and remobilization. In the uptake stage, at least six transporters participate in nitrate uptake in Arabidopsis (NPF6.3/NRT1.1, NPF4.6/NRT1.2, NRT2.1, NRT2.2, NRT2.4, and NRT2.5) [2,3]. In the assimilation process, nitrate is reduced to ammonium by nitrate and nitrite reductases, and then the ammonium is assimilated into amino acids by glutamine synthetase (GS), glutamine antinotransferase (GOGAT), and asparagine synthetase (AS) [4]. The N remobilization process comprises protein degradation and amino acid transport. Several genes involved in the ubiquitin-26S proteasome pathway [5] and encoding amino acid transporters [6,7] have also been shown to be associated with these processes. These studies have provided important clues for understanding the mechanism of N metabolism. However, present knowledge about the complicated N regulatory network remains incomplete, and the molecular basis governing the genetic variation of N-use efficiency (NUE) among crop cultivars remains unclear.
Pepper (Capsicum annuum L.) is one of the most important vegetable crops of the Solanaceae family and is grown worldwide as a food, medicine source and ornamental plant [8]. It is well documented that, within a certain range, there is a positive correlation between the relative growth rate of pepper and the N concentration in the soil. However, inappropriate N utilization can lead to undesirable growth, yield, and quality [9,10]. In recent years, in areas with fertile soil, farmers have been using a higher amount of N fertilizer than is required, which has resulted in low NUE and serious environmental pollution [11]. By contrast, in some areas where the soil quality is poor, the need for a high input of N fertilizer represents a great economic burden on pepper growers [12]. Therefore, improving NUE is urgently needed for the development of sustainable pepper production. Genetic variation in NUE has been reported for different crops such as rice [13], barley [14], wheat [15], rapeseed [16], maize [17], and cotton [18], but molecular knowledge about the genetic variation of NUE is still very poor. Several genes responsible for improved NUE have been identified in rice, including DEP1 [19], OsNRT1.1B [20], OsNRT2.3b [21], ARE1 [22], OsNRT1.1A [23], and OsNPF6.1 [24]. Unfortunately, studies on the molecular regulation mechanism of N metabolism in pepper are scarce, let alone those studies focusing on NUE variation.
The next-generation high-throughput RNA sequencing technology (RNA-Seq) is a powerful tool for revealing genome-wide changes under biotic/abiotic stresses and can provide system level information regarding the N metabolism network. RNA-Seq analysis has been applied to the transcriptome analysis of low-N response of a single plant genotype, such as those of cucumber [25], maize [26], wheat [27], physic nut [28], and rice [29]. A large number of candidate genes involved in low-N response were detected [25,[27][28][29]. Furthermore, the potential regulatory roles of IncRNAs in response to N stress have also been investigated [26]. However, it is difficult to reveal variations in NUE using only one genotype. Therefore, comparative transcriptome analysis of genotypes with different low N tolerances has become more recognized as a tool for understanding NUE [30][31][32][33][34][35]. A high abundance of transcripts related to high affinity nitrate transporters (NRT2.2, NRT2.3, NRT2.5, and NRT2.6) in the N-stress tolerant sorghum genotypes [31] and an energysaving assimilation pattern in N-stress tolerant Tibetan wild barley genotype have been revealed [32].
The focus of this study was to identify low-N-responsive genes that were differentially expressed between low-N-tolerant and low-N-sensitive genotypes after long-term low-N stress, with the aim of providing more information on NUE variation. First, we examined the low-N tolerance of 100 pepper cultivars. Subsequently, we analyzed the genome-wide gene expression changes of two pepper cultivars with contrasting low-N tolerance under low-N stress. Lastly, the different low-N-responsive genes between the two cultivars were intensively analyzed using RNA-Seq and various bioinformatics methods, revealing potential new candidate genes for improving the low-N tolerance of pepper cultivars.

Screening of Pepper Cultivars with Different Low-N Tolerance
The low-N tolerance of 100 pepper cultivars was analyzed hydroponically in a greenhouse at the Chongqing Key Laboratory of Adversity Agriculture. Yamazaki nutrient solution for pepper was used as a sufficient-N solution containing 1. 50 PO 4 were 20% that of the sufficient-N solution, and the other nutrients were the same as the sufficient-N solution. Furthermore, K 2 SO 4 , KH 2 PO 4 , and CaCl 2 were added in moderation to the low-N solution to avoid K, P, and Ca deficiency. The rest of the seedlings grew under N-sufficient condition as controls. Each cultivar had three replicates (15 plants/replicate) for each N condition (sufficient-N or low-N). The culture solution was refreshed every seven days at pH 6.2-6.4. All plants were grown in a greenhouse at 28 • C (day) and 25 • C (night) with a relative humidity of 70-80% and a photoperiod of 14/10 h, 300 µmol m −2 s −1 . The leaf areas, shoot dry weights, root dry weights, and plant dry weights were evaluated 28 days later.

RNA Preparation and Sequencing
Seeds from the 750-1 (low-N-tolerant) and ZCFB (low-N-sensitive) pepper cultivars were germinated, and the seedlings were cultured under two different N conditions (sufficient-N and low-N) as described above. Four weeks later, the youngest fully expanded leaves and roots of the seedlings grown in sufficient-N solution and low-N solution were collected separately at 10:00-11:00 A.M., flash frozen in liquid N, and then stored at −80 • C. Samples were collected from 45 independent plants from three replicates (15 plants/replicates). Total RNA was extracted using an RNA Plant Kit (Aidlab Biotech, Beijing, China). The quantity and quality of total RNA was examined on a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Madison, WI, USA) and an Agilent 2100 Bioanalyzer (Waldbronn, Germany), respectively. cDNA libraries were constructed and sequenced on a BGISEQ-500 platform at BGI (Shenzhen, China).

Transcriptome De Novo Assembly, Gene Functional Annotation, and Differentially Expressed Genes (DEGs) Analyses
The raw reads of transcriptome sequencing were filtered using SOAPnuke v.15.2 (BIG, Shenzhen, China), and a set of clean reads was obtained. After de novo assembling, the clean reads were mapped to the pepper reference genome [36]. Genes were functionally annotated based on the NCBI non-redundant (Nr) [37], Gene Ontology (GO) [38], and the Kyoto Encyclopedia of Genes and Genomes (KEGG) databases [39]. For gene expression analysis, the numbers of matched reads were calculated and then normalized to RPKM by RSEM v.1.2.12. Significant differential expression genes (DEGs) were identified as those with a fold-change ≥ 2.0 and an FDR ≤ 0.001. DEGs were clustered using GO-Term Finder software, and pathway enrichment analysis was performed based on terms from the KEGG database.

Quantitative Real-Time PCR (qRT-PCR) Analysis
To verify the reliability of RNA-Seq data, total RNA was extracted from the leaves and roots of biological triplicates from the two sequenced genotypes, and the template cDNA samples were prepared using an iScrip First Strand Synthesis System Kit (Bio-Rad Laboratories, Hercules, CA, USA). Primers for each PCR reaction were designed to have a melting temperature of 58-62 • C and to produce a PCR product between 100 and 200 bp (Table S1). The pepper actin gene (AY572427) was used as an internal control. qRT-PCR reactions were performed using a CFX96 TM Real-Time System (Bio-Rad Laboratories, Hercules, CA, USA). The reaction conditions were as follows: 95 • C for 30 s, 40 cycles of 95 • C for 5 s, and 60 • C for 30 s. The 2 −∆∆Ct method was used to calculate the expression levels of genes.

Screening for Pepper Cultivars with Contrasting Low-N Tolerance
To identify pepper cultivars with contrasting low-N tolerance, 100 pepper cultivars were screened based on various metrics. To avoid the effects of natural variation in individual embryos and/or endosperm size on initial growth from different genotypes, the relative value (low-N treated sample/control) of each character was used in the following analyses. Our results showed that peppers from these 100 cultivars had a wide range of relative values of leaf areas, shoots, roots, and plant dry weights, and the coefficient of variation (CV) of all of the relative values was above 15%, demonstrating the high variation in tolerance to low-N among cultivars (Table 1 and Table S2). Cluster analysis was then carried out based on the relative values of leaf areas, shoot dry weights, root dry weights, and plant dry weights. The 100 pepper cultivars could be grouped into three clusters, consisting of five (Group I), six (Group II), and 89 (Group III) cultivars ( Figure 1). There were significant differences in each characteristic between the three groups ( Table 2). The mean relative values of each characteristic in Group I were significantly higher than those of Group II, and the mean relative values of each characteristic in Group III were between those of Groups I and II, indicating that Group I was a low-N-tolerant group, while Group II was a low-N-sensitive group. 750-1 had the highest values for the relative leaf area as well as shoot and plant dry weights (Table S2) and was classified into Group I. Thus, it was considered low-N tolerant. ZCFB had the lowest relative values of shoot, root and plant dry weights (Table S2) and was classified into Group II. This cultivar was therefore considered low-N sensitive. Cluster analysis of the low-nitrogen tolerance from 100 pepper cultivars.Cluster analysis was carried out based on the relative values of leaf area and shoot, root, and plant dry weights. Agglomerative hierarchical clustering was performed using the hclust algorithm in R, with dst = dist (dat1, method = "euclidean") and hclust (dst, method = "average"). Group I, low-N-tolerant group; Group II, low-N-sensitive group; the low-nitrogen tolerance of Group III was between that of Groups I and II.

RNA-Seq Analysis of Two Pepper Cultivars with Contrasting Low-N Tolerances
To better understand the molecular mechanism governing the genetic variation of N-use efficiency (NUE) in pepper cultivars, the global gene expression profile of two pepper genotypes (750-1 and ZCFB) with contrasting low-N tolerance was next analyzed using RNA-Seq. Samples were taken 28 days after low-N treatment, and the following samples were used for sequencing: TCL1 (750-1-leaf control), TCR2 (750-1-root control), TTL3 (750-1-leaf low-N treatment), TTR4 (750-1-root low-N treatment), SCL5 (ZCFB-leaf control), SCR6 (ZCFB-root control), STL7 (ZCFB-leaf low-N treatment), and STR8 (ZCFBroot low-N treatment). Sequence data are available from the NCBI SRA database with accession number PRJNA700470. Over 65 million clean reads were generated from each library. A total of 91.75-93.55% of these clean reads could be well aligned to the pepper reference genome, at least 69% could be uniquely mapped to the reference genome (71.88% on average). The number of genes with at least one mapped transcript were 29,167,29,817,28,970,29,981,29,330,29,814,29,078, and 30,335 for libraries TCL1, TCR2, TTL3, TTR4, SCL5, SCR6, STL7, and STR8, respectively (Table 3).   To validate our RNA-Seq data, the expression levels of eight genes (LOC107841845, LOC107855963, LOC107858815, LOC107867964, LOC107870519, LOC107875600, LOC107875602, and LOC107875603) in the leaves and roots of 750-1 and ZCFB were next confirmed by qRT-PCR. The expression trends of the eight genes were consistent with our results from RNA-Seq in the eight libraries ( Figure 3).

GO and KEGG Enrichment Analysis of DEGs
Next, GO analysis was carried out using the identified DEGs in the three main GO categories: biological process, cellular component, and molecular function ( Figure 4A). KEGG pathway analysis showed that the DEGs involved in N metabolism were enriched in all of the four libraries ( Figure 4B). The DEGs involved in nitrogen metabolism were then intensively analyzed. The DEGs associated with nitrate transportation, assimilation, and remobilization processes were detected in each genotype ( Figure 5 and Figure S1A-D). The genes of the nitrate transporter family were up-regulated, especially in the roots, while N assimilation systems were suppressed. The genes participating in N remobilization were also enhanced. Overall, the expression patterns of the DEGs involved in N metabolism were similar in both genotypes. Thus, further analysis is needed to find clues governing different NUEs between the two genotypes examined.

Different DEGs between Two Cultivars
To examine the different DEGs between the two cultivars, comparisons were made between the DEGs of 750-1 and ZCFB ( Figure 6A,B). The absolute values of log2 foldchanges were set as equal to or greater than five. Furthermore, the DEGs involved in N metabolism, which were analyzed previously, were excluded from these analyses. The results indicated that 110 unique DEGs from 750-1 could be detected in leaves, and most were down-regulated ( Figure 6A), while 142 unique DEGs from ZCFB were detected in roots, and most were up-regulated ( Figure 6B). Most of the unique DEGs are described in the following section ( Figure 7A-J).

Carbon and Lipid Metabolism
In the carbon metabolism process, genes associated with photosynthesis were downregulated in the leaves of 750-1 ( Figure 7B). Two genes (LOC107843574 and LOC107857137) of the tricarboxylic acid cycle (TCA) were upregulated in the leaves and roots of ZCFB. In the lipid metabolism process, DEGs involved in lipid synthesis and degradation were inhibited in the leaves of 750-1, while those in the roots of ZCFB were significantly enhanced, including members of the 3-ketoacyl-CoA synthase (LOC107845580, LOC107847631, and LOC107868466), fatty acid desaturase (LOC107851810), and phospholipase (LOC107859298) families ( Figure 7C).

Secondary Metabolism
The genes responsible for the biosynthesis of phenylpropanoid, sesquiterpenoid, diterpenoid, flavonoid, anthocyanin, and other secondary metabolites showed altered expression under N starvation in both genotypes, including members of the cytochrome P450 family (LOC107844024, LOC107860971, LOC107849698, LOC107871074, LOC107867086, LOC107845094, LOC107877698, LOC107856388, and LOC107874613) and anthocyanin biosynthesis (LOC107845761 and LOC107866550) ( Figure 7D). The flavonoid biosynthesis process consists of the largest number of DEGs followed by the phenylpropannoid biosynthesis process. All the unique DEGs in the leaves of 750-1 were repressed, while most of the unique DEGs in the roots of ZCFB were up-regulated.

Cell Wall Synthesis and Structure
In the leaves, all of the unique DEGs in this category were repressed in both of the genotypes, including xyloglucan 6-xylosyltransferase 2 (LOC107842350), which is responsible for cell wall synthesis in ZCFB, as well as Glucan edo-1,3-beta-glucosidase (LOC107879143) and beta-D-xylosidase (LOC107854300), which functions in cell wall degradation in 750-1 ( Figure 7E). Most of these genes in the roots of the two genotypes were up-regulated, including probable pectate lyase 5 (LOC107848186 and LOC107848183) in 750-1 and putative pectinesterase 11 (LOC107869412) and polygalacturonase-like (LOC107856362 and LOC107869773) in ZCFB.

Transport
In 750-1, all of the unique DEGs engaged in transport were down-regulated in the leaves, including genes associated with osmosis (LOC107865680) and stoma movement (LOC107844000), while in the roots, except for the gene encoding sodium-dependent phosphate transporter protein 1 (LOC107857112), all of the genes were up-regulated, including metal (LOC107869634) and sugar transport (LOC107873879) ( Figure 7F). In ZCFB, most of the genes were up-regulated in both of the leaves and roots, including genes associated with channel (LOC107855209), endocytosis (LOC107848190), sodium/metabolite and urea transport (LOC107871496 and LOC107870432) in the leaves, and gene encoding channel proteins (LOC107840638) and proteins involved in endocytosis (LOC107848190, LOC107840435, and LOC107845831) in the roots.

Signal Transduction and Transcription Factor
For hormone synthesis and regulation, all of the genes were inhibited in the leaves of both genotypes ( Figure 7G), including those involved in ethylene synthesis (LOC107840369) and response (LOC107865651) in 750-1 and those associated with the auxin response (LOC107840512 and LOC107853594) in ZCFB. In 750-1, the genes of the LRR receptor-like serine/threonine-protein kinase family (LOC107869144, LOC107869171, and LOC107869351) and those involved in Ca2+ signaling (LOC107847049 and LOC107844087) were up-regulated in the leaves and roots, respectively. The genes of transcription factors (TFs) that showed changes in expression mainly fell into the MYB, bHLH, NAC, WRKY, BEE, and ERF families ( Figure 7H).

Stress Response and Disease Resistance
In 750-1, most of the unique DEGs related to stress response and disease resistance were observable in the leaves, and except for protein SRG1-like (LOC107859182) and NBS-LRR root-knot nematode resistance protein (LOC107864520), all the DEGs of the two pathways were down-regulated ( Figure 7I,J). Meanwhile, in ZCFB, DEGs involved in stress response and disease resistance were detected in both the leaves and roots, and most of the genes in the leaves were down-regulated, in contrast to in the roots.

Discussion
In this study, a large number of low-N-responsive genes were identified from low-Ntolerant (750-1) and low-N-sensitive (ZCFB) pepper genotypes after long-term N deficiency, and the number of low-N-responsive genes in ZCFB was far greater than those in 750-1 ( Figure 2). The expression profiles of the genes associated with nitrate transport, assimilation, and remobilization were similar between the 750-1 and ZCFB ( Figure 5). However, many unique DEGs were detected between the two pepper genotypes after long-term N deficiency, implying different molecular mechanisms of resistance to N deficiency in these two genotypes.

Unique DEGs in the Low-N-Tolerant Genotype
In 750-1, the unique DEGs involved in primary metabolism, secondary metabolism, stress response, and disease resistances were down-regulated in the leaves after long-term N deficiency. However, the relative biomass of 750-1 was significantly higher than that of the other cultivars (Table S2), indicating that the growth of 750-1 was repressed less than other cultivars under N starvation. The relative values of root dry weights and leaf areas of 750-1 were much higher than those of other cultivars (Table S2), implying a larger absorption area of N and stronger photosynthesis in 750-1. Interestingly, genes of somatic embryogenesis receptor kinase 2 (SERK2)-like (LOC107857123) and LRR receptorlike serine/threonine-protein kinase (LOC107869144 and LOC107869171) were up-regulated in the leaves of 750-1 ( Figure 7G). Both SERK2 and LRR receptor-like serine/threonineprotein kinase belong to the receptor-like kinase (RLK) family. SERK2 is essential for male microsporogenesis in Arabidopsis [40]. It has been postulated that at an early low-N stress stage, processes such as absorption, transportation, and assimilation of N might be highly enhanced in 750-1 under low-N conditions, as well as primary metabolism, which could result in improved biomass accumulation in 750-1 relative to other cultivars. After longterm N limitation, 750-1 might slowly decrease vegetative growth and enter reproductive stage earlier to accelerate its life cycle and seed handing down.

Unique DEGs in the Low-N Sensitive Genotype
Lipids are a major subcellular component, the biosynthesis and composition of which are influenced by N [41,42]. As essential components of the membrane, lipids play an important role in endocytosis, in which plasma membrane lipids and associated proteins are internalized in vesicles that fuse with the endosomal system [43]. In our study, a number of DEGs involved in lipid metabolism were up-regulated in the roots of ZCFB ( Figure 7C). Meanwhile, several genes associated with endocytosis were also enhanced ( Figure 7F). In plants, endocytosis plays an important role not only in transporting membrane proteins, lipids, and extracellular molecules into the cell but also in nutrient delivery, toxin avoidance, and pathogen defense [44,45]. To our knowledge, few studies have detected the involvement of endocytosis in N starvation. Here, the detection of endocytosis-related genes being involved in N starvation could provide new insights into N metabolism in plants.
Secondary metabolites have no direct role in plant growth, but their roles in stress defense have garnered much interest in recent years [46,47]. In the present study, genes involved in the synthesis of secondary metabolites, such as phenylpropanoid and flavonoid, were significantly up-regulated after low-N stress in ZCFB, as well as those encoding ATP-binding cassette (ABC) transporters that are involved in the transport of plant secondary metabolites ( Figure 7F) [48,49]. We also found that genes involved in the oxidation/antioxidation response ( Figure 7H) and disease resistance ( Figure 7I) were upregulated in ZCFB, as well as the expression of genes belonging to the NAC, ERF, MYB, and bHLH families ( Figure 7H), which play essential roles in stress response. Considering the low relative biomass of ZCFB (Table S2), enhanced secondary metabolism and stress response could suggest that there was a diversion of materials from primary metabolism to the defense system under N deficiency. Indeed, previous studies have shown that N limitation results in the reduction of plant growth; however, some secondary metabolites were also shown to accumulate [25,50]. Moreover, genes involved in the tricarboxylic acid cycle (TCA) were up-regulated in ZCFB ( Figure 7B), implying that to survive N deficiency, ZCFB may require more energy to support enhanced biological processes, such as endocytosis, secondary metabolism, and transcriptional regulation.

P and N Crosstalk
Phosphorus (P) is acquired by plants primarily in the form of inorganic phosphate (Pi) by Pi transporters (PHT). The Nitrogen limitation adaptation (NLA) gene is involved in adaptive responses to low-N stress in Arabidopsis [51] and is an important regulator of PHTs under N deficiency [52][53][54][55]. Over-expression of the rice Pi transporter gene OsPT2 enhances tolerance to Pi deficiency, as well as N2 fixation and ammonium assimilation [56,57]. These studies indicate the intimate crosstalk between the N limitation and P deficiency response pathways. In the present study, Sodium-dependent phosphate transport protein 1 (LOC107857112) was repressed in the roots of 750-1, and PHT1;2 (LOC107866602) and PHT1;9 (LOC107866766) were up-regulated in the roots of ZCFB under low-N stress ( Figure  7F), providing new information regarding the low-N induced crosstalk between N and P.

Conclusions
Improving NUE is urgently needed for the development of sustainable pepper production. However, present knowledge about the molecular basis governing the genetic variation of N-use efficiency (NUE) among pepper cultivars remains unclear. In this study, transcriptome of two pepper genotypes with contrary low-N tolerance (750-1, tolerant; ZCFB, sensitive) were compared using RNA-Seq and various bioinformatics methods after N starvation for 28 days. The results showed that the transcriptomic responses to low-N stress differed considerably between 750-1 and ZCFB, especially in genes that were not directly involved in N metabolism. The unique low-N-responsive genes between the two genotypes provide new insights for a comprehensive understanding of genotypic variation in NUE.

Data Availability Statement:
The data supporting the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no competing interest.