Transcriptome analysis of green and purple fruited pepper provides insight into novel regulatory genes in anthocyanin biosynthesis

Plant materials

Four pepper varieties, HN191, HN192, HN005 and EJT were grown in the experimental field of Chengdu Academy of Agriculture and Forestry Sciences (Chengdu, Sichuan, China). Erjingtiao (EJT) is a famous local pepper variety in Sichuan Province with a long cultivation history. It is characterized by its thin skin, thick flesh, delightful fragrance, moderate spiciness, and high nutritional value. This pepper is an essential ingredient in authentic Sichuan cuisine. The unripe fruits are green, while the ripe ones turn red. It is a very representative green fruit pepper material. HN191, HN192, and HN005 are inbred lines developed by our research group. These lines exhibit purple unripe fruits that transform into red mature fruits. In 2020, they were officially recognized by the Sichuan Provincial Non-Major Crops Committee.

Fruits from HN191 and EJT were harvested at 10, 20, 30, and 60 DAA, while fruits from HN192 and HN005 were harvested at 30 and 60 days after anthesis (DAA). For each variety and developmental stage, we collected three biological replicates and stored them at −80 °C for further experiments.

Assessment of anthocyanin

Anthocyanin mixtures from 30 DAA fruits of four pepper varieties (HN191, HN192, HN005, and EJT, respectively) were extracted and measured according to the agricultural standard of China (NY/T 2640-2014). The anthocyanin standard included delphinidin, cyanidin, petunidin, pelargonidin, peonidin, and malvidin, which were dissolved in methanol with 10% hydrochloric acid.

Microscope observation

Fresh fruits of HN191 at 30 DAA were cut to create free-hand sections and observed under a microscope (BX41; Olympus, Tokyo, Japan).

RNA extraction and sequencing

The pepper fruit samples that had been stored were removed from the −80 °C ultra-low temperature freezer and promptly placed into a foam box containing liquid nitrogen. Subsequently, the samples were taken out from the foam box and ground into powder using sterilized mortar and pestle in liquid nitrogen. The powdered samples were then transferred into 1.5 ml RNase-free microtubes (Corning Incorporated, Corning, NY, USA) immediately. Afterwards, total RNA was extracted using the RNAprep Pure extraction kit (Tiangen Biotech Co, Ltd., Beijing, China) in accordance with the manufacturer’s instructions. The purity and integrity of the RNA samples were evaluated using the Agilent Bioanalyzer 2100 system (Agilent Co., Ltd., Beijing, China). The sequencing library of all samples were constructed, and sequenced on the DNBSEQ-T7 platform in China National GeneBank (CNGB).

Data analysis

The adaptor and low-quality reads were first removed using Fastp (Chen et al., 2018). The resulting clean reads were then mapped to the pepper genome (Genbank assembly accession: GCA_000512255.2) using HISAT2 software (Kim, Langmead & Salzberg, 2015; Kim et al., 2014). To calculate the FPKM (fragments per kilobase of transcript per million mapped reads) for each gene, StringTie and Ballgown software were employed (Pertea et al., 2016).

The correlation between different samples was analyzed using the Pearson Correlation Coefficient (PCC) and Principal Component Analysis (PCA), which were implemented in the Factoextra and FactoMineR packages in R software.

For the differential expression analysis between samples, DESeq2 software was utilized (Varet et al., 2016). Genes with —log2Fold Change— ≥ 1 and a False Discovery Rate (FDR) < 0.05 were considered as differentially expressed genes (DEGs). The iTAK program was used to predict transcription factors from the DEGs (Zheng et al., 2016). Additionally, heatmaps of selected genes were generated using Tbtools software (Chen et al., 2020).

qRT-PCR analysis

The qPCR assay was configured following the recommendations of ‘The MIQE guidelines’ (Bustin et al., 2009). Ten structural genes and eleven regulatory genes (transcription factors, TFs) were selected to validate the RNA-seq results. Primers for qRT-PCR were designed using the Primer3 online software (http://bioinfo.ut.ee/primer3-0.4.0/) and synthesized by Sangon Biotech Co., Ltd., Shanghai, China (Table S1). For the cDNA synthesis, 1 µg of total RNA was reverse transcribed using the PrimeScript™ RT reagent Kit (Takara Bio Inc., San Jose, CA, USA), following the manufacturer’s protocol. Quantitative RT-PCR was performed on a CFX96 Real-Time PCR system (Bio-Rad Laboratories Inc., Hercules, CA, USA) using TB Green^® Premix Ex Taq™ II (Takara Bio Inc., Beijing, China). The consumables used include RNase-free tips and 8 Strip PCR tubes from Axygen^® Brand Products (Corning Incorporated, Corning, NY, USA). The quantification was performed in triplicate using 25 µL reactions. Each reaction included 12.5 µl of TB Green Premix Ex Taq II, 1 µl of each primer (10 µM), 8.5 µl RNase-free water, and 2 µl of 1:5 diluted cDNA. The PCR amplification conditions consisted of an initial denaturation at 95 °C for 30 s, followed by 40 cycles of denaturation at 95 ° C for 5 s and annealing at 60 °C for 30 s. A melting curve was obtained at the end of each PCR by gradually increasing the temperature to 95 °C (increment rates of 0.5 °C/s) after cooling to 65 °C for 5 s. Each gene was analyzed on the same amplification for all samples, so inter-run calibration was not necessary. The data obtained were analyzed by the Bio-Rad CFX Manager software (version 3.0), which generated the raw quantification cycle (Cq) values for each reaction. The relative expression of the selected genes was normalized to a pepper actin gene (GenBank ID: T459_30033) and calculated using the 2-ΔΔCT method. Further qPCR details are supplied in a MIQE checklist table (Table S2).

Phenotype analysis

In this study, four pepper varieties (HN191, HN192, HN005, and EJT) were selected as research materials. The unripe fruits of HN191, HN192, and HN005 have a purple color at 10, 20, and 30 days after anthesis (DAA), while the unripe fruits of EJT are green (Fig. 1). However, all four varieties have red ripe fruits at 60 DAA. The result of HPLC analysis showed that delphinidin was the only anthocyanin present in the three purple-fruited varieties (HN191, HN192, and HN005), while the levels of cyanidin, petunidin, pelargonidin, peonidin, and malvidin were below the detection threshold. No anthocyanins were detected in the green-fruited variety EJT. Furthermore, when we examined the fruits under a microscope, we observed that anthocyanins were only accumulated in the exocarp of HN191 fruits at 30 DAA (Fig. 2).

Global transcriptome results

We performed RNA-seq to investigate the dynamics of the transcriptome in fruit samples obtained from HN191 and EJT at four stages (10, 20, 30, and 60 DAA), as well as fruit samples from HN192 and HN005 at two stages (30 and 60 DAA). We analyzed three independent biological replicates of fruits for each variety and stage, resulting in a total of 36 samples. The extracted RNA exhibited a yield ranging from 14 to 52 µg and an RNA Integrity Number (RIN) ranging from 7.2 to 9.1, indicating relatively high RNA quality (Table S3).

A total of 7.25 billion clean reads, representing 725.48 GB of clean bases, were generated from all samples (Table S4) and were mapped to the pepper genome (Kim et al., 2014) at an average rate of 94.03% (range: 92.80%–94.80%). Power analysis on the transcriptome data was conducted using the RNASeqPower package (Hart et al., 2013). Based on a sample size of three (the number of biological replicates used in this study), the statistical power among different groups ranged from 0.85 to 0.91, with an average value of 0.87 (Table S5).

For each gene, the FPKM values were calculated, and genes with FPKM >1 were considered expressed. In EJT and HN191 at 10, 20, 30, and 60 DAA, as well as HN192 and HN005 at 30 and 60 DAA, a total of 19,303–21,826 genes were found to be expressed (Fig. 3A). The expression of 9%–11% of genes was very high (FPKM > 20) in the four pepper varieties at different developmental stages (Fig. 3B). Genes with high expression (FPKM 5 to ≤ 20), moderate expression (FPKM 1 to ≤ 5), and low expression (FPKM 0.1 to ≤ 1) accounted for 20%–23%, 34%–37%, and 29%–35% of the total genes, respectively.

Comparison of global transcriptome between samples

The global differences in transcriptomes during fruit development in HN191, HN192, HN005, and EJT varieties were investigated. The mean FPKM values of three biological replicates were used to calculate the Pearson correlation coefficients (PCC) for all samples (Fig. 4A). The PCC values between purple and green unripe fruits ranged from 0.69 to 0.88 (HN191-10 VS EJT-10, 0.88; HN191-20 VS EJT-20, 0.80; HN191-30 VS EJT-30, 0.69; HN192-30 VS EJT-30, 0.88; HN005-30 VS EJT-30, 0.88). However, for each variety, the PCC values between unripe fruits (10, 20, 30 DAA) and ripe fruits (60 DAA) were relatively low, measuring below 0.6.

Additionally, principal component analysis (PCA) was performed (Fig. 4B). The distances between samples in the plot indicate the similarity of their transcriptional programs. Samples of unripe fruits, both green and purple, were located on the left side of the plot, indicating that the transcriptional activity of purple and green unripe fruits is similar despite the difference in fruit color. On the other hand, samples of ripe fruits clustered on the right side, suggesting that unripe and ripe fruits are more likely to exhibit different transcriptomes and functions/activities, which aligns with the Pearson correlation coefficient results.

DEG analysis

In order to investigate genes associated with anthocyanin biosynthesis in pepper, our focus was on analyzing the differentially expressed genes (DEGs) between five comparisons of purple and green fruit varieties at 10, 20, and 30 days after anthesis (DAA). These comparisons included HN191 and EJT at 10 DAA, HN191 and EJT at 20 DAA, HN191 and EJT at 30 DAA, HN192 and EJT at 30 DAA, and HN005 and EJT at 30 DAA.

In total, we identified 14,734 upregulated genes, which included 674 genes encoding transcription factors (TFs), and 12,752 downregulated genes, which included 664 TF-encoding genes (Fig. 5A). The highest number of differentially expressed genes (7,012) was observed between HN192 and EJT at 30 DAA (HN192-30 VS EJT-30), followed by HN005 and EJT at 30 DAA (HN005-30 VS EJT-30). Conversely, the smallest number of differentially expressed genes (4,260) was found between HN191 and EJT at 10 DAA (HN191-10 VS EJT-10).

Common DEGs between purple and green fruits

The key distinguishing factor among these five comparisons is the fruit color, with green fruits being compared to purple fruits. Venn diagrams of the DEGs in each of the five comparisons were plotted using Venny2.1 (Fig. 5B). Among these comparisons, a total of 503 DEGs were identified that were common across all five, indicating their potential involvement in the biosynthesis of anthocyanin in pepper fruit.

Structure genes identified from common DEGs

Structural genes encode enzymes that directly participate in the biosynthesis of anthocyanin. From the pool of 503 common genes, we identified 7 structural genes involved in the anthocyanin biosynthesis pathway, namely CHS (T459_15321), F3H (T459_06228), F3′5′H (T459_28730), DFR (T459_06365), ANS (T459_00408), UFGT (T459_26178), and 3RT (MSTRG.12525), as shown in Fig. 6. The expression levels of these genes were relatively higher in purple fruits (HN191 at 10, 20, 30DAA, HN192 at 30 DAA, HN005 at 30DAA) compared to green fruits (EJT at 10, 20, 30 DAA). Moreover, their expression in purple fruits sharply decreased at 60 DAA, aligning with the observed pattern of anthocyanin accumulation. An additional important structural gene, CHI (T459_15095), was not part of the pool of 503 common DEGs because its fold change between HN191 and EJT at 10 DAA did not meet the threshold of 2. Nevertheless, in the remaining four comparisons, the expression levels of CHI in purple fruits were more than two times higher than in green fruits. For further reference, the sequences of these 8 structural genes can be found in Table S6.

TF prediction from common DEGs

Regulatory genes are responsible for encoding transcription factors that play a crucial role in the biosynthesis of anthocyanin (Lloyd et al., 2017). To identify these transcription factors, we utilized the iTAK software to predict them from the 503 common DEGs. As a result, a total of 24 transcription factors were identified (Table S7). We then generated a heatmap (Fig. 7) based on the FPKM values of these 24 transcription factors. Based on their expression patterns, these transcription factors can be divided into two distinct groups. The first group consisted of seven transcription factors (from T459_25295 to T459_26036) that exhibited relatively low expression in green fruits (EJT at 10, 20, 30 DAA), but showed higher expression in fruits of the three purple varieties (HN191 at 10, 20, 30 DAA, HN192 at 30 DAA, HN005 at 30 DAA). Conversely, the second group included 17 transcription factors (from T459_30722 to T459_18726) that displayed the opposite expression pattern. These findings suggest that the seven transcription factors in the first group may be closely associated with the biosynthesis of anthocyanin in pepper fruit.

qPCR validation

To validate the reliability of the RNA-seq data, we selected and analyzed 10 structural genes and 11 regulatory genes (transcription factors) associated with anthocyanin biosynthesis using qPCR.

The qPCR results demonstrated that the expression patterns of these selected genes were largely consistent with those observed in the RNA-seq data, except for one gene T459_08508 (Fig. 8, Table S8, and Table S9). The Pearson correlation coefficient between the RNA-seq and qRT-PCR data ranged from 0.54 to 0.99, with an average value of 0.88. This suggests that the RNA-seq data accurately reflect the abundance of transcript levels.