Genome-wide sequence identification and expression analysis of ARF family in sugar beet (Beta vulgaris L.) under salinity stresses

Materials

The research materials used in this experiment for salinity stress were “O68” strain, which was the salt-tolerant sugar beet strain selected by the laboratory (Shi et al., 2008). In the sugar beet seedling stage, the sugar beet seeds were first soaked in running water for 12 h and then sown into the wet sponge which was placed in the incubator for 2 days (dark, 24 °C). Sugar beet seedlings were taken out after they germinated and placed in the culture pot at 22 ± 1 °C with 16 h/d supplemental lighting, followed by 8h/d darkness, and cultured with clear water until they grow a third pair of true leaves, then they were treated with salinity stress. The salt stress conditions were that the nutrient solution was replaced with a 300 mM NaCl solution for 24 h, the remaining conditions were unchanged, and the control plants were set at the same time without salt treatment. After the stress treatment, leaves and roots were collected. And they were wrapped in tin foil, quickly placed in liquid nitrogen for precooling. The frozen samples were uniformly stored in the ultra-low temperature refrigerator (−80 °C).

BvARF identification in sugar beet

The whole genome data of sugar beet was published by Max Planck Institute for Molecular Genetics in 2014 (http://bvseq.molgen.mpg.de/index.shtml). And the assembly RefBeet-1.1 of KWS2320 was downloaded for bioinformation research. It includes functional annotation, genome sequence, transcript sequences, and protein sequences etc. ARF proteins conserved domain (PF06507) on Pfam were used as a seed sequence to search the whole genome data, and e-value  <  1e⁻⁵ was set on HMMER (http://www.hmmer.org/). Pfam online tool performed a conserved domain analysis of candidate BvARF protein sequences (Finn et al., 2016), screening protein sequences containing the B3 domain (PF02362), AUX_RESP domain (PF06507) and AUX/IAA (PF02309) domains to remove redundant proteins, and remainders were considered to be BvARF proteins. DNAMAN7.0 was used for multiple sequences alignment of the full-length coding sequence of BvARF proteins, and the conserved domains in the BvARF protein sequences were identified by Weblogo2.8.2 (http://weblogo.berkeley.edu/logo.cgi).

Bioinformatics analysis of BvARF family

Molecular weight and isoelectric point of BvARF proteins were predicted by ExPASY (https://web.expasy.org/protparam/). MapInspect was used to map the position of the BvARF genes on chromosomes. Exon and intron structures of the BvARF genes were analyzed by GSDS (http://gsds.cbi.pku.edu.cn/Gsds_about.php) (Hu et al., 2015a). ClustalX was used to carry out a multiple sequence alignment of BvARF proteins (Larkin et al., 2007). A phylogenetic tree of BvARF family and AtARF family and a phylogenetic tree of BvARF family were constructed by the neighbor-joining with bootstrap replicate set to 1000, and the other parameters default. And the protein sequences of AtARF for comparison were provided in File S1. BvARF proteins were submitted to MEME (http://meme-suite.org/tools/meme) in the order of homology and subjected to motif analysis (Bailey et al., 2006), in which the number of expected motifs was set to 20, and the rest parameters were all default values. Subcellular localization of BvARF proteins was predicted by using CELLO (http://cello.life.nctu.edu.tw/). Gene Ontology (GO) categories of genes were obtained from PlantTFDB (http://planttfdb.cbi.pku.edu.cn/), and GO level 2 of BvARF genes were visualized with WEGO (http://wego.genomics.org.cn/). And the GO annotation numbers of BvARF were provided in File S2.

Expression analysis of BvARF under salinity stress

Total RNA was isolated from the aforementioned sugar beet samples using MiniBEST Plant RNA Extraction Kit (TaKaRa, Japan) and reverse transcribed into cDNA using High Capacity cDNA Reverse Transcription Kit (applied biosystems, USA). The patterns of expression of these genes under normal growth conditions and in response to salinity stresses were analyzed by quantitative real-time PCR (qRT-PCR). Primers for qPCR were designed using Primer 5 and NCBI Primer-BLAST, synthesized by Huada Gene Co., Ltd., and the primer sequences are shown in Table 1. The ICDH (Isocitrate dehydrogenase) gene was used as an internal control.

For quantitative analysis using the CFX96 Real-Time System (BIO-RAD, USA), and iTaq Universal SYBR Green Supermix kit (BIO-RAD, USA) was used with ICDH as the housekeeping gene. The reaction system was 10 µL, in which the dye was 5 µL, the ddH2O was 3.4 µl, the forward primer and reverse primers were each 0.4 µL, and the cDNA was 0.8 µL, and each treatment was repeated three times. Data analysis was performed by calculating 2^−ΔΔT, and the relative expression amount of each gene was expressed by mean ± standard deviation.

Identification of BvARF family members

Seventeen candidate BvARF proteins were obtained by using PF06507 as a seed sequence to perform HMMER alignment in the sugar beet genome database. And the sequences of 17 candidate BvARF proteins were analyzed by Pfam for conserved domain, and the proteins containing no specific domains of ARF were removed. A total of 16 BvARF proteins were identified, as shown in Fig. 1.

It can be seen from the figure that all BvARF have a C-terminal B3 binding domain and an intermediate AUX_RESP domain, while the AUX/IAA domain exists only in 11 BvARF except for opag, okdq, ghtj, jrpi, and eqms.

From the comparison results of DNAMAN7.0, the B3 DNA binding domain is highly conserved and occasionally modified. The AUX_RESR domain is not as conserved as B3, the AUX/IAA domain is the least conserved and in short length, as shown in Fig. 2.

In the 16 BvARF proteins, gcik.ti and gcik.t2 were both translation products of transcript variants of BvARF gene gcik. And the efdx.t1 and efdx.t2 were both translation products of transcript variants of BvARF gene efdx. Therefore, a total of 14 BvARF gene loci were obtained and indicated by the name in the sugar beet genome-wide database.

The 16 BvARF proteins were blastp aligned with the NCBI database and all found to be members of the BvARF family, as shown in Table 2. No more new BvARF was identified. Descriptions in NCBI suggest that jrpi, opag and zzhs may have isoforms.

Analysis of physicochemical properties of BvARF proteins

Sequence analysis showed that the average length of the coding region of ARF gene was 2,390 bp (1,734–3,429 bp), the average number of amino acids encoding proteins was 796 (577–1,142), and the average molecular weight predicted by ExPASY was 88.7 kDa (64,060.51–127,711.63 Da), the isoelectric point averages 6.13 (5.21–7.9), as shown in Table 3.

Chromosomal localization of BvARF gene

The location analysis showed that 11 BvARF genes were distributed in 8 chromosomes except chromosome 3. There are 2 BvARF on chromosomes 1, 8, and 9, and one on chromosomes 2, 4, 5, 6, and 7, as shown in Fig. 3. The location information of other three genes hgze, jrpi, orwr is unknown in the database, that may be caused by these BvARF genes were located in the scaffolds which could not be assembled in the chromosomes.

Phylogenetic relationships and gene structures analysis of BvARF

To survey the evolutionary relationships between BvARF proteins from sugar beet and previously-reported ARF proteins from the dicot Arabidopsis we constructed an unrooted phylogenetic tree with 16 BvARF and 23 Arabidopsis ARF protein sequences. The results are shown in Fig. 4, based on existing descriptions of ARF conserved domains in Arabidopsis (Hagen & Guilfoyle, 2002), and all 39 ARF proteins were divided into six categories.

The full-length sequence of the 16 BvARF proteins screened by MEGA7 was used to make a rootless tree by the neighbor-joining method, and the bootstrap procedure was tested with 1000 bootstrap repetitions. The results are shown in Fig. 5, which are also divided into six categories: Class I (AtARF1/2-like), Class II (AtARF3/4-like), Class III (AtARF5-like), Class IV (AtARF7/19-like), Class V (AtARF6/8-like) and Class VI (AtARF10/16/17-like). Four members were assigned to Class I (qzmp, orwr, kddw and yzaj) 2 to Class II (opag and efdx), Class III only contains hgze, 2 to Class IV (qfwi and zzhs), 2 to Class V (okdq and gcik), and 3 to Class VI (ghtj, jrpi and eqms). The intron/exon structure analysis of the BvARF genes can provide valuble imformation concerning evolutionary relationships among taxa. Analysis of these sequences revealed that all BvARF genes contained introns in the coding sequence, and the number of exons varied from 3 to 15. In general, most of the BvARF members within a given group possessed a similar exon/intron structure in terms of intron numbers and exon length.

Motifs analysis and subcellular localization prediction of BvARF proteins

For the motifs analysis of the BvARF proteins, the expected number of motifs was set to 20. The motifs were sorted according to the e-value from small to large, and the protein sequences were ranked according to the reliability of the predicted motifs. The results are shown in Figs. 6 and 7. Some motifs were common to all BvARF proteins, such as motifs 1, 2, 3, and 5, and some motifs were specific to certain proteins, such as motifs 13, 18, and 20. Comparing the results of the motifs analysis with the conserved domains of the previous Pfam, it was found that the motif 1, 2 constituted the B3 domain, the motif 6, 8, 11, 12 constituted AUX_RESP, and the motif 4, 9, 14 constituted AUX/IAA. Previous studies have shown that ARF activates or inhibits target genes in relation to amino acids in the middle region (Tiwari, 2003; Ulmasov, Hagen & Guilfoyle, 1999), and it is worth noting that phantom 17 has a dense glutamine (Q) composition. It has shown in previous research that when glutamines gather in the middle region of the ARF proteins were generally considered to have activation domain activity (Shen et al., 2015). It plays a transcriptional activation role in the transcriptional regulation of auxin, and further predicts that qfwi, zzhs, okdq, gcik may have transcriptional activation which were all in Class IV and Class V.

The subcellular localization prediction of CELLO on BvARF proteins is shown in Table 4. The results showed that the most of BvARF proteins were only scored in the nucleus, that is, they were all more reliably located in the nucleus. The scores of hgze and jrpi were high in cytoplasm and nucleus. Subcellular localization results indicate that almost all BvARF proteins play a regulatory role in the nucleus.

Gene ontology annotation

To survey the functions of the BvARF, GO annotation was obtained from PlantTFDB to construct GO graphs. As shown in Fig. 8, all BvARF (14, 100%) are involved in cell part, cell and organelle. The results also showed that the BvARF were involved in diverse biological processes and predominantly participated in metabolic process, biological regulation, regulation of biological process, cellular process(12, 85.7%), binding(11, 78.6%) and response to stimulus(8, 57.1%). The GO terminologies of BvARF were relatively concentrated.

Expression profiles of BvARF genes in sugar beet under normal growth conditions and in response to salinity stress

In recent years, it has reported that auxin acts as integral part of plant response to salinity stress. Thus we compare the levels of expression of BvARF in condition with and without salt treatment. Expression patterns of 14 BvARF in leaves and root were analyzed by CFX96 Real-Time System to clarify the roles of BvARF genes in salt tolerance. The results are shown in Figs. 9 and 10. In the leaves, except for jrpi, yzaj and eqms, other BvARF genes all showed different degrees of up-regulation or down-regulation. Most of them showed a down-regulation, and in which the degree of 6 BvARF were significantly reduced (ghtj, okdq, opag, zzhs and efdx). And a few BvARF genes were up-regulated, with only upward trend of gcik was obvious. In the root, after treatment with salinity stress, except for qfwi, orwr and eqms, other BvARF genes all showed different degrees of up-regulation or down-regulation. Most of the BvARF genes showed an upward trend, and in which 6 BvARF showed an obvious upward trend (hgze, okdq, yzaj, kddw, gcik and efdx). And 5 BvARF genes showed a downward trend, but the degree of down-regulation was not significant (ghtj, qzmp, qfwi, jrpi and eqms).