Genome-wide identification and expression pattern analysis of lipoxygenase gene family in turnip (Brassica rapa L. subsp. rapa)

Identification of the turnip LOX genes

The genome, CDS, and protein sequence files of turnip were downloaded from the Brassicaceae database (http://brassicadb.cn/#/) (Chen et al., 2021). The Hidden Markov Model file for lipoxygenase (PF00305) was downloaded from the Pfam database (https://pfam.xfam.org/) (Mistry et al., 2020), and the turnip genome was searched using HMMER 3.0 software and candidate LOX proteins with an e value < 1.0E−05 were screened (Wheeler & Eddy, 2013) and submitted to the conserved domains database (CDD, https://www.ncbi.nlm.nih.gov/cdd/?term=) to confirm the presence of lipoxygenase and the presence of PLAT/LH2 structural domain (Lu et al., 2020). After excluding sequences without the lipoxygenase and PLAT/LH2 structural domains, the non-redundant sequences were named according to their chromosomal position order. ExPASy (https://web.expasy.org/protparam/) was used to analyze the basic physicochemical properties of LOX proteins. The chloroplast transit peptides and the subcellular localization were predicted using TargetP 2.0 Server (https://services.healthtech.dtu.dk/service.php?TargetP-2.0) and WoLF PSORT (https://psort.hgc.jp/), respectively. Global sequence alignment was performed using the EMBOSS Needle tool (https://www.ebi.ac.uk/Tools/psa/emboss_needle/) to determine the similarity and identity between the LOX members. The conserved motifs of LOX (the number of motifs is set to 20) were analyzed using the MEME suite (https://meme-suite.org/meme/tools/meme) (Bailey et al., 2015), and the gene structure was mapped using TBtools software (Chen et al., 2020).

Phylogenetic analysis

The published LOX protein sequences of Arabidopsis (Umate, 2011), radish (Wang et al., 2019), watermelon (Liu et al., 2020), cotton (Shaban et al., 2018), sorghum (Shrestha, Pant & Huang, 2021), buckwheat (Ji et al., 2020), banana (Liu et al., 2021), and poplar (Chen et al., 2015) were retrieved from (https://www.arabidopsis.org/), BRAD (http://brassicadb.cn/#/), CuGenDB (http://cucurbitgenomics. org/organism/1), CottonGen (https://www.cottongen.org/), NCBI (https://www.ncbi.nlm.nih.gov/), MBKBASE (http://www.mbkbase.org/Pinku1/), Banana Genome Hub (https://banana-genome-hub.southgreen.fr), and Phytozome (http://www.Phytozome.net/), respectively. Multiple sequence alignment was performed using MUSCLE software (Edgar, 2021) and phylogenetic trees were constructed by the neighbor-joining method using MEGA X software (Partial deletion; 95% site coverage cutoff, and 1,000 bootstrap replicates) (Kumar et al., 2018) and visualized using iTOL (https://itol.embl.de/) (Letunic & Bork, 2021).

Chromosome localization and gene duplication analysis

Whole-genome duplication and tandem repeat relationships between BrrLOX genes and genetic covariance between species were detected using MCScanX software (https://github.com/wyp1125/MCScanX) (Wang et al., 2012). Chromosome location was obtained from the gene feature format (GFF) gene annotation files and Visualized gene co-linear relationships using Circos software (http://circos.ca) (Krzywinski et al., 2009). Then, the synonymous substitution rate (Ka) and the non-synonymous substitution rate (Ks) were calculated for gene pairs using KaKs_Calculator 2.0 software (Wang et al., 2010).

Analysis of cis-acting elements in the promoters of turnip LOXs

A 2,000 bp sequence upstream of the start codon of each turnip LOX gene was extracted from the turnip genome database, and the cis-acting elements of the promoter were predicted using PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002).

Expression analysis of RNA-seq data from turnip fleshy roots

The raw transcriptome data of the three developmental periods (each with two independent biological replicates) of turnip fleshy roots were downloaded from NCBI (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA273340/) (Li et al., 2015), the transcriptome data were remapped to the turnip reference genome using HISAT2 (Pertea et al., 2016), and the TPM values (Transcripts Per Kilobase of exon model per Million mapped reads) of turnip LOX were calculated using the featureCounts tool (Liao, Smyth & Shi, 2014), and finally the TPM values were normalized by a logarithmic scale with base 2 and visualized in a heat map using the TBtools software (Chen et al., 2020).

Plant materials

We performed a transcriptome analysis of ‘Qiamagu’ turnip seedlings showing the transcriptome changes caused by four treatments (mock control; 100 µM MeJA spray; 100 mM NaCl watering; 100 µM MeJA spray + 100 mM NaCl watering) on leaves of four-leaf stage seedlings (Table S1), with three biological replicates of each treatment. The plant material used for the stress and hormone treatments was the Xinjiang cultivar ‘Qiamagu’. Turnips were grown at 25 °C in a light incubator with 16 h light/8 h dark photoperiod and 60% relative humidity. Seedlings of uniform size and two weeks old were selected for stress and hormone treatments. To examine the effect of methyl jasmonate (MeJA), salicylic acid (SA), and abscisic acid (ABA) treatments on turnip LOX, leaves were sampled at 0, 4, 8, 12, and 24 h after spraying 100 µM of MeJA, SA, and ABA, respectively. To investigate the effect of abiotic stress on turnip LOX, turnip plants at the four-leaf stage were watered with 100 mM NaCl and 20% PEG6000 to simulate salt stress and drought stress, respectively, and exposed to 4 °C and 40 °C to simulate low and high temperatures, respectively. The leaves were sampled at 0, 4, 8, 12, and 24 h after treatment, immediately frozen in liquid nitrogen, and stored at −80 °C for further use.

Real-time quantitative PCR

Total RNA was extracted from the leaf samples using RNAprep pure Plant Kit (Tiangen, Beijing, China) and reverse-transcribed to generate the first-strand cDNA using RT reagent Kit with gDNA Eraser (Takara, Beijing, China). The cDNA was diluted and used for qRT-PCR. Primers for BrrLOX family genes used in qRT-PCR were designed using Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi) with turnip β-Actin as the internal reference gene (Ye et al., 2012). All primers used in this study are listed in Table S2. Relative expression levels of the BrrLOX genes were calculated following the 2^−ΔΔCt method. Three biological replicates were set up for each time point per treatment, with each repetition consisting of nine plants.

Weighted gene co-expression network analysis (WGCNA)

Functional annotation of turnip CDS sequences using EggNOG-mapper (version 2.1, http://eggnog5.embl.de/#/app/emapper) (Huerta-Cepas et al., 2017) and WGCNA(version 1.69) analysis of successfully annotated genes in transcripts (Zhang & Horvath, 2005) (Mock-control; 100 µM MeJA spray; 100 mM NaCl watering; 100 µM MeJA spray + 100 mM NaCl watering). In the analysis, the parameters were set as follows. an optimal soft threshold (power) value of 7; minModuleSize of 30, and mergeCutHeight of 0.25. Finally, all co-expression networks associated with BrrLOXs were extracted, and the edges with weights below 0.2 were filtered out. We visualized the network connections using Cytoscape software (version 3.8.0, https://cytoscape.org/) (Shannon et al., 2003). Further, the functional annotation of BrrLOXs and co-expressed genes was performed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) using ggplot2 to identify the top 10 GO annotations and the top 20 KEGG functional annotations (Kanehisa et al., 2016).

Identification of turnip LOX gene family members

A total of 15 BrrLOX candidate genes were identified in the turnip genome (Table 1, Table S3) and named sequentially from BrrLOX1 to BrrLOX15 based on their chromosomal location. The distribution of BrrLOX genes on the chromosomes is shown in Fig. 1. The BrrLOX genes were distributed on 6 out of the 10 turnip chromosomes, of which A07 contained seven LOX genes, A02 contained three, A01 contained two, and A05, A08, and A09 contained only one each.

The CDS lengths of the 15 predicted BrrLOX family genes ranged from 2,499 to 2760 bp, respectively, and encoded BrrLOXs proteins from 832 to 919 amino acids in length. The predicted molecular weight of BrrLOX proteins ranged from 94.97 to 104.76 kDa, with a theoretical isoelectric point range of 5.18 to 7.9. The chloroplast transit peptide is present in 11 members of BrrLOXs. Subcellular localisation predicted by WoLF PSORT shows that nine of the BrrLOXs proteins are localised in the chloroplast and the remaining six are localised in the cytoplasm. The protein sequence alignment revealed a sequence similarity of 51.30%–100% and a sequence identity of 35.70%–100% among the BrrLOX members (Fig. 2). Specifically, BrrLOX7 and BrrLOX8 showed the highest similarity and identity (100%). The lowest similarity (51.3%) was between BrrLOX4 and BrrLOX6, BrrLOX6 and BrrLOX13; BrrLOX4 and BrrLOX6 showed the lowest identity (35.70%).

Phylogenetic relationship of LOX family members in different species

We constructed phylogenetic trees based on 128 LOX protein sequences from 9 different species, including turnip (Brassica rapa L. subsp. Rapa), radish (Raphanus sativus L.), Arabidopsis (Arabidopsis thaliana), cotton (Gossypium hirsutum L.), sorghum (Sorghum bicolor L.Moench), banana (Musa acuminata), buckwheat (Fagopyrum tataricum L.), poplar (Populus trichocarpa) and watermelon (Citrullus. lanatus subsp. vulgaris), using the neighbor-joining algorithm in MEGA X software (Fig. 3 and Table S4) (Kumar et al., 2018). The LOX protein sequences were classified into 9-LOX, 13-LOX type I, and 13-LOX type II. Of the 15 BrrLOXs identified, BrrLOX1, BrrLOX2, and BrrLOX6 belong to the 9-LOX subfamily, while the other 12 to 13-LOX type II. The BrrLOX gene family analyzed lacked the 13-LOX type I, a result also found in Arabidopsis thaliana and other Brassica species (Brassica rapa, Brassica oleracea and Raphanus sativus) (Wang et al., 2019).

Conserved motifs and gene structure of turnip LOXs

The protein motifs of the turnip BrrLOX proteins were analyzed using MEME suite. The analysis revealed 16 conserved motifs (motif2–15, motif18, motif19) widely present in all BrrLOX proteins (Supplementary Figure_S1). The BrrLOXs showed some degree of specificity in motif distribution, with all LOXs except BrrLOX5 containing motif17 (Fig. 4B). Meanwhile, all LOXs except BrrLOX10 have motif1, and motif1 contains a histidine-rich (His) 38 amino acid residue motif consisting of His-(X)4-His-(X)4-His-(X)17-His-(X)8-His. In addition, BrrLOX13 has an additional motif15 at the N-terminus, and BrrLOX14 has an additional motif8 at the N-terminus. These differences in motif arrangement may be responsible for the functional differentiation of the turnip LOX proteins. Our analysis of the gene structure showed that the BrrLOX genes consist of 5–8 introns, with BrrLOX5 having the fewest introns (five) and BrrLOX2, BrrLOX4, BrrLOX6, and BrrLOX13 having the most (eight) (Fig. 4C).

Gene duplication and collinearity analysis of the BrrLOX genes

In the evolution of plants, gene duplication regulates the acquisition of new genes and the creation of new genetic traits (Zhang, 2003). The present study explored the gene duplication events of the LOX family of turnips using the MCScanX software (Wang et al., 2012). In the BrrLOX gene family, seven tandem replication gene pairs (BrrLOX7/BrrLOX8, BrrLOX9/BrrLOX10/BrrLOX11/BrrLOX12) and four whole gene replication (WGD) or fragment replication (SD) gene pairs (BrrLOX2/BrrLOX6, BrrLOX4/BrrLOX13, BrrLOX5/BrrLOX14, BrrLOX14/BrrLOX15) were identified. Furthermore, Ka and Ks were calculated to estimate the evolutionary selection pressure on the duplicate gene pairs of turnip LOX genes.

The Ka/Ks values for all 11 LOX duplicate gene pairs were less than one (Table 2), indicating that LOX family genes are mainly influenced by purifying selection during evolution. Analysis of the collinearity of LOX genes in Arabidopsis, radish, and turnip by interspecies collinearity (Fig. 5) identified that 8 members of the 15 BrrLOX genes were collinear with 5 AtLOX genes, giving rise to a total of 11 interspecies collinear gene pairs, mainly on chromosome 2 of turnip. Ten members were collinear with nine RsLOX genes (Fig. 5; Table S5), yielding a total of 17 interspecies collinear gene pairs, mainly on chromosome 2. These collinear gene pairs provide a basis for functional studies on the BrrLOX genes.

Analysis of cis-elements in BrrLOX gene promoters

We retrieved a 2,000 bp sequence upstream of the start codon of the BrrLOX genes and performed a cis-regulatory element prediction analysis. The analysis identified 38 cis-acting elements classified into four main categories: phytohormone responsive elements, growth and metabolic responsive elements, stress responsive elements, and light responsive elements (Fig. 6). The main phytohormone response elements included TGA-element associated with auxin, GARE-motif, P-box, and TATC-box with gibberellin, ABRE with abscisic acid, CGTCA-motif and TGACG-motif with methyl jasmonate, and TCA-element with salicylic acid, suggesting multiple phytohormones regulate the BrrLOX genes. The primary growth and metabolism response elements included circadian, GCN4-motif, Ry-element, and O2-site involved in circadian control, endosperm expression, seed-specific regulation, and zein metabolism regulation, suggesting the roles of BrrLOX genes in plant growth and metabolism. The major stress response elements included TC-rich repeats, LTR, ARE, and GC-motif related to defense and stress, low-temperature response, anaerobic induction, and anoxic-specific inducibility. Light response elements, mainly Box 4, G-box, and GT1-motif, exist in all BrrLOX promoters regions.

Weighted gene co-expression network analysis (WGCNA) of BrrLOXs

There were 42,877 functional genes in the turnip CDS sequence that were successfully annotated by EggNOG-mapper (version 2.1) and a co-expression network was constructed using WGCNA based on these functional genes and four transcriptomic data. Co-expression networks associated with BrrLOX genes were visualized using Cytoscape software (version 3.8.0), and a total of seven co-expression modules, named modules A to G, were obtained (Fig. 7). The largest of these modules is module B, with 5,426 genes co-expressed with BrrLOX2 and BrrLOX5. The smallest is module D with 383 genes co-expressed with BrrLOX6. Furthermore, GO enrichment analysis showed that genes co-expressed with BrrLOXs were mainly enriched in the secondary metabolic process, cellular response to lipid, response to a bacterium, sulfur compound metabolic process, and carbohydrate derivative metabolic process in the biological process category (Fig. 8A). Gene annotations were mainly enriched in the cellular components in the thylakoid, extracellular region, plastid thylakoid, chloroplast thylakoid, and plastid membrane. From molecular functions, genes were primarily enriched in tryptophan metabolism, phenylpropanoid biosynthesis, carbon fixation in photosynthetic organisms, MAPK signaling pathway-plant, and 2-oxocarboxylic acid metabolism. According to KEGG functional enrichment analysis, these genes co-expressed with BrrLOXs were mainly enriched in tryptophan metabolism, phenylalanine biosynthesis, carbon fixation in photosynthetic organisms, MAPK signaling pathway-plant, and 2-oxocarboxylic acid metabolism (Fig. 8B).

Transcriptome analysis of the LOX family genes of turnip

To investigate the expression of turnip LOX family genes during fleshy roots development, we reanalyzed the expression levels of BrrLOX genes using publicly available RNA sequence data from three fleshy roots development stages, including ES (early stage before cortex splitting), CSS (cortex-splitting stage) and RTS (secondary root-thickening stage) (Fig. 9; Table S6), in which BrrLOX1 and BrrLOX13 were highly expressed at all three stages; BrrLOX11 expression was higher in ES and RTS than in CSS, while BrrLOX9 and BrrLOX10 expression levels were higher in RTS than in ES and CSS. BrrLOX12 was highly expressed in RTS. However, BrrLOX3 and BrrLOX15 were expressed at 3 periods with little or no expression.

Expression analysis of BrrLOX family genes in turnip under abiotic stress and exogenous hormone treatment

We further performed qRT-PCR analysis to detect the changes in the transcript levels of LOX genes in turnip leaves under different abiotic stresses (salt stress, drought stress, cold stress, and heat stress) (Fig. 10A and Table S7) and exogenous hormone treatment (MeJA, SA, and ABA) (Fig. 10B and Table S7). Under NaCl treatment, the expression of most genes (BrrLOX1-6, BrrLOX9, BrrLOX11-15) was downregulated, while that of BrrLOX7 and BrrLOX8 was upregulated at 4 h and 8 h after exposure; BrrLOX10 was upregulated at 12 h. Under drought stress, 10 LOX genes were upregulated, among which BrrLOX9 was most significantly expressed at 12 h. Under low-temperature treatment, BrrLOX3-5, BrrLOX11, BrrLOX14, and BrrLOX15 genes were upregulated at all time points, with BrrLOX15 showing a significant high expression at 8 h; BrrLOX13 was downregulated at all time points. Under high-temperature treatment, BrrLOX2, BrrLOX5, and BrrLOX11 were upregulated at all time points, with the most significant expression for BrrLOX15 at 4 h.

Furthermore, the BrrLOX genes showed different expression patterns under various hormonal treatments (Fig. 10B). BrrLOX1 gene expression was downregulated under the MeJA treatment and upregulated under the ABA treatment. BrrLOX4, BrrLOX7-12, and BrrLOX14 were upregulated with three exogenous hormones, with the highest effect for BrrLOX9 at 24 h.