Genome-wide characterization, expression analyses, and functional prediction of the NPF family in Brassica napus

Background NITRATE TRANSPORTER 1/PEPTIDE TRANSPORTER (NRT1/PTR) family (NPF) members are essential transporters for many substrates in plants, including nitrate, hormones, peptides, and secondary metabolites. Here, we report the global characterization of NPF in the important oil crop Brassica napus, including that for phylogeny, gene/protein structures, duplications, and expression patterns. Results A total of 199 B. napus (BnaNPFs) NPF-coding genes were identified. Phylogenetic analyses categorized these genes into 11 subfamilies, including three new ones. Sequence feature analysis revealed that members of each subfamily contain conserved gene and protein structures. Many hormone−/abiotic stress-responsive cis-acting elements and transcription factor binding sites were identified in BnaNPF promoter regions. Chromosome distribution analysis indicated that BnaNPFs within a subfamily tend to cluster on one chromosome. Syntenic relationship analysis showed that allotetraploid creation by its ancestors (Brassica rapa and Brassica oleracea) (57.89%) and small-scale duplication events (39.85%) contributed to rapid BnaNPF expansion in B. napus. A genome-wide spatiotemporal expression survey showed that NPF genes of each Arabidopsis and B. napus subfamily have preferential expression patterns across developmental stages, most of them are expressed in a few organs. RNA-seq analysis showed that many BnaNPFs (32.66%) have wide exogenous hormone-inductive profiles, suggesting important hormone-mediated patterns in diverse bioprocesses. Homologs in a clade or branch within a given subfamily have conserved organ/spatiotemporal and hormone-inductive profiles, indicating functional conservation during evolution. qRT-PCR-based comparative expression analysis of the 12 BnaNPFs in the NPF2–1 subfamily between high- and low-glucosinolate (GLS) content B. napus varieties revealed that homologs of AtNPF2.9 (BnaNPF2.12, BnaNPF2.13, and BnaNPF2.14), AtNPF2.10 (BnaNPF2.19 and BnaNPF2.20), and AtNPF2.11 (BnaNPF2.26 and BnaNPF2.28) might be involved in GLS transport. qRT-PCR further confirmed the hormone-responsive expression profiles of these putative GLS transporter genes. Conclusion We identified 199 B. napus BnaNPFs; these were divided into 11 subfamilies. Allopolyploidy and small-scale duplication events contributed to the immense expansion of BnaNPFs in B. napus. The BnaNPFs had preferential expression patterns in different tissues/organs and wide hormone-induced expression profiles. Four BnaNPFs in the NPF2–1 subfamily may be involved in GLS transport. Our results provide an abundant gene resource for further functional analysis of BnaNPFs. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-020-07274-7.

Brassica napus is a significant source of human-edible vegetable oil and animal protein feed; thus, it is an essential oil crop, extensively cultivated in Asia, North America, and Europe. Given essential roles in plant nitrate, di−/tripeptide, hormone, potassium, chloride, and secondary metabolite transports, NPFs have been systematically identified and analyzed in many species, including Arabidopsis [20], Oryza sativa (rice) [24], Triticum aestivum (wheat) [25], and Malus domestica (apple) [26] at the genome-wide level. Identifying and analyzing this gene family in the B. napus genome will provide a solid foundation for exploring its potential roles in transporting nitrate, hormones, and GLS, among others.
This study identified NPFs in the B. napus genome, accompanied by comprehensive analysis of their gene and protein structural features, chromosomal location, classification, promoter regulation network, and genomic duplication mechanism. Further, we performed systematic expression profile analysis of this gene family in diverse tissues across different developmental stages in Arabidopsis (79 tissues) and B. napus (50 tissues). Additionally, expression patterns of NPF gene family in B. napus under five exogenous hormone inductions (IAA, auxin; ABA, abscisic acid; GA 3 , gibberellic acid; 6-BA, cytokinin; and ACC, ethylene) were assessed, based on the RNA-Seq dataset. Moreover, expression patterns of 12 candidate NPFs of the NPF2-1 subfamily in one high-and one low-GLS B. napus variety, as well as their expression profiles under hormone induction, were assessed using qRT-PCR. Our study provides an abundant gene resource for further functional analysis of NPFs in B. napus.

Protein characteristics and intron pattern diversity
Based on multiple alignment analysis of the 155 fulllength BnaNPFs with relatively complete coding regions, the protein sequence feature was further explored. TMs and other conserved protein domains were predicted using the HMMER software (http://www.ebi.ac.uk/Tools/ hmmer) [36].
At the nucleic acid sequence level, we further analyzed the intron insertion site, number, and phase of candidate BnaNPFs by using the Gene Structure Display Server (GSDS) 2.0 (http://gsds.gao-lab.org/) [38]. Our results showed that all 155 BnaNPFs contained 1-16 introns, and 86.45% (134/155) of BnaNPFs had 3-5 introns (Additional file 6: Figure S4). Notably, three of these introns were highly conserved in almost all BnaNPFs in terms of insertion sites and phases; one intron was inserted ahead of the PTR2 domain, and two introns were inserted within the PTR2 domain (one in the third TM and another between the sixth and seventh TMs) (Additional file 6: Figure S4). This finding suggests that these three introns may be necessary for the function of BnaNPFs. Moreover, apart from these three introns, the other introns were commonly conserved within each subfamily or clade, but were less conserved among distinct subfamilies (Additional file 6: Figure S4). Furthermore, we found that the intron insertion sites and phases of BnaNPFs and AtNPFs were highly conserved in each clade or subfamily (Additional file 6: Figure S4), indicating conserved structural features during their evolution.
Overall, the conserved protein and gene sequences strongly support our subfamily division based on phylogenetic analysis.
In summary, our results reveal that expression of BnaNPFs may be regulated by various kinds of hormones, abiotic stresses, and TFs.

Chromosomal location and syntenic relationship in BnaNPFs
The distribution of BnaNPFs on B. napus chromosomes was analyzed based on genomic annotation information obtained from the GENOSCOPE database (http://www. genoscope.cns.fr/brassicanapus/) [27]. As shown in Fig. 3a, most of the 199 BnaNPFs were mapped on the 19 chromosomes; however, the exact locations of 6 genes in A n subgenome and 29 genes in the C n subgenome were unclear (Additional file 1: Table S1). The numbers of BnaNPFs in A n (95) and C n (104) subgenomes were similar. However, the distribution of BnaNPFs on different chromosomes was uneven. For example, A03, A04, and C01 contained only three genes, while A07 had up to 22 genes (Fig. 3a). Notably, BnaNPFs belonging to the same subfamily tended to cluster on several chromosomes: 39.44% (28/71) of NPF5-1 subfamily members were distributed on the A02, A07, A09, and C02 chromosomes (Fig. 3a). Similar trends in the NPF gene family were observed in Arabidopsis, B. rapa, and B. oleracea. In Arabidopsis, all members of NPF2-2 subfamily (AtNPF2.1-AtNPF2.7) were clustered on the 03 chromosome, and 46.67% (7/15) of NPF5-1  subfamily (AtNPF5.10-AtNPF5.16) members were clustered on the 01 chromosome (Additional file 9: Figure  S5a); Similarly, members of NPF5-1 subfamily were distributed mainly on the A07 chromosome in B. rapa (Additional file 9: Figure S5b) and C07 chromosome in B. oleracea (Additional file 9: Figure S5c). These results suggest that NPFs in the same subfamily tend to assemble as gene clusters, and this trend may be conserved in plants.
The collinearity of NPFs in B. oleracea, B. rapa, and B. napus genomes was analyzed using the CoGe tool (https://genomevolution.org/CoGe/) [42] to explore the expansion mechanism of BnaNPFs. Our results show that 133 of 199 BnaNPFs in the B. napus genome had a syntenic relationship ( Fig. 3b; Additional file 10: Table S5). All 133 genes had a collinear relationship with BraNPFs, while 127 BnaNPFs had a collinear relationship with BolNPFs. We further speculated that 49 of the 133 genes (36.84%) were inherited from B. rapa, and 28 genes (21.05%) were inherited from the B. oleracea genome, based on the syntenic relationship between the descendant and its ancestors. Given that B. napus evolved by hybridization between B. oleracea and B. rapa~7500 years ago, it was evident that allopolyploidy (57.89%) heavily contributed to the rapid expansion of NPFs in B. napus. Moreover, gene loss following allopolyploidy was biased; the NPFs inherited from B. rapa were inclined to be retained. Furthermore, 39.85% (56/133) of genes originated from other duplication events within the B. napus genome, including 30 genes from segmental exchange (SE), 21 genes from segmental duplication (SD), and 5 genes from homologous exchange (HE) events. These results proved that small-scale duplication events (including HE, SE, and SD) also contributed to the massive expansion of NPFs in B. napus, especially the SE and SD events. Notably, of the 21 genes that underwent SD events, 15 were derived from B. rapa, while the remaining 6 were inherited from B. oleracea; this indicates that the genes from B. rapa tended to undergo SD in B. napus. Regarding the HE event, three of the five HE genes were from the A n subgenome, which replaced the genes in the C n subgenome. This finding confirmed that the A n subgenome replaced more of the C n subgenome after allopolyploidy and featured more dominantly in each chromosome [43]. Three pairs of putative tandem duplication (TD) genes (BnaNPF2.26/BnaNPF2.27, BnaNPF4.7/ BnaNPF4.9, and BnaNPF5.22/BnaNPF5.23) were observed, based on their chromosome distribution and sequence similarity.
Overall, our results indicate that allopolyploidy and small-scale duplication events (including SE, SD, and HE) are the primary driving force for the rapid expansion of BnaNPFs in B. napus, and that those derived from B. rapa tended to be retained during evolution.

Comparative expression analysis of AtNPFs and BnaNPFs across plant development
As gene expression pattern is an essential clue as to its function, in order to explore gene expression patterns as well as expression and function similarity between different species, we analyzed and compared global expression profiles of AtNPFs and BnaNPFs in different tissues and organs at distinct developmental stages. We used public expression datasets of Arabidopsis (http://bar.utoronto. ca/efp/cgi-bin/efpWeb.cgi) [44] and B. napus (BioProject ID PRJNA358784).
In B. napus, with the exception of 67 BnaNPFs having no detectable expression values (FPKM < 1) that were excluded from analysis (Additional file 11: Table S6), most (132/199) of the remaining genes had preferential expression profiles in the 50 tissues of seven organs (root, stem, leaf, hypocotyl, flower, silique pericarp, and seed) at six developmental stages ( Fig. 4 and Additional file 11: Table  S6). For instance, members of the NPF1 subfamily had higher transcriptional levels in root, stem, hypocotyl, flower, and silique pericarp; members of NPF7 were highly expressed in flowers, silique pericarp, and seeds; and members of NPF2-1 were mainly expressed in flower and seed tissues (Fig. 4). In general, expression patterns were conserved in each subfamily or each clade within a subfamily, but were quite different across different subfamilies, suggesting the expression differentiation trend of this gene family. For example, expression patterns of NPF1, NPF2-2, NPF3, and NPF6-1 subfamilies were similar in each subfamily, while the expression profile of the NPF2-1 subfamily was classified into three conserved patterns that were consistent with the three major clades in this subfamily. Additionally, we found that 40% (6/15) of the BnaNPFs expressed explicitly in seeds belong to the NPF2-1 subfamily, and 33.33% (5/15) belong to the NPF4 subfamily, suggesting essential roles for these two subfamilies in seed development.
In Arabidopsis, consistent with the situation in B. napus, most of the AtNPFs had preferential expression patterns in the organs investigated (Additional file 12: Figure S6). Members of the NPF2-2 subfamily (AtNPF2.3, AtNPF2.4, AtNPF2.5, and AtNPF2.7) were preferentially expressed in roots; AtNPF4.1 and AtNPF4.5 in the NPF4 subfamily were mainly expressed in seeds; and AtNPF2.10 and AtNPF2.11 in the NPF2-1 subfamily had higher expression levels in roots, stems, leaves, and flowers. Notably, the expression patterns of homologs in both species were generally conserved. Members of NPF2-1 in B. napus and Arabidopsis were preferentially expressed in flower and seed organs, and members of NPF2-2 in these two species were preferentially expressed in roots. Given that genes with similar expression patterns may share similar functions, the homologs may have similar/conserved functions in Arabidopsis and B. napus.

Expression profiles of BnaNPFs under different hormone inductions
As mentioned (Fig. 2), many hormone-responsive CREs were observed in the promoter regions of candidate BnaNPFs, suggesting possible roles for plant hormones in BnaNPF expression. Therefore, we analyzed the expression profiles of BnaNPFs under five exogenous hormone treatments (IAA, ABA, 6-BA, GA 3 , and ACC) in B. napus seedling roots, based on the RNA-seq data (BioProject ID: PRJNA608211).
Results showed that 32.66% (65/199) of the BnaNPFs were upregulated by one or more types of hormones (Fig. 5). With the exception of the NPF5-2 subfamily, which had no detectable expression level, the expression patterns of the other 10 subfamilies were induced at different levels under hormone treatments. NPF1 subfamily members were positively induced by ACC induction; 42.86% (3/7) of NPF3 subfamily members were upregulated by ABA and ACC treatments; and 90.91% (10/11) of NPF6-1 subfamily members had higher expression  Table S6 levels under the five hormone treatments. In contrast, 35.18% (70/199) of BnaNPFs were downregulated by these five hormone treatments (Fig. 5) (Fig. 5), implying their functional conservation.
Overall, the expression of many BnaNPFs was sensitive to exogenous hormone induction, suggesting that the essential roles of this gene family in diverse plant processes may be regulated by hormone-mediated patterns.
Expression of BnaNPFs in high-and low-GLS content B. napus varieties GLS is a class of important secondary metabolites found in Brassicaceae that have distinctive benefits for plant defense and human nutrition (such as inhibiting carcinogen activation) [45,46]. Recently, several Arabidopsis NPFs, including AtNPF2.10, AtNPF2.11, and AtNPF2.9, were shown to be involved in GLS transport [20,47]. Phylogenetic analysis showed that these proteins and their B. napus homologs (20 proteins, BnaNPF2.12-BnaNPF2.30 and BnaNPF2.40) were clustered into the NPF2-1 subfamily with conserved sequence features (Fig. 1) and expression profiles (Fig. 4), implying that the B. napus homologs may have similar roles in GLS transport. To confirm the possible roles of these BnaNPFs in GLS transport, we further compared the expression profiles of 12 GLS-coding genes between low-and high-GLS content B. napus varieties, Zhongshuang 11 (ZS11; Fig. 6a and Additional file 13: Table S7) and Zhongyou 821  Fig. 6b and Additional file 13: Table S7), using qRT-PCR [48]. Eight genes in this subfamily that have no detectable expression by RNA-Seq analysis may be pseudogenes, and thus were excluded from this analysis.
Consistent with RNA-seq results (Fig. 4), all 12 candidates were preferentially expressed in a few organs at different developmental stages; however, one gene (BnaNPF2.18) was not expressed in any of the samples investigated. Moreover, most of the other 11 candidates exhibited different expression patterns in the root, stem, leaf, flower, silique, and seed organs between the two varieties ( Fig. 6a,b). For example, the three homologs of the AtNPF2.9 gene (BnaNPF2.12, BnaNPF2.13, and BnaNPF2.14) were highly expressed in seed tissues of ZY821 but had no detectable expression levels in the seed tissues of ZS11. Similarly, although the expression levels were relatively lower, two homologs of AtNPF2.11 (BnaNPF2.26 and BnaNPF2.28) were expressed in seed tissues of ZY821 but were not expressed in those of ZS11. In contrast, the homologs of AtNPF2.10 (BnaNPF2. 19 and BnaNPF2.20) were highly expressed in seed tissues of ZS11 but not in those of ZY821. It was previously reported that AtNPF2.9 is the typical indole-specific GLS transporter gene, whereas AtNPF2.11 is the transporter gene for both indole and aliphatic GLS [20,37]. Given the fact that GLS content in ZS11 seeds is significantly lower than in ZY821, our results imply that these five genes are involved in GLS transport in B. napus seed tissues, especially the three more highly expressed in ZY821 seed tissues.
Notably, NPFs in most subfamilies can generally transport more than one type of substrate (Table 1). For example, NPF2-1 subfamily members can transport four types of substrates: nitrate, GLS, GA, and MeJA ( Fig. 1 and Table 1); in Arabidopsis, AtNPF6.3 of the NPF6-1 subfamily transports nitrate as well as auxin [19]. Moreover, many studies have found that the role of NPFs in transporting diverse substrates generally demonstrates hormone-mediated characteristics. For example, AtNPF6.3 of the NPF6-1 subfamily was highly induced by IAA treatment under low nitrogen conditions [72]. AtNPF2.4 was repressed by ABA treatment and then played a role in chlorine transport [21]. AtNPF2.10 was upregulated by MeJA treatment, which then accelerated the transport of gibberellin [54]. AtNPF3.1 was upregulated by ABA to promote the transport of gibberellin [55]. AtNPF5.2 was regulated by SA, MeJA, and ABA treatments against biotic and abiotic stresses [60]. In this study, many CREs involved in hormone responses, such as SA-(103/199 genes), ABA-(93/199 genes), and MeJA-responsive CRE (123/199 genes), were found in a series of BnaNPF promoters (Fig. 2), suggesting their potential hormone-inducing characteristics. Accordingly, the expression of 32.66% of the BnaNPFs (65/199 genes) were regulated by one or more types of hormone inductions (ABA, IAA, 6-BA, GA 3 , and ACC) (Fig. 5). Consistent with previous work [72], NPF6-1 subfamily proteins in B. napus (e.g., BnaNPF6.5, BnaNPF6.7, and BnaNPF6.7) were also highly induced by IAA in our study. Additionally, we revealed that the genes involved in GLS transport are induced by IAA, GA 3 , and 6-BA treatments (Fig. 6). Together, these results support the hypothesis that hormones have an essential role in substrate transport by NPFs.
Given that the role of NPFs in transporting many substrates is crucial for plant development and stress response, genome-wide analyses of the NPF gene family have been performed in many plant species. However, the classification of this gene family is not yet uniform. For example, this gene family was divided into 10 supergroups and 32 groups (subfamilies) based on phylogenetic analysis of 20 plant genomes [73]. In contrast, other research divided this gene family into 8 subfamilies (NPF1-NPF8) based on similar analysis in 33 plant genomes (including Physcomitrella patens and Selaginella moellendorffii) [7]. Subsequent studies generally followed the criteria of the latter division [25,26,74]. Recently, in apple, the NPF2 subfamily was further divided into two groups  [26], implying a new classification trend. In this study, we found that three of the eight previously demonstrated subfamilies (NPF2, NPF5, and NPF6) [7] should be divided into two subfamilies with high bootstrap values: NPF2-1/NPF2-2, NPF5-1/NPF5-2, and NPF6-1/ NPF6-2, respectively (Fig. 1). To confirm this result, we further expanded our dataset to include NPFs from O. sativa, Populus trichocarpa, Z. mays, B. rapa, B. oleracea, and Glycine max (Additional file 2: Table S2). Phylogenetic analysis of the NPFs from these different species highly supported that of BnaNPFs (Additional file 15: Figure S7). Moreover, the gene and protein structures (Additional file 4: Figure S2, Additional file 5: Figure S3, and Additional file 6: Figure S4) and expression patterns (Figs. 4 and 5) of BnaNPFs in each subfamily supported our classification, as well. Interestingly, all the NPFs involved in secondary metabolite transport known to date belong to the NPF2-1 subfamily (Table 1). We confirmed that 267 proteins belonging to this subfamily exist in 31 angiosperms, though not in the lower plants P. patens and S. moellendorffii [7] (Additional file 16: Table S9). Given that currently known members of the NPF2-1 subfamily across different plant species, including lower plants (Table 1), are mainly involved in nitrate transport, we speculated that the secondary metabolite transport feature of this subfamily was newly evolved in a given lineage or species in angiosperms during their evolution, indicating the specific subfunctionalization trend of this gene family.

Conclusions
In this study, 199 BnaNPFs were identified in the B. napus genome and divided into 11 subfamilies having conserved gene and protein structures within each subfamily or clade. The allopolyploidy produced by its ancestors and the small-scale duplication events in B. napus acted as the primary driving forces for the massive expansion of this gene family in B. napus. Genes derived from B. rapa were retained after the allopolyploidy event during B. napus evolution. Most of the BnaNPFs were likely to be preferentially expressed in a few tissues or organs, and these expression profiles were commonly conserved in each subfamily or in each clade within a subfamily. Hormone inductions regulated the expression of many BnaNPFs. Five genes (BnaNPF2.12, BnaNPF2.13, BnaNPF2.14, BnaNPF2.26, and BnaNPF2.28) in the NPF2-1 subfamily may be involved in GLS transport in B. napus, mediated by IAA, ACC, GA 3 , or 6-BA.

Chromosomal location and collinearity synteny analysis
The chromosome information of the 199 candidate BnaNPFs was obtained from the GENOSCOPE database. The collinear relationship of the candidate NPFs in B. oleracea, B. rapa, and B. napus genomes was assessed using the CoGe online software (https://genomevolution. org/CoGe/) [42] with default parameters. Duplication events of candidate NPFs were defined according to the method used in our previous report [76]. Based on cross-genome collinearity analysis, the species with the maximum orthologous blocks/most closer colinear relationships (including the NPF orthologous gene pairs) are considered the progenitors of BnaNPFs. HE, SE, and SD events were distinguished from each other based on chromosomal homology and colinear relationship (orthologous gene pairs in orthologous blocks) of the A n (derived from B. rapa) and C n (derived from B. oleracea) subgenomes and their respective progenitor genomes (B. rapa and B. oleracea) in all possible combination pairs. The chromosome map of candidate BnaNPFs was drawn by using the MapChart software with default parameters [77].

Expression profile analysis of NPFs in Arabidopsis and B. napus
The expression profile of AtNPFs (including root, stem, leaf, apex, flower, and seed) was obtained from AtGenExpress (http://weigelworld.org/resources.html) [43]. The RNA-seq data, including 50 tissues of seven B. napus variety ZS11 organs (root, stem, leaf, flower, seed, and silique) at six developmental stages (seed germination, seedling, budding, initial flowering, full-bloom, and seed maturation), were obtained from NCBI (BioProject ID PRJNA358784). Expression profiles of candidate BnaNPFs in ZS11 seedling roots under five exogenous hormone treatments (IAA, GA 3 , 6-BA, ABA, and ACC) were extracted from our RNA-seq dataset (BioProject ID: PRJNA608211). BnaNPFs with FPKM ≥1 were retained, and the FPKM values of candidates were log2transformed for visualization by the R package [78]. The heatmap was combined with hierarchical clustering methods of the log2-transformed RNA-Seq data [79].

Plant material and hormone treatment
The seeds of a high-GLS content B. napus variety (Zhongyou 821; ZY821) and a low-GLS content variety (Zhongshuang 11; ZS11) were grown in Beibei (Chongqing, China) with standard agronomic procedures to analyze the temporal and spatial expression patterns of BnaNPFs involved in GLS transport. The root (Ro-s), stem (St-s), and leaf (Le-s) tissues at the seedling stage; the root (Ro-f), stem (St-f), leaf (Le-f), and flower (Fl-f) tissues at the flowering stage; and the root (Ro-ss), stem (St-ss), leaf (Le-ss), siliques 3 days after pollination (Si-3d), siliques 7 days after pollination (Si-7d), siliques 15 days after pollination (Si-15d), the seed 15 days after pollination (Se-15d), and seed 45 days after pollination (Se-45d) at the mature stage were collected from both ZS11 and ZY821. All tissues were immediately frozen in liquid nitrogen and stored at − 80°C for RNA isolation.
To further analyze the expression patterns of the GLS transporter genes under five exogenous hormone inductions (ABA, IAA, GA 3 , 6-BA, and ACC), ZS11 seeds were germinated and cultivated at 25°C under a 16/8 h (day/ night) photoperiod in an artificial climate chamber. At the three-leaf stage, seedlings were transferred to Hoagland solution and were further cultured to the five-leaf stage. Seedlings were then treated with Hoagland's liquid medium, containing 10 μmol/L IAA, 50 μmol/L ABA, 75 μmol/L 6-BA, 120 μmol/L GA 3 , and 60 μmol/L ACC. Root tissues were harvested at 0, 1, 3, 6, 12, and 24 h after each treatment. All samples were quickly frozen in liquid nitrogen and then stored at − 80°C for RNA isolation.
Expression analysis of NPFs in B. napus using qRT-PCR The expression profile of 12 putative GLS transporters encoding BnaNPFs of the NPF2-1 subfamily in different tissues and under five exogenous hormone inductions was analyzed via qRT-PCR, using BnaActin7 (GenBank accession no. AF024716) and BnaUBI (GenBank accession no. NC027770) as the reference genes. Primer pairs for qRT-PCR analysis were designed using Primer Premier 5 (Additional file 17: Table S10).
Total RNA was extracted from each of the samples with the EASYspin Total RNA Extraction Kit (Biomed, Beijing). The concentration and quality of total RNA for each sample were confirmed through gel electrophoresis analysis and NanoDrop 2000 spectrophotometer measurement.
Potentially contaminating DNA was eliminated by DNase I (Promega, USA). The cDNA of each sample was synthesized using the M-MuLV RT kit (Takara Biotechnology, Japan). The real-time PCR analysis (qRT-PCR) was performed using the SYBR-Green PrimeScript RT-PCR Kit (Takara Biotechnology, Japan) with the CFX Connect™ Real-Time System (Bio-Rad, Chongqing, China). The parameters of qRT-PCR were as follows: 95°C for 3 min (initial denaturation), followed by 40 cycles of 95°C for 10s (denaturation) and 58°C for 30s (annealing). Each PCR was validated in three independent repeat experiments. The qRT-PCR results were calculated using the 2 -ΔΔct method [47]. Expression values were log2-transformed and visualized with the R package [78].