Phylogenomics-Based Reconstruction and Molecular Evolutionary Histories of Brassica Photoreceptor Gene Families

Photosensory proteins known as photoreceptors (PHRs) are crucial for delineating light environments in synchronization with other environmental cues and regulating their physiological variables in plants. However, this has not been well studied in the Brassica genus, which includes several important agricultural and horticultural crops. Herein, we identified five major PHR gene families—phytochrome (PHY), cryptochrome (CRY), phototropin (PHOT), F-box containing flavin binding proteins (ZTL/FKF1/LKP2), and UV RESISTANCE LOCUS 8 (UVR8)—genomic scales and classified them into subfamilies based on their phylogenetic clustering with Arabidopsis homologues. The molecular evolution characteristics of Brassica PHR members indicated indirect expansion and lost one to six gene copies at subfamily levels. The segmental duplication was possibly the driving force of the evolution and amplification of Brassica PHRs. Gene replication retention and gene loss events of CRY, PHY, and PHOT members found in diploid progenitors were highly conserved in their tetraploid hybrids. However, hybridization events were attributed to quantitative changes in UVR8 and ZTL/FKF1/LKP2 members. All PHR members underwent purifying selection. In addition, the transcript expression profiles of PHR genes in different tissue and in response to exogenous ABA, and abiotic stress conditions suggested their multiple biological significance. This study is helpful in understanding the molecular evolution characteristics of Brassica PHRs and lays the foundation for their functional characterization.


Introduction
Plants have a series of unique photoreceptors (PHRs) with different wavelength absorption spectra and biochemical properties to delineate light environments precisely and regulate their physiological variables [1]. In the model plant Arabidopsis thaliana, five classes of PHRs, namely, phytochrome (PHY), cryptochrome (CRY), phototropin (PHOT), F-box containing flavin binding proteins (ZEITLUPE (ZTL)/LOV KELCH PROTEIN1 (LKP1), FLAVIN-BINDING KELCH REPEAT F-BOX 1 (FKF1), and LOV KELCH PROTEIN2 (LKP2)) and UV RESISTANCE LOCUS 8 (UVR8), control photoreception and corresponding light signaling networks. Plants optimize the light absorption through these PHRs, in synchronization with other environmental cues, as it is essential to balance photochemistry and photosynthesis reaction rates for the adaptive responses [2]. Of these, PHYs are useful in perceiving red/far-red lights (600-750 nm), while CRY, PHOT, and ZTL/FKF1/LKP2 are for detecting blue/UV-A light (320-500 nm) [3]. Similarly, UVR8 is a cytosol and nuclear-localized photoreceptor responsible for UV-B absorption and protecting plants from harmful effects of UV-B light (280-320 nm) [4]. The UVR8-mediated protective effect could partially come from induced flavonoid biosynthesis [5]. Besides enabling a photosynthetic lifestyle, these light-activated PHRs help plants guide various critical growth and developmental processes, including seed germination, flowering time, and the circadian Therefore, to get insight into the evolution and expansion of five main PHR gene families-PHYs, CRYs, PHOTs, UVR8, and ZTL/FKF1/LKP2-were comprehensively and systematically analyzed in three diploids and two tetraploids of Brassica species. BLASTP and domain architecture of Arabidopsis homologues were considered for Brassica PHR gene family identification. The phylogenetic clustering of Brassica PHR gene members was used for family and subfamily classifications. Further, the gene structure, gene synteny between PHR members in a single genome, gene duplicates including segmental and tandem, chromosomal distribution pattern, and gene losses after WGT events, were analyzed to understand the evolutionary past and/or expansion histories of five major PHR gene families in Brassica. In addition, the transcript expression profiles of PHR genes across different tissue types (root, leaf, bud, silique, and callus) and in response to exogenous ABA application, abiotic (light, drought, salt, and cold), and biotic stress conditions (Pectobacterium carotovorum infection) suggested PHR genes play important roles in plant growth and development and stress responses. The presence of ABA-, stress-associated putative promoter motifs along with light response cis-elements further validate their multisensory functionality. The present study unearthed several new members of PHR genes, and the molecular evolutionary characteristics of PHR will be useful in further studies focusing on functional characterization of photoreceptor genes in Brassica species.

History of Replication Retention and Gene Loss Events in Photoreceptor Genes in Brassica Species
We performed gene loss and replication retention analysis to elucidate the evolution of the photoreceptor genes in Brassica. We obtained quantitative changes in the number of genes in subfamilies based on the phylogenetic reconstruction (Table 1). For diploid Brassica, one A. thaliana gene should theoretically correspond to three genes, but there are different amounts of gene loss (fractionation) ranging from zero to three. The whole-genome triplication (WGT) event in the Brassica genus contains three diploids B. rapa (AA, n = 10), B. oleracea (CC, n = 9), B. nigra (BB, n = 8), and allopolyploids like B. napus (AACC, n = 19) and B. juncea (AABB, n = 18) theoretically should result in more gene copies than that of A. thaliana orthologues. However, the systematic identification and classification of Brassica PHR genes identified 0-2 gene copies in diploids and 0-4 in tetraploids. For instance, CRY1, 2, and 3 subfamilies had one copy in diploids and 2 in their tetraploids. However, the new CRY3-like subfamily had a maximum of 2 and 4 in di-and tetraploid species. Interestingly, the CRY3 subfamily was lost in B. rapa, or it might have evolved as a CRY3like subfamily. In case of the PHY gene family, PHYD is not present in both A and B genomes of di-and tetraploid species, while all other subfamily members, including new PHYA-like, were observed (Table 1).   -1  1  1  2  3   Total  3  3  3  3  6  7 Numbers represent the number of Arabidopsis homologues observed in each species at the subfamily level. The "L" indicates gene loss and the number after "L" represents the number of gene copies lost after the triploidization event.
In total, 5-6 PHYs in diploids and 10-11 PHYs in tetraploids were identified. Among ZTL/FKF1/LKP2 members, only LKP2 and ZTL/FKF1/LKP2-like members were observed in Brassica species. Another family, PHOT, comprised of three subfamilies, namely PHOT1, 2, and PHOT2-like, which were present in all species. All the parental diploids have one member, PHOT1, like AtPHOT1. Unlike PHOT1, PHOT2 has retained duplicates in B. rapa and B. oleracea, albeit B. nigra remained with a single copy. The UVR8 phylogenetic classification added three new subfamilies (UVR8-like 1, 2, and 3) to this gene family. Of these, the members of UVR8 like-1 in B. oleracea, UVR8-like 3 in B. nigra and in B. juncea were not found in this study. The key observation is that the diploid B. nigra and tetraploid B. juncea unusually share the same number of gene copies for UVR8. B. napus had the largest UVR8 family consisted of eight members.
Whole-genome duplication (WGD), segmental/tandem duplication, inversions, and translocations are driving forces of genome evolution that increase genome plasticity. Based on the gene copy number and their distribution pattern in a single genome, they can be classified as singletons, dispersed duplicates, proximal duplicates, tandem duplicates, and segmental/WGD duplicates [24]. According to [24], we determined the evolutionary events/gene types in Brassica PHR gene families. In this study, 9 gene pairs each from B. rapa and B. oleracea, 4 pairs of B. nigra, 46 pairs of B. napus, and 12 pairs of B. juncea were the results of segmental duplication, thus altogether accounting for 94.73% of the total duplication events in PHR gene families (Table S3; Figure 6). On the other hand, one pair of tandem duplication each from B. rapa, B. oleracea, and B. juncea, and two tandem duplications from B. napus were observed. It is worth noting that tandem duplication was observed only in ZTL/FKF1/LKP2 gene families in Brassica, with the exception of B. nigra, where no tandem duplicates was observed. The tandem duplicates share 78 to 92 AA similarities and are adjacent paralogues on the same/single chromosomes. The calculated Ka/Ks ratio for the duplicates ranging from 0.03 to 0.40 shows that the aforementioned genes were under strong purifying selection. The lowest Ka/Ks ratio (0.03) was observed for BjuPHOT2-like 2: BjuPHOT2-like 1, while another gene pair (BnUVR8 like 2: BnUVR8 As demonstrated by [25], we also investigated the effect of hybridization events by comparing the number of PHR genes lost between parental diploids and their hybrids/tetraploids in this study. Gene losses or divergence in CRY, PHY, and PHOT of B. napus and B. juncea were directly reflected in parental diploids. For example, two copies were lost in the CRY1 subfamily of B. rapa and B. oleracea, which accounted for 4 lost copies in tetraploid B. napus. This trend was also true for PHY and PHOT members. However, the expansion or gene losses in UVR8 and ZTL/FKF1/LPK2-like subfamily of B. napus and B. juncea are not originally from their diploid parents and possibly the effects of hybridization events.

Domain Architecture and Physiochemical Properties of PHRs
The Brassica LKP2 protein might have lost one or two Kelch_1 domain at the third and fourth positions at C-terminal ends during evolution, while PAS_9 and F-box-like domains were present as in Arabidopsis LKP2 homologues (Table S4, Figure 5B CRY2-like (~575AA) were shorter in AA length than CRY2 subfamily members (~623AA); hence the molecular weight and the hydrophilic coefficient of CRY2-like were smaller than CRY2 while the pI ranges between 5.56 to 5.99. Both CRY2 and CRY2-like were localized to the chloroplast (Table S5). Unlike AtCRY3, which was localized to chloroplast and mitochondria, Brassica CRY3 members were predicted to localize only in the chloroplast, while all the CRY3-like members were likely to present in chloroplast and mitochondria. In general the Brassica CRY family has 40 members distributed across 5 species, BnCRY3-like 4 being the largest in terms of length (813 AA), while the shortest one was also from the same subfamily, i.e., BnCRY3-like 1 (516 AA). All Brassica PHY members, including newly reported PHYA-like members, were likely destined to work in the nucleus with~1131 AA and pI ranging from 5.55 to 6.29. similarly, ZTL/FKF1/LKP2 members localized to the nucleus while the pI of ZTL/FKF1/LKP2-like members was higher than others. The PHOT members are likely to be present in cell membranes at cellular levels. Of these three members, BraPHOT1, BjuPHOT1-1, and BjuPHOT2-1 were likely to be present in the nucleus in addition to cell membranes. The mean AA length is~911 AA (lowest-727, highest-996), and the pI varies between 6.19 to 8.7. Unlike other PHR genes, the UVR8 members were found with single to multiple cellular localization signals, most of which have nucleus localization signals followed by cell wall and cytoplasm. Distinctly, BolUVR8-like 3 was predicted to present in the Golgi apparatus in addition to the cytoplasm and nucleus.

Gene Structure and Promoter Motifs Analysis
The structural comparison of two or more gene of a single family will aid in understanding the evolutionarily structural changes/gene diversification that occurred in those genes. The distribution patterns of the 45 of 143 Brassica PHR exon/introns (e.g., number of exons, order and their types (symmetrical and asymmetrical)) were similar to that of Arabidopsis homologues, indicating these genes were evolutionarily conserved. This includes three UVR8 genes (BolUVR8, BnUVR8-1, and BniUVR8), all genes of the ZTL/FKF1/LKP2 family except BraLKP2-1, nine members of PHY (BniPHYA, BnPHYA-1, BolPHYA, BnPHYB-1, BolPHYD, BnPHYD, BniPHYA, BnPHYA-1, and BolPHYA), four members of PHOT (BolPHOT1, BniPHOT1, BolPHOT2-1, and BniPHOT2), and eight of CRY (BnCRY2-2, BolCRY2, BnCRY2-1, BniCRY2, BolCRY2-like, BnCRY2-like 1, BnCRY2-like 2, and BnCRY1-1). At the levels of gene family, the intron/exon organization of PHY, ZTL/FKF1/LKP2 were fairly conserved over the members of other PHR gene families. In terms of species, B. napus genes were largely conserved over other species. It is worth noting that the PHR genes of B. rapa were the least conserved among others. This could be majorly contributed to BraLKP2-1, BraUVR8, BraUVR8-like 1, BraUVR8-like 2, and BraUVR8-like 3 having 4, 5, 8, 2, and 5 exons, respectively, against 2,12,12,12,12 of their orthologues. Similarly, the structure of BraPHYA-like looks different from others with 10 exons, including 3 asymmetrical exons. BraPHYC has 6 instead of 3 exons. Brassica PHOTs have 19-22 exons except for BraPHOT1, BraPHOT2-2, and BraPHOT2-like had 4, 1, and 4 exons, respectively. BraCRY3-like 1 has just an exon against 12 of AtCRY3 and Brassica genes belonging to this family, while another BraCRY1 had 6 against 4. Among CRY subfamilies, CRY2 was largely conserved over others. Based on the presence of the number of introns, the PHR genes can be classified into two types: a gene that contains no or one intron, and the gene that contains multiple introns. BraCRY3-like 1 has no introns, while BraUVR8-like 2, BraPHOT2-2, and all of the ZTL/FKF1/LKP2, with the exception of BraLKP2-1, contain just one intron. Most other genes contain 2-22 introns, and the PHOT gene family has the highest number of introns among other PHR gene families. PHOTs had symmetrical exon (0,0 class) at the 3 end with exception of BraPHOT2-2 and BraPHOT2-like. Conversely, except for BraLKP2-1, all of the ZTL/FKF1/LKP2 genes have asymmetrical exons (0,1 or 0,2) at the 3 end. The tandem duplicated gene pairs (Figure 7) have structural similarities between them and are conserved for exon/intron organization and intron phases, with the exception of B. rapa genes.
The promoter motifs/cis-regulatory elements of a gene can determine their regulatory pathways/downstream gene expression in response to perceived internal and external signal cues in plants. Therefore, we used in silico promoter analysis for Brassica PHR genes. The results showed they are enriched with several motifs, including ABA, light, plastid, circadian, growth and development, metabolites, photosynthesis, and stress-related. By abundance, the ZTL/FKF1/LKP2 gene family has the largest number (~96.63 motifs) of light response motifs, followed by PHY, PHOT, CRY, ZTL, and UVR8. Interestingly the ranking of ZTL/FKF1/LKP2 was also true to ABA response motifs. Metabolites, stress, growth, and development-related motifs were relatively higher in PHOT members than in other PHR genes. Similar analysis revealed that PHY has the highest number of motifs associated with plastids and photosynthesis. In comparison with UVR8-like genes, UVR8 genes had a markedly larger mean number of ABA and stress-related motifs. Similarly, the promoters of UVR8, and PHY had higher average for ABA and stress-related motifs than UVR8-like and PHY-like subfamilies. In contrast, CRY-like subfamilies average number of ABA and stress motifs were higher than their classical subfamilies. However, PHOT and PHOT-like subfamilies have almost similar motif compositions.

Expression Profiling of PHR Genes in Different Tissue and Stress Conditions
Due to several ABA-, stress-associated promoter motifs in PHR genes, we investigated the microarray-based expression profiles of 12 B. rapa PHR genes in the presence of ABA hormones or stress conditions like drought, salinity, cold, and P. carotovorum infections at different time intervals. True to the enriched promoter motifs, most of the genes were differentially expressed under ABA and stress conditions. BraUVR8-like 2, BraPHYA, BraPHOT1, BraLKP2, BraUVR8, and BraZTL/FKF1/LKP2-like genes were upregulated by those treatments at one or more intervals. BraZTL/FKF1/LKP2-like was upregulated during abiotic stress conditions, especially under cold stress highest upregulation was observed. BraPHYA was the only gene that showed consistent upregulation under P. carotovorum infection while most other genes were downregulated. Under salinity, BraPHYA, BraPHOT2-1, BraCRY1, BraPHOT1, and BraPHYB were consistently downregulated. By contrast, none of the genes was consistently downregulated under ABA presence.

Discussion
Plants utilize PHRs to monitor almost all facets of light (e.g., intensity, direction, duration, and wavelength of light) to optimize their growth and development. Plant PHRs consisted of specific photoreceptive domains that are useful in perceiving, interpreting, and transducing light signals to nuclei, where transcriptional reprogramming takes place, which ultimately offers developmental plasticity in plants [26][27][28]. Before this study, there was no systematic study done for the identification of PHR genes in Brassica species. Therefore, in the present study, we used systematic approaches at genomic scales to iden- Genomic replication often results in the expansion of the gene family [20,21,25]. Therefore, to understand the quantitative changes that occurred in the PHR gene family following WGT in Brassica, we obtained the quantitative changes in the number of gene copies for each subfamily. Following WGT, every diploid Brassica species theoretically should contain three Arabidopsis homologues and six copies in tetraploid Brassica species. However, none of the Brassica PHR members contains the aforementioned copies in any of the tetraploids or their progenitors, suggesting PHR genes were lost on a large scale, possibly due to functional redundancy between the replicates after the triploidization. The large-scale loss/genome contraction was also noted for other gene families like AHLs, EXPs, and flowering-time genes in Brassica [25,29,30]. In particular, heavy gene losses have occurred for CRY3 and PHYD subfamilies in all the Brassica species. In PHYD, there were no copies found in B. rapa (A genome), B. nigra (B genome) and B. juncea (A and B genome). In contrast, a single copy was noted in both B. oleracea (C genome) and B. napus (A and C genome), suggesting PHYD is possibly required only in species with the C genome. An interesting recent segmental duplication occurred in B. napus (for BnUVR8 like 2: BnUVR8 like 2-1, BnUVR8 like 2: BnUVR8 like 2-2) was 3.7 MYA. Regardless of the species, segmentally duplicated gene pairs of PHY undergo a strong purifying selection than that of segmental/tandem duplicated gene pairs found in ZTL/FKF1/LKP2, CRY, UVR8, and PHOT members. This strong purifying selection would be the reason that PHYA, PHYB, PHYC, and PHYE subfamilies do not retain replicates and maintain single-member subfamilies in diploid genomes. Moreover, in tetraploids, the expression pattern between duplicated gene pairs such as BnPHYB-1: BnPHYB-2, BnPHYE-1: BnPHYE-2, and BnPHYA-like 1: and BnPHYA-like 2 in response to differential light qualities and in different tissue are similar (Figure 7) indicating the possible functional redundancy between them, which also can be another reason for their strong negative selection. Although there were several tandem duplicates found in Brassica species [35], only ZTL/FKF1/LKP2 was found with tandem duplicates among all the PHR gene families in Brassica species in this study.
In conclusion, this is the first comprehensive and systematic study to identify 144 photoreceptor genes in Brassica species. This comprehensive study reported several new subfamilies, including CRY2-like, CRY3-like, PHYA-like, PHOT2-like, UVR8-like 1, 2, 3, and ZTL/FKF1/LKP2-like based on their phylogeny with Arabidopsis homologues. The molecular evolution of PHR showed indirect expansion and gene losses at a large scale. In particular, PHYD was completely lost in the A and B genomes and retained only in the C subgenome. In comparison, the UVR8 was markedly expanded than other PHR gene families in Brassica. Segmental duplication at large and tandem at small scales are responsible for expanding PHR gene families. All of the PHR genes undergo purifying selection, yet relatively stronger selection was observed in PHY genes. The presence of several important cis-elements in PHR promoter motifs and the expression of some PHR genes in the presence of exogenous ABA, differential light, stress conditions, and tissue-specific expression patterns suggest that PHRs are important in multiple biological processes. Collectively, this study provides basic resources for functional characterization of PHR genes, which will be helpful in Brassica crop-improvement programs to harness the desired agronomic traits.

Phylogenetic Analysis and Classification of Photoreceptor Gene Family/Subfamily Members
The multiple BLASTP hits with CDs of each gene family and Arabidopsis queries were aligned using the MUSCLE alignment tool of MEGA11 [44]. The alignment options were kept as default, with UPGMA as the cluster method for each gene family. The unrooted phylogenetic tree was constructed using the neighbor-joining tree method with 1500 bootstrap replications, and Poisson as substitution model. The other parameters were default. The classification was based on the phylogenetic clustering of Brassica sequences with Arabidopsis subfamily members/homologues. A subfamily's single-copy gene gets the name of Arabidopsis homolog present in the same clusters with the exception of the initials of its species of origin. Additionally, a serial number was placed after the gene name as per the phylogenetic ranking if multiple gene copies are present in a subfamily. The phylogenetic outliers were designated as novel/new subfamily, and the gene nomenclature is partly obtained from possibly closest phylogenetic clusters.

Gene Duplication, Synteny, Selection Pressure Analysis
Tbtools v1.09873 was used to investigate gene duplication events and the syntenic relationship between members of a particular Brassica photoreceptor gene family. Briefly, the protein sequences of B. rapa, B. napus, B. oleracea, B. juncea, and B. nigra downloaded from the BRAD database were individually self-BLASTED by the "Blast several Sequences to a Big Database." The parameters set were Outfmt: Table; NumofThreads: 2; E-value: 1 × 10 −5 ; NumofHits: 5 and NumofAligns: 5. Meanwhile, the Gff gene annotation file for each genome was preprocessed/simplified with "File Merge For MCScanX" with Merge mode set to "GtfGff2SimGxf." MCScanX was executed by "Quik Run MCScanX Wrapper" [24] for each genome separately to find collinearity genes and gene types (tandem/segmental/proximal/dispersed) at the genomic scale. The resultant gene collinearity information fed to "File Merge For MCScanX" (merge mode: collinear) to identify gene pairs/duplicated genes of desired gene family. The linked regions from gene pairs and respective merged Gff files were utilized to identify synteny through "File Transformat for MicroSynteny Viewer." The chromosome features of all five species were derived from individual gff annotation files, and the syntenic relationships of selected genes were visualized through "Advanced Circos" of Tbtools. The "simple Ka/Ks Calculator (NG)" was utilized to infer the synonymous (Ks) and nonsynonymous (Ka) nucleotide substitutions in tandem or segmentally duplicated genes and their ratios (Ka/Ks) to measure the selection pressures (Ka/Ks > 1 (positive) or Ka/Ks < 1 (purifying)) during evolution. According to [45], the divergence time (MYA, million years ago) of duplicated genes was calculated with the formula T = Ks/2r, where r is 1.5 × 10 −8 synonymous substitutions per site per year, and it is the rate of divergence for nuclear genes from plants.

Promoter Motif Analysis
To get the putative promoter regions of each PHR gene family of Brassica, the genome sequences and gff annotation files of respective Brassica species downloaded from the BRAD database were used. To extract 2 kb upstream bases of each gene reported in this study, "Gtf/Gff3 Sequence extract" under GTF/GFF3 manipulate command of Tbtools v1.09873 was utilized. GFF3 and genome of Brassica species were separately fed as input; parameters were set to extract 2 kb bases upstream of the transcription start site of each gene and requested to retain only upstream bases. The resultant putative promoter sequences of selected genes were separated, and random manual verification was done to confirm the promoter sequences of desired genes. For cis-acting element analysis, the online tool New PLACE (https://www.dna.affrc.go.jp/PLACE/?action=newplace) (accessed on 9 May 2022) [46] was used. The abundance of putative promoter motifs belonging to key biological processes was represented by heatmap using Tbtools.

Gene Structure Analysis
The gene structure of each PHR family was visualized by the "Gene Structure View (Advanced)" of Tbtools [37]. The tree layout was set to the origin, and the Intron phase number was included to know the distribution of symmetrical and asymmetrical exons for understanding the evolutionarily structural changes.

Expression Profiling of Photoreceptor Genes
The plant material, growth conditions, treatments, and cDNA synthesis for microarray experiments were originally from a previous study [47]. Briefly, the leaves of B. rapa cv. Chiifu inbreed seedlings (3-week-old) grown in soil-containing plastic pots, maintained at standard growth conditions (16 h/8 h photoperiods, 25 • C temperature, and 70% relative humidity), were sprayed with abscisic acid (ABA) (100 µM)/NaCl (250 mM; salt)/ polyethylene glycol 6000 (4% w/v; drought)/Pectobacterium carotovorum suspension (6.15 log 10 colony-forming units/mL) and whole plants were sampled at different intervals. Similarly, for inducing cold stress, plants were kept at 4 • C and sampled. Plants grown at standard growth conditions without chemical/stress agents/mocked samples (e.g., 0.2% Tween 20 and 0.1 DMSO for ABA) were designated as control and collected simultaneously with respective treated samples. Samples in triplicates for each treatment were used for total RNA extraction, cDNA synthesis, and subsequent analysis with microarray. The microarray hybridization reaction and data analysis were performed by GGBIO (http://www.ggbio.com) according to the manufacturer's instructions (NimbleGen Inc., Madison, WI, USA; GenePix scanner 4000B (Axon, Scottsdale, AZ, USA)). The relative expression changes (Log2FC) of PHR transcripts in comparison with controls (0 h) were calculated and presented as heatmap. Additionally, the normalized expression values of PHRs (FPKM values) in Brassica napus var. napus (cv. zhongshuang 11) treated with differential LED light qualities (red (100%), blue (100%), and compound lights (75R:25B or 75B:25R)) were extracted from the Brassica Expression DataBase (BrassicaEDB), v1.0 [48] and log2 [FPKM] values of BnPHRs were used to generate a heatmap with Tbtools. Also, the expression changes of BnPHRs across five different tissue types (root, calli, silique, bud, and leaf) were extracted from the public database (https://bnaomics.ocri-genomics.net/ tools/exp-view/tissue.php) (accessed on 4 June 2022) and presented as a heatmap.

Conflicts of Interest:
The authors declare no conflict of interest.