Excavation of Genes Responsive to Brassinosteroids by Transcriptome Sequencing in Adiantum flabellulatum Gametophytes

Brassinosteroids (BRs) are a class of polyhydroxysteroid plant hormones; they play important roles in the development and stress resistance of plants. The research on BRs has mainly been carried out in angiosperms, but in ferns—research is still limited to the physiological level and is not in-depth. In this study, Adiantum flabellulatum gametophytes were used as materials and treated with 10−6 M brassinolide (BL). The differentially expressed genes (DEGs) responsive to BRs were identified by transcriptome sequencing, GO, KEGG analysis, as well as a quantitative real-time polymerase chain reaction. From this, a total of 8394 DEGs were screened. We found that the expressions of photosynthetic genes were widely inhibited by high concentrations of BL in A. flabellulatum gametophytes. Moreover, we detected many BR synthase genes, except BR6ox2, which may be why castasterone (CS) rather than BL was detected in ferns. Additionally, we identified (for the first time) that the expressions of BR synthase genes (CYP90B1, CYP90C1, CYP90D1, CPD, and BR6ox1) were negatively regulated by BL in fern gametophytes, which indicated that ferns, including gametophytes, also needed the regulatory mechanism for maintaining BR homeostasis. Based on transcriptome sequencing, this study can provide a large number of gene expression data for BRs regulating the development of fern gametophytes.


Introduction
Brassinosteroids (BRs) are a class of polyhydroxysteroid plant hormones. Nowadays, more than 70 analogs have been found, among which brassinolide (BL) is the most active form [1][2][3]. BRs play important roles in regulating plant growth and development, such as promoting cell elongation, regulating plant fertility, and affecting photomorphogenesis [4,5]. The BR signal is transduced by the receptor-like protein kinases BRI1/BAK1-mediated pathway, in which BRI1 and BAK1 are the receptor and co-receptor of BRs, respectively [6,7]. So far, the research on BRs has mainly been carried out in angiosperms [8], but in ferns, research has been limited to the physiological level and is not in-depth. For example, BL with concentrations of 10 −7 M and 10 −6 M can inhibit spore germination of Pteridium aquilinum and gametophyte growth of A. flabellulatum, respectively [9,10], but the related genes are still unknown.
Castasterone (CS) is the immediate precursor of BL [11]. At present, CS, rather than BL, has been found in 12 species of Polypodiopsida (true ferns), such as P. aquilinum, Osmunda japonica, and Deparia japonica [12]. In angiosperms, some monocots, such as Oryza sativa and Zea mays, can only synthesize CS [12], but some dicots, such as Arabidopsis thaliana and Solanum lycopersicum (tomato), can synthesize BL from CS [13,14]. The studies show that the synthesis of CS from campesterol (CR) requires many enzymes, but CS can further synthesize BL only when BR6ox2 (C-6 oxidation) is present [15][16][17]. The above results indicate that BR6ox2 is presumably not present in ferns. Additionally, plants have the regulatory mechanisms for maintaining BR homeostasis [18]. For example, the expressions of BR synthase genes (CYP90B1 (DWF4), CYP90C1 (ROT3), CYP90D1, CPD, BR6ox1, and BR6ox2) are negatively regulated by BRs in A. thaliana [19,20]. Although many BR synthase genes are present in ferns [12], it is still unknown whether their expressions are negatively regulated by BRs.
Transcriptome sequencing has been widely used in the screening of differentially expressed genes (DEGs) and pathway analysis. Huang et al. [21], for example, treated Gerbera hybrida petals with BL, then used transcriptome sequencing, gene ontology (GO), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis to identify the DEGs, which were annotated to multiple plant hormone signaling pathways. Wu et al. [22] employed the methods above and discovered that BL increased the contents of total glucosinolates and sulforaphane by upregulating the gene expressions related to glucosinolate core pathways in Brassica oleracea. BL is the most commonly used drug in the study of BR regulating plant development, and it is also commonly used in O. sativa, Z. mays, P. aquilinum, and other plants that only synthesize CS [9,[23][24][25][26][27][28]. In view of this, 10 −6 M BL was applied to the A. flabellulatum gametophytes. The genes responsive to BL were identified and the BR synthase genes negatively regulated by BL were explored through transcriptome sequencing, GO, KEGG analysis, as well as a quantitative real-time polymerase chain reaction (qRT-PCR). So, this study can provide a large number of gene expression data for BR regulating the development of fern gametophytes.

Experimental Materials and Drugs
A. flabellulatum spores were collected from the rubber forest in Ma'an Mountain, Danzhou City, Hainan, China (109.519 • E, 19.501 • N), then stored at 4 • C in a dry, closed container. BL was purchased from Solarbio Company (CAS:78821-43-9), prepared into a 10 −3 M mother solution with 75% ethanol, and stored at 4 • C.

Culture of Adiantum flabellulatum Gametophytes
A. flabellulatum spores were immersed in 0.1% HgCl 2 to sterilize for 5 min and soaked in sterile water 5 times; they were then inoculated on 1/4 MS solid medium and cultured under light conditions (illumination was 500 Lx, the temperature was 24 ± 5 • C). After germination, young gametophytes with diameters of 2 mm were selected and transferred to a new 1/4 MS solid medium with different BL concentrations (0, 10 −12 , 10 −11 , 10 −10 , 10 −9 , 10 −8 , 10 −7 , and 10 −6 M), then cultured for 60 d under the illumination of 1000 Lx (24 h/d, 24 ± 2 • C). Nine biological replicates were set for each of the treatments. After photography, the perimeters and the projected areas of gametophytes in each treatment were measured by Image J (v1.8.0, Bethesda, Rockville, MD, USA), then the means and standard deviations were counted, and statistical significance analyses were performed by a least significant difference (LSD) test (SPSS 25.0, Armonk, NY, USA).

RNA Extraction, Library Construction, Sequencing, and Data Analysis
The gametophytes of A. flabellulatum grew slowly. In order to obtain enough sequencing materials, RNA extraction was performed from gametophytes cultured for 60 d in 0 M and 10 −6 M BL. Each treatment was set up with 3 biological replicates for a total of 6 samples. The CTAB method was used to extract the total RNA of each sample, then mRNAs were enriched by oligo (dT) magnetic beads. The enriched mRNAs were fragmented, and the cDNA libraries were constructed by reverse transcription. The BGISEQ-500 platform was used for sequencing. After removing raw reads with low quality, joint contamination, and unknown base N content greater than 5% by SOAPnuke (v1.4.0, Shenzhen, China), the clean reads were obtained [29]. The above steps were completed by the Beijing Genomics Institute (BGI, Beijing, China). Moreover, the clean reads of the 6 samples were uploaded to the SRA database of NCBI with the accession number PRJNA786072.

Quality Evaluation of Transcriptome Sequencing
The clean reads were aligned to the full-length transcripts (PRJNA733457) of A. flabellulatum gametophytes by Bowtie2 (v2.2.5, Baltimore, MD, USA) [30,31], and then a sequencing saturation analysis was performed. The expression levels of genes (calculated in fragments per kilobase of transcript per million mapped reads, FPKM) in each sample were counted by RSEM software (v1.2.8, Madison, WI, USA) [32]. According to the methods of Cai et al. [33], Pearson's correlation coefficient analysis, the principal component analysis (PCA), and a cluster analysis were performed on all samples.

Screening of Differentially Expressed Genes, GO, and KEGG Analysis
After removing the isoforms with average FPKMs < 0.5 in both treatments, the differentially expressed genes (DEGs) were screened under the conditions of |log 2 FoldChange| ≥ 1 and Q-Values < 0.001. According to the methods of Cai et al. [30], all DEGs were annotated and classified by GO (http://www.geneontology.org/ (accessed on 24 April 2022)) and KEGG (http://www.genome.jp/kegg/ (accessed on 24 April 2022)) analysis [30]. After that, enrichment of the analysis was performed by using the Phyper package in R software (Auckland, New Zealand). GO and KEGG annotations of the DEGs were completed by BGI.

Quantitative Real-Time Polymerase Chain Reaction
According to the methods of Cai et al. [33], 18 genes were selected for qRT-PCR (Supplementary Table S1). Three parallel tests were performed for each sample. The isoform_38906 (Actin) was used as the reference gene and the 2 -∆∆Ct method was used for relative quantification [34].

Quality Evaluation of Transcriptome Sequencing
Eight concentrations of BL were set up to observe their effects on the growth of A. flabellulatum gametophytes. The results showed that the perimeters and the projected areas of A. flabellulatum gametophytes were inhibited significantly when the concentrations of BL were more than 10 −7 M (Supplementary Figure S1A-C), which were consistent with the results by Wang et al. [10]. However, when the concentrations of BL were between 0 and 10 −8 M, the gametophytes changed irregularly (Supplementary Figure S1A). Therefore, transcriptome sequencing was performed on six samples from 0 M and 10 −6 M BL treatments (represented by CK and BL, respectively). The results showed that each sample generated at least 6.32 Gb data, with an average of 6.34 Gb data. The clean reads were obtained after quality filtering was performed on the raw reads, the acquisition rates were greater than 96%, wherein Q20 > 97% and Q30 > 92% (Supplementary Table S2). A sequencing saturation analysis showed that the increase of gene identification ratios in each sample tended to be flat when the read numbers were greater than 50 × 100 K, which indicated the sequencing data reached saturation (Supplementary Figure S2A). The cluster and correlation analyses were performed on the six samples. The results showed that three biological replicates in each treatment were clustered into one branch (Supplementary Figure S2B), and the Pearson correlation coefficients between the samples of each treatment were greater than or equal to 0.97 ( Figure 1A). Moreover, the principal component analysis (PCA) showed that CK and BL treatments were separated on PC1 ( Figure 1B). In conclusion, the six samples were of high quality with reliable sequencing data, which could be further analyzed.
Genes 2022, 13, x FOR PEER REVIEW 4 of 13 sample tended to be flat when the read numbers were greater than 50 × 100 K, which indicated the sequencing data reached saturation (Supplementary Figure S2A). The cluster and correlation analyses were performed on the six samples. The results showed that three biological replicates in each treatment were clustered into one branch (Supplementary Figure S2B), and the Pearson correlation coefficients between the samples of each treatment were greater than or equal to 0.97 ( Figure 1A). Moreover, the principal component analysis (PCA) showed that CK and BL treatments were separated on PC1 ( Figure 1B). In conclusion, the six samples were of high quality with reliable sequencing data, which could be further analyzed.

Screening of DEGs
A total of 210,415 isoforms were detected in the six samples, and the medians of expression levels were greater than 0.99 (Supplementary Table S3, Figure 2A). After removing the isoforms with average FPKMs of less than 0.5 in two treatments, 8394 differentially expressed genes (DEGs) were screened under the conditions of|log2 FoldChange| ≥ 1 and Q-Values < 0.001. Among them, 5081 (60.53%) DEGs were upregulated and 3313 (39.47%) DEGs were downregulated in the BL treatment (Supplementary Table S4, Figure 2B,C). The DEGs above were used for further analysis.

Screening of DEGs
A total of 210,415 isoforms were detected in the six samples, and the medians of expression levels were greater than 0.99 (Supplementary Table S3, Figure 2A). After removing the isoforms with average FPKMs of less than 0.5 in two treatments, 8394 differentially expressed genes (DEGs) were screened under the conditions of |log 2 FoldChange| ≥ 1 and Q-Values < 0.001. Among them, 5081 (60.53%) DEGs were upregulated and 3313 (39.47%) DEGs were downregulated in the BL treatment (Supplementary Table S4, Figure 2B,C). The DEGs above were used for further analysis.

GO Analysis of DEGs
The 8394 DEGs were aligned to the GO database, among which 5275 DEGs were annotated, accounting for 62.84% (Supplementary Table S5). GO terms were classified into three categories, including the biological process (BP), cellular component (CC), and molecular function (MF). Among them, MF annotated the most DEGs (4184). In BP, 2146 DEGs were annotated to the "cellular process", accounting for the highest proportion, followed by the "metabolic process" (1751). In CC, the DEGs were mainly classified into "cellular anatomical entity" (3495). In MF, most DEGs were annotated to "binding" (2635), followed by "catalytic activity" (2570). Additionally, in most GO terms, the upregulated DEGs in the BL treatment were in the majority (Supplementary Table S6, Figure 3A).

GO Analysis of DEGs
The 8394 DEGs were aligned to the GO database, among which 5275 DEGs were annotated, accounting for 62.84% (Supplementary Table S5). GO terms were classified into three categories, including the biological process (BP), cellular component (CC), and molecular function (MF). Among them, MF annotated the most DEGs (4184). In BP, 2146 DEGs were annotated to the "cellular process", accounting for the highest proportion, followed by the "metabolic process" (1751). In CC, the DEGs were mainly classified into "cellular anatomical entity" (3495). In MF, most DEGs were annotated to "binding" (2635), followed by "catalytic activity" (2570). Additionally, in most GO terms, the upregulated DEGs in the BL treatment were in the majority (Supplementary Table S6, Figure 3A). According to the GO enrichment results of DEGs, we showed the top 20 terms with low Q-values. Among them, "DNA binding" had the largest gene number (439), and the lowest Q-value, while "cilium" had the highest rich ratio (0.20). In addition, we found

GO Analysis of DEGs
The 8394 DEGs were aligned to the GO database, among which 5275 DEGs were annotated, accounting for 62.84% (Supplementary Table S5). GO terms were classified into three categories, including the biological process (BP), cellular component (CC), and molecular function (MF). Among them, MF annotated the most DEGs (4184). In BP, 2146 DEGs were annotated to the "cellular process", accounting for the highest proportion, followed by the "metabolic process" (1751). In CC, the DEGs were mainly classified into "cellular anatomical entity" (3495). In MF, most DEGs were annotated to "binding" (2635), followed by "catalytic activity" (2570). Additionally, in most GO terms, the upregulated DEGs in the BL treatment were in the majority (Supplementary Table S6, Figure 3A). According to the GO enrichment results of DEGs, we showed the top 20 terms with low Q-values. Among them, "DNA binding" had the largest gene number (439), and the lowest Q-value, while "cilium" had the highest rich ratio (0.20). In addition, we found According to the GO enrichment results of DEGs, we showed the top 20 terms with low Q-values. Among them, "DNA binding" had the largest gene number (439), and the lowest Q-value, while "cilium" had the highest rich ratio (0.20). In addition, we found gene numbers in "thylakoid", "photosynthesis", and "photosystem" were all large, which were 350, 275, and 242, respectively (Supplementary Table S7, Figure 3B).

KEGG Analysis of DEGs
The 8394 DEGs were aligned to the KEGG database, among which 3146 DEGs were annotated, accounting for 37.48% (Supplementary Table S8). KEGG pathways were classified into five categories-"metabolism", "genetic information processing", "environmental information processing", "cellular processes", and "organismal systems". Among them, "metabolism" annotated the most DEGs (1941) followed by "genetic information processing" (796). In "metabolism", 1629 DEGs were annotated to "global and overview maps", which had the largest gene number, followed by "carbohydrate metabolism" (614). In "genetic information processing", 301 DEGs were annotated to "translation" and 272 DEGs were annotated to "transcription". In "environmental information processing", the DEGs were mainly annotated to "signal transduction" (349). In "cellular processes", 212 DEGs were annotated to "transport and catabolism", and in "organismal systems", 182 DEGs were annotated to "environmental adaptation" (Supplementary Table S9, Figure 4A). sified into five categories-"metabolism", "genetic information processing", "environmental information processing", "cellular processes", and "organismal systems". Among them, "metabolism" annotated the most DEGs (1941) followed by "genetic information processing" (796). In "metabolism", 1629 DEGs were annotated to "global and overview maps", which had the largest gene number, followed by "carbohydrate metabolism" (614). In "genetic information processing", 301 DEGs were annotated to "translation" and 272 DEGs were annotated to "transcription". In "environmental information processing", the DEGs were mainly annotated to "signal transduction" (349). In "cellular processes", 212 DEGs were annotated to "transport and catabolism", and in "organismal systems", 182 DEGs were annotated to "environmental adaptation" (Supplementary Table S9  According to the KEGG enrichment results of DEGs, we showed the top 20 pathways with low Q-values. Among them, the pathways related to photosynthesis were "photosynthesis-antenna proteins", "photosynthesis", and "carbon fixation in photosynthetic organisms", and the DEG numbers were 113, 155, and 93, respectively. In addition, we found two pathways related to plant hormones, which were "plant hormone signal transduction" and "BR biosynthesis", and the DEG numbers were 203 and 33, respectively (Supplementary Table S10, Figure 4B).

The Expressions of Photosynthesis-Related Genes Were Inhibited
The KEGG analysis showed that the pathways of "photosynthesis-antenna proteins", "photosynthesis", and "carbon fixation in photosynthetic organisms" were enriched significantly, and the DEG numbers were all large. Therefore, we analyzed the DEG expressions in the above pathways.
The "photosynthesis-antenna proteins" pathway includes the LHCI complex and LHCII complex. The LHCI complex was composed of five subunits, (Lhca1-Lhca5) and the LHCII complex was composed of seven subunits (Lhcb1-Lhcb7). The analysis showed that, except for Lhca5 and Lhcb7, DEGs were detected in the rest of the subunits, and the According to the KEGG enrichment results of DEGs, we showed the top 20 pathways with low Q-values. Among them, the pathways related to photosynthesis were "photosynthesis-antenna proteins", "photosynthesis", and "carbon fixation in photosynthetic organisms", and the DEG numbers were 113, 155, and 93, respectively. In addition, we found two pathways related to plant hormones, which were "plant hormone signal transduction" and "BR biosynthesis", and the DEG numbers were 203 and 33, respectively (Supplementary Table S10, Figure 4B).

The Expressions of Photosynthesis-Related Genes Were Inhibited
The KEGG analysis showed that the pathways of "photosynthesis-antenna proteins", "photosynthesis", and "carbon fixation in photosynthetic organisms" were enriched significantly, and the DEG numbers were all large. Therefore, we analyzed the DEG expressions in the above pathways.
Except for PsbY, all DEGs in the rest of the subunits were downregulated in the BL treatment (Supplementary Table S12, Figure 6A,B). There were 76 DEGs in the "Calvin cycle", among them, RBCS, GAPA, TKT, SBPASE, PRK, and RPE were downregulated in the BL treatment, while RPIA was upregulated (Supplementary Table S13, Supplementary Figure S3). In addition, almost all of the DEGs in "chlorophyll biosynthesis" as well as most of the DEGs in "carotenoid biosynthesis" were downregulated in the BL treatment (Supplementary Tables S14 and S15, Supplementary Figures S4 and S5).

The Expressions of BR Synthase Genes Were Negatively Regulated by BL
According to the results of the KEGG enrichment, the pathways of "plant hormone signal transduction" and "BR biosynthesis" were significantly enriched. Therefore, the expressions of encoding genes in the two pathways were analyzed. The results showed that in "BR signal transduction", multiple members were detected, except BSU1. It is noteworthy that the encoding genes of co-receptor BAK1, transcription factor BZR1/2, and response factor TCH4 were upregulated by BL (Supplementary Table S16, Supplementary Figure S6). Upregulated expression of TCH4 indicated that BR signaling was activated continuously [35]. In addition, we detected many BR synthase genes CYP90B1, CYP724B1, CPD, CYP90C1, CYP90D1, Det2, D2, BR6ox1, and CYP92A6 in A. flabellulatum gametophytes, except BR6ox2. Among them, CYP90B1, CYP90C1, CYP90D1, CPD, and BR6ox1 were downregulated in the BL treatment (Supplementary Table S17, Figure 7A,B). So, the five kinds of BR synthase genes were negatively regulated by BL in A. flabellulatum gametophytes.

Discussion
In angiosperms, sprinkling with 10 −8 M BL can promote the net photosynthetic rate of S. lycopersicum leaves [36]. After pretreating seeds of Pisum sativum with 10 −7 M BL, the chlorophyll content of its seedlings significantly increased [37]. However, for the bes1-D seedlings of A. thaliana (BR signaling was constitutively activated), the chlorophyll content significantly decreased [38]. In spore plants, the chlorophyll and carotenoid contents significantly increased by applying 10 −8 M BL to Chlorella vulgaris and 10 −10 -10 −5 M BL to the green alga Acutodesmus obliquus [39,40]. In the green alga A. obliquus, the promotional role of 10 −6 M BL was the greatest [40]. In gene expression, the expressions of RBCS and RBCL in the "Calvin cycle" were upregulated when 10 −7 M BL was sprayed on Cucumis sativus leaves [41], while the expressions of GLK1 and GLK2 (transcription factors responsible for chloroplast development and photosynthetic gene activation), as well as Lhcb (such as Lhcb1. 4 and Lhcb2.2) were downregulated by BL and/or in bes1-D seedlings of A. thaliana [38]. In this study, almost all of the DEGs in "photosynthesis-antenna proteins", "photosynthesis", and "chlorophyll biosynthesis", as well as most of the DEGs in the "Calvin cycle", and "carotenoid biosynthesis" were downregulated after 10 −6 M BL was applied to A. flabellulatum gametophytes. In conclusion, BL is involved in regulating photosynthesis and its related gene expressions in angiosperms and spore plants (including fern gametophytes).
Whether BL exists in plants will vary according to the plant species. Nowadays, CS rather than BL is detected in liverwort (Marchantia polymorpha), moss (Physcomitrella pat-

Discussion
In angiosperms, sprinkling with 10 −8 M BL can promote the net photosynthetic rate of S. lycopersicum leaves [36]. After pretreating seeds of Pisum sativum with 10 −7 M BL, the chlorophyll content of its seedlings significantly increased [37]. However, for the bes1-D seedlings of A. thaliana (BR signaling was constitutively activated), the chlorophyll content significantly decreased [38]. In spore plants, the chlorophyll and carotenoid contents significantly increased by applying 10 −8 M BL to Chlorella vulgaris and 10 −10 -10 −5 M BL to the green alga Acutodesmus obliquus [39,40]. In the green alga A. obliquus, the promotional role of 10 −6 M BL was the greatest [40]. In gene expression, the expressions of RBCS and RBCL in the "Calvin cycle" were upregulated when 10 −7 M BL was sprayed on Cucumis sativus leaves [41], while the expressions of GLK1 and GLK2 (transcription factors responsible for chloroplast development and photosynthetic gene activation), as well as Lhcb (such as Lhcb1.4 and Lhcb2.2) were downregulated by BL and/or in bes1-D seedlings of A. thaliana [38]. In this study, almost all of the DEGs in "photosynthesis-antenna proteins", "photosynthesis", and "chlorophyll biosynthesis", as well as most of the DEGs in the "Calvin cycle", and "carotenoid biosynthesis" were downregulated after 10 −6 M BL was applied to A. flabellulatum gametophytes. In conclusion, BL is involved in regulating photosynthesis and its related gene expressions in angiosperms and spore plants (including fern gametophytes).
Whether BL exists in plants will vary according to the plant species. Nowadays, CS rather than BL is detected in liverwort (Marchantia polymorpha), moss (Physcomitrella patens), Lycopodiopsida (lycophytes) (Selaginella moellendorffii and Selaginella uncinata), and 12 species of Polypodiopsida (true ferns), including P. aquilinum, O. japonica, Lygodium japonicum, etc. [12]. In plants, the existence of BR6ox2 determines whether CS can be further converted into BL [15][16][17]. It is speculated that BR6ox2 is not present in ferns [12]. The results of our study support this conjecture because BR6ox2 was not detected in A. flabellulatum gametophytes. Although it is indisputable that CS is present in spore plants, the content is generally lower than that in angiosperms. For example, the CS content in spore plants (maximum is 38 pg g −1 fresh weight) is lower than 29% of that in A. thaliana (134 pg g −1 fresh weight), 26% in O. sativa (151 pg g −1 fresh weight), and 1.5% in Z. mays (2580 pg g −1 fresh weight) [12]. This may be related to the low expressions of BR synthase genes (Supplementary Table S17). Additionally, our studies showed that the expressions of five kinds of BR synthase genes (CYP90B1, CYP90C1, CYP90D1, CPD, and BR6ox1) were negatively regulated by BL, which indicated that ferns (including gametophytes) also needed the regulatory mechanism for maintaining BRs homeostasis.

Conclusions
In this study, the genes of fern gametophytes responsive to BRs were identified for the first time by transcriptome sequencing. We found that the expressions of photosynthetic genes were widely inhibited by high concentrations of BL in A. flabellulatum gametophytes. Moreover, many BR synthase genes were detected, except BR6ox2, which may be the reason why CS rather than BL is present in ferns. Finally, the expressions of BR synthase genes were negatively regulated by BL in A. flabellulatum gametophytes.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/genes13061061/s1, Figure S1: A. flabellulatum gametophyte growth response to different BL concentrations; Figure S2: Sequencing saturation analysis and cluster analysis of 6 samples; Figure S3: Most of the DEGs in Calvin cycle were down-regulated in BL treatment; Figure S4: Almost all of the DEGs in chlorophyll biosynthesis were down-regulated in BL treatment; Figure S5: The DEGs down-regulated in carotenoid biosynthesis were in the majority in BL treatment; Figure S6: Expressions of the DEGs in BR signaling; Table S1: Information of primers used in qRT-PCR; Table S2: Summary of transcriptome data for the 6 samples; Table S3: The expression profiles of 210,415 isoforms in the 6 samples; Table S4: 5081 and 3313 genes were up-regulated and down regulated in BL treatment, respectively; Table S5: GO annotation of 8394 DEGs; Table S6: GO functional classification of 5275 annotated DEGs; Table S7: GO enrichment of 5275 annotated DEGs; Table S8: KEGG annotation of 8394 DEGs; Table S9: KEGG pathway classification of 3146 annotated DEGs; Table S10: KEGG enrichment of 3146 annotated DEGs; Table S11: The expression levels of DEGs involved in photosynthesis -antenna proteins pathway; Table S12: The expression levels of DEGs involved in photosynthesis pathway; Table S13: The expression levels of DEGs involved in Calvin cycle; Table S14: The expression levels of DEGs involved in chlorophyll biosynthesis; Table S15: The expression levels of DEGs involved in Carotenoids biosynthesis; Table S16: The expression levels of DEGs involved in BR signaling; Table S17: The expression levels of DEGs involved in BR biosynthesis; Table S18: Processing and analysis of qRT-PCR data.