Digital Gene Expression Analysis of Corky Split Vein Caused by Boron Deficiency in ‘Newhall’ Navel Orange (Citrus sinensis Osbeck) for Selecting Differentially Expressed Genes Related to Vascular Hypertrophy

Abstract

Corky split vein caused by boron (B) deficiency in ‘Newhall’ Navel Orange was studied in the present research. The boron-deficient citrus exhibited a symptom of corky split vein in mature leaves. Morphologic and anatomical surveys at four representative phases of corky split veins showed that the symptom was the result of vascular hypertrophy. Digital gene expression (DGE) analysis was performed based on the Illumina HiSeq™ 2000 platform, which was applied to analyze the gene expression profilings of corky split veins at four morphologic phases. Over 5.3 million clean reads per library were successfully mapped to the reference database and more than 22897 mapped genes per library were simultaneously obtained. Analysis of the differentially expressed genes (DEGs) revealed that the expressions of genes associated with cytokinin signal transduction, cell division, vascular development, lignin biosynthesis and photosynthesis in corky split veins were all affected. The expressions of WOL and ARR12 involved in the cytokinin signal transduction pathway were up-regulated at 1st phase of corky split vein development. Furthermore, the expressions of some cell cycle genes, CYCs and CDKB, and vascular development genes, WOX4 and VND7, were up-regulated at the following 2nd and 3rd phases. These findings indicated that the cytokinin signal transduction pathway may play a role in initiating symptom observed in our study.

About the Authors:

Cheng-Quan Yang

* E-mail: [email protected] (CQY); [email protected] (SAP)

Affiliation: College of Horticulture and Forestry Sciences, Huazhong Agricultural University, Key Laboratory of Horticultural Plant Biology, Ministry of Education, Wuhan, China

Yong-Zhong Liu

Affiliation: College of Horticulture and Forestry Sciences, Huazhong Agricultural University, Key Laboratory of Horticultural Plant Biology, Ministry of Education, Wuhan, China

Ji-Cui An

Affiliation: College of Horticulture and Forestry Sciences, Huazhong Agricultural University, Key Laboratory of Horticultural Plant Biology, Ministry of Education, Wuhan, China

Shuang Li

Affiliation: College of Horticulture and Forestry Sciences, Huazhong Agricultural University, Key Laboratory of Horticultural Plant Biology, Ministry of Education, Wuhan, China

Long-Fei Jin

Affiliation: College of Horticulture and Forestry Sciences, Huazhong Agricultural University, Key Laboratory of Horticultural Plant Biology, Ministry of Education, Wuhan, China

Gao-Feng Zhou

Affiliation: College of Horticulture and Forestry Sciences, Huazhong Agricultural University, Key Laboratory of Horticultural Plant Biology, Ministry of Education, Wuhan, China

Qing-Jiang Wei

Affiliation: College of Horticulture and Forestry Sciences, Huazhong Agricultural University, Key Laboratory of Horticultural Plant Biology, Ministry of Education, Wuhan, China

Hui-Qing Yan

Affiliation: College of Horticulture and Forestry Sciences, Huazhong Agricultural University, Key Laboratory of Horticultural Plant Biology, Ministry of Education, Wuhan, China

Nan-Nan Wang

Affiliation: College of Horticulture and Forestry Sciences, Huazhong Agricultural University, Key Laboratory of Horticultural Plant Biology, Ministry of Education, Wuhan, China

Li-Na Fu

Affiliation: College of Horticulture and Forestry Sciences, Huazhong Agricultural University, Key Laboratory of Horticultural Plant Biology, Ministry of Education, Wuhan, China

Xiao Liu

Affiliation: College of Horticulture and Forestry Sciences, Huazhong Agricultural University, Key Laboratory of Horticultural Plant Biology, Ministry of Education, Wuhan, China

Xiao-Mei Hu

Affiliation: College of Horticulture and Forestry Sciences, Huazhong Agricultural University, Key Laboratory of Horticultural Plant Biology, Ministry of Education, Wuhan, China

Ting-Shuai Yan

Affiliation: College of Horticulture and Forestry Sciences, Huazhong Agricultural University, Key Laboratory of Horticultural Plant Biology, Ministry of Education, Wuhan, China

Shu-Ang Peng

* E-mail: [email protected] (CQY); [email protected] (SAP)

Affiliation: College of Horticulture and Forestry Sciences, Huazhong Agricultural University, Key Laboratory of Horticultural Plant Biology, Ministry of Education, Wuhan, China

Introduction

B is essential for higher plants [1], and soluble B content in soil and irrigation water is a crucial for determining the yield of agricultural products [2]. B deficiency (BD) has been reported in 132 plants in more than 80 countries [3]. It is well known that B deficiency will lead to different phenotypes in vascular plants, for example, the male sterility in wheat and rice [4], [5], the reduction of root biomass and nodules in legume [6], [7], the deformation of young leaves and the hypertrophy of petioles in castor bean [8], and the cracked bend veins of leaves in mulberry [9].

Citrus is an agriculturally and economically important fruit tree in the world. It is very sensitive to low B content in the soil. The occurrence of B deficiency has been reported in the major citrus producing countries, such as Spain, America, Brazil and China [3]. In China, large south and east areas contain extremely low level of soluble B (<0.25 mg kg−1) [10], [11]. Being a predominant region of naval orange production, Ganzhou an area of Jiangxi province in south China plays an important role in the country's citrus production. However, corky split veins of leaves usually occurred in main local cultivars, ‘Newhall’ navel orange (Citrus sinensis Osb.) [12], which declines the vigour of tree rapidly after the fruit set, and eventually affects fruit yield and quality in the coming years [13].

The primary function of B in higher plants is to form borate esters with apiose residues of rhamnogalacturonan II (RG-II) [14]. The formation of the complex is essential for cell wall structure and function [15] since it contributes significantly to the control of cell wall porosity [16] and tensile strength [17]. However, this participation does not seem to explain physiological phenomena of B deficiency documented in many other studies, for example, B deficiency increases the membrane permeability [18], [19] and the activity of proton-pumping ATPase in sunflower [20], causes an accumulation of phenolics through the stimulation of the enzyme phenylalanine-ammonium lyase (PAL) in tobacco and squash [21], [22] and alters amino acid profiles in white lupin [23]. Molecular genetic research have identified two types of B transporters: NIPs (nodulin-26-like intrinsic proteins), boric acid channels, and BORs, B exporters, which are both important for the uptake of B by roots, xylem loading and B distribution among leaves under B deficiency condition [2], [24], [25], [26], [27], [28], [29], [30], [31], [32]. In addition, other B-related proteins, for example, WRKY6, a low-B-induced transcription factor, is reported to be essential for Arabidopsis normal root growth under low-B condition [33], and the pathogenesis-related (PR) proteins from the PR-10 family are reported to be highly induced in low-B nodules during the legume-rhizobia interaction [34]. Recently, a quantitative trait locus (QTL) analysis for seed yield and yield-related traits under low and normal B conditions have been carried out in Brassica napus [35]. Despite the increasing studies on B deficiency in plants, little have been documented on corky split vein and the effect of the stress on leaf vein.

The purpose of this study was to get insight into the molecular mechanisms of corky split vein of ‘Newhall’ navel orange which was sandy cultured in a greenhouse under the ccondition of B deficiency. A genome-wide analysis of gene expression profiling at four phases during a long-term B deficiency treatment was performed in the experiment. Our results yielded numbers of DEGs related to the formation of corky split veins.

Materials and Methods

Plant materials and B treatments

‘Newhall’ navel orange (C. sinensis Osb.), grafted on Trifoliate orange (Poncirus trifoliata (L.) Raf.) was used in the experiment. The virus-free plants were obtained from the National Indoor Conservation Center of Virus-free Germplasms of Fruit Crops at Huazhong Agricultural University in Wuhan. According to the published method [36], [37], [38], the one- year plants grown in a greenhouse were washed with deionized water to remove surface contaminants, then were transplanted to the black pots containing B-free medium composed of quartz sand: perlite (1∶1, v/v). The plants were irrigated every other day with a modified Hoagland's No.2 nutrient solution containing either 0 (B deficiency, BD) or 0.25 mg l−1 (control, CK) B [38], [39]. The treatments started at the beginning of April, 2010 and terminated in May, 2011 when the visible symptom (corky split vein) of B deficiency appeared. Based on the morphological character of corky split vein development, both CK and BD veins collected at 28th March, 2011 (CK1 and BD1), 7th April, 2011 (CK2 and BD2), 16th April, 2011 (CK3 and BD3) and 2nd May, 2011 (CK4 and BD4) were used for the anatomical observation, DGE analysis, and quantitative real-time PCR (qRT-PCR). In addition, leaves collected on 2nd May were analyzed for the total B concentration.

Total B concentration measurement

Total B concentration was measured by the method of previous research [38]. 0.50 g of each sample was dry-ashed in a muffle furnace at 500°C for 6 h, followed by dissolution in 0.1 N HCl, and B was determined by ICP-AES.

Light microscopy for cross-section observation of veins

Microscopy observation was performed on fresh lateral veins (LVs) according the method previously described [40]. The 3rd LV with the lenghth of 0.5 centmeter (cm) near the petiole were collected (Figure S1, asterisked). LVs were immerged in the solution of formalin-acetic acid-alcohol (FAA) and vacuumized for 30 min. After fixation, the samples were routinely dehydrated in an ethanol series, and in the end they were transparentized in dimethylbenzene solution, embedded in paraffin, and sectioned at 6 µm by using microtome (Lecia RM-2255, Germany). The sections were transferred onto Superfrost Plus slides and dried for five days at 40°C–45°C, and then dyed by safranine-fast green.

Digital gene expression profiling

The B treatment was terminated when the visible symptom appeared as described above. Leaf samples of four phases covering the whole gradient of variation in morphology and anatomy of corky split veins were collected for the following RNA isolation. The total RNA were isolated from veins (cutted as Figure S1, whole) of the four stages, and were named as libraries BD1, CK1, BD2, CK2, BD3, CK3, BD4, CK4 respectively. RNA extraction was performed according to the manufacturer's instructions of TRIzol reagent (TaKaRa, Dalian, China).

Quality and quantity analysis of total RNA, library construction, and sequencing were carried out at Huada Genomics Co., Ltd., Shenzhen, China. After the total RNA was extracted from the samples, mRNA was enriched by using the oligo (dT) magnetic beads. By adding the fragmentation buffer, the mRNA was interrupted to short fragments (about 200 bp), and then the first strand cDNA was synthesized by random hexamer-primer using the mRNA fragments as templates. Buffer, dNTPs, RNase H and DNA polymerase I were added to synthesize the second strand. The double strands cDNA was purified with QiaQuick PCR extraction kit and washed with EB buffer for ending repair and the addition of single nucleotide A (adenine). Finally sequencing adaptors were ligated to the fragments. The fragments were purified by agrose gel electrophoresis and enriched by PCR amplification. The library products were ready for sequencing analysis via Illumina HiSeq™ 2000. Then the raw reads data available as fasta files were deposited in the NCBI Sequence Read Archive (SRA) database (http://www.ncbi.nlm.nih.gov/Traces/sra_sub/sub.cgi?view=submissions) under eight corresponding run accession numbers, SRR824732 (CK1), SRR824755 (CK2), SRR825214 (CK3), SRR825191 (CK4), SRR825192 (BD1), SRR825193 (BD2), SRR825194 (BD3), and SRR825195 (BD4).

To get the clean reads, the dirty raw reads were removed from raw data before data analysis by filtering reads with adaptors, reads in which unknown bases were more than 10%, and low quality reads in which the percentage of the low quality bases of quality value ≤5 were more than 50%. Subsequently, the proportion of clean reads in raw reads of the eight libraries was classified.

The clean reads of eight libraries were mapped to reference sequences [41], [42], [43] (http://www.phytozome.net/search.php?method=Org_Csinensis) using SOAPaligner/soap2 [44] and the NCBI database (http://www.ncbi.nlm.nih.gov/). Mismatches no more than 2 bases were allowed in the alignment. The reads mapped to reference sequences from multiple genes were filtered. Subsequently, sequencing saturation analysis and the randomness assessments were carried out.

Identification of differentially expressed genes

The gene expression level was calculated by using RPKM [45] method (Reads Per kb per Million reads), and the formula is shown as follows:C represents the number of reads that uniquely aligned to a gene, N represents the total number of reads that uniquely aligned to all genes, and L represents the number of bases on a gene.

All libraries of clean reads were normalized to RPKM value to obtain normalized gene expression level. A strict algorithm [46] was performed to identify DEGs between two samples. Relative change threshold in pairwise comparisons across four phases was performed the absolute value of log2(BD-RPKM/CK-RPKM). The threshold with a FDR (False Discovery Rate) ≤0.001 and the absolute value of log2 Ratio ≥1 was used to judge the significance of gene expression difference.

For Gene Ontology (GO) enrichment analysis and pathway enrichment analysis, all DEGs were mapped to GO terms in the database (http://www.geneontology.org/) and pathway terms in KEGG database (http://www.genome.jp/kegg/), respectively.

Quantitative real-time PCR analysis

To verify the DGE results, qRT-PCR was used for the investigation of DEGs expression levels. The qRT-PCR analysis protocol was performed as previously described [28]. The experiment was conducted with three replicates for each sample. Citrus actin gene was used as a normalizer and the relative expression levels of genes were presented by 2−ΔΔCT (gene/actin). Primers of actin and DEGs for qRT-PCR were shown in Table S1.

Results

B deficiency caused leaf vein tissue hypertrophy

Corky split veins were observed in mature leaves of the B-deficient plants (Figure 1D, E) at approximately one year after treatments, and the result of determination on B concentrations of leaves showed that B level was extremely lower in BD treatment (8.86 mg kg−1DW) than that in CK treatment (125.98 mg kg−1DW) (Figure 2).

[Figure omitted. See PDF.]

Figure 1. Defining developmental and anatomical characters in mature leaves.

[A–E] Morphological performance of CK leaf [A], BD1 leaf [B], BD2 leaf [C], BD3 leaf [D], and BD4 leaf [E]. [F–J] Light micrograph of lateral vein of CK [F], BD1 [G], BD2 [H], BD3 [I], and BD4 [J].

https://doi.org/10.1371/journal.pone.0065737.g001

[Figure omitted. See PDF.]

Figure 2. B concentration of whole leaf.

Columns and bars represented the means and standard error of B concentration respectively (y-axis). Three replicates were used for leaf B concentration measure. DW, dry weight.

https://doi.org/10.1371/journal.pone.0065737.g002

Morphological (Figure 1A–E) and anatomical (Figure 1F–J) characters were performed at the four stages in corky split leaf vein development to further create inventories of gene expression. As no obvious difference in morphology and anatomy among the CK veins of four stages (CK1–CK4) was observed (data not shown), thus, only one morphological or anatomical photograph was chosen to be shown as CK (Figure 1A, F). Correspondingly, the BD veins of the four stages which covered the whole gradient of variation in morphology and anatomy of corky split veins were selected for the following transcriptomic analysis. These stages included a phase (BD1, Figure 1B) of being similar with CK (Figure 1A), a phase (BD2, Figure 1C) of veins protruding from the blade, a phase (BD3, Figure 1D) of slight corking veins and a phase (BD4, Figure 1E) of seriously corky split veins.

Light micrographs of lateral veins exhibited a continuous variation of hypertrophy in vascular tissue compared with CK (Figure 1F–J). At the BD3 phase, the epidermis was destroyed by the proliferating vascular cylinder for the first time in terms of anatomical observation (Figure 1I).

Analysis of DGE libraries and reads mapping

Being operated in ice, veins were separated from leaves (Figure S1, whole). Illumina HiSeq™ 2000 platform was used to perform high throughput sequences analysis on the eight citrus leaf vein libraries to investigate the transcriptomic response to B deficiency in the vein. More than 7.1 million total reads per library were obtained and above 99.21% of total reads were identified as clean reads (Figure S2) before mapping them to the reference database (Table 1).

[Figure omitted. See PDF.]

Table 1. Major characteristics of eight libraries.

https://doi.org/10.1371/journal.pone.0065737.t001

The total clean reads of eight libraries were mapped to a reference gene database of C. sinensis, in which 46147 unigenes were included. The randomness assessments showed that the reads in every position of reference gene distributed evenly and demonstrated highly similar tendencies in the eight libraries (Figure S3) which indicated that the randomness of RNA fragmentation was sufficient for subsequent bioinformatics analysis of gene expression. Over 5.3 million clean reads per library were successfully mapped to the reference database, and more than 22897 mapped genes per library were simultaneously obtained. In addition, 55.53%–58.07% of the clean reads were mapped to the unigene database perfectly, 18.73%–19.45% of the clean reads were mapped to the unigene database with ≤2 bp mismatch and 28.49%–31.89% of the clean reads were mapped to the unigene database with a unique match (Table 1).

To estimate whether the sequencing depth was sufficient for the transcriptome coverage, the sequencing saturation was analyzed in the eight libraries (Figure S4). With the number of reads increasing, the number of detected genes was increasing in the eight libraries. However, when the number of reads reached 3 million reads or higher, the growth rate of detected genes became flatten, which showed that the number of detected genes tended to be saturated.

Changes in global gene expression in B-deficient leaf veins

Global gene expression of each sample of the four B-deficient phases (BD1–BD4) was assayed using the corresponding control phase (CK1/CK2/CK3/CK4) as a common reference. The rigorous algorithm method and the relative change threshold described above were applied to judge the significant difference of gene expression. The results showed that 7202 genes had the threshold in at least one of the pairwise comparisons, which were declared to be differentially expressed during the B deficiency course (Table 2). 1387 genes were differentially expressed between the BD1 and CK1 veins. Among these genes, 683 (49.24%) were up-regulated and 704 (50.76%) were down-regulated in response to B deficiency. Most of genes were up-regulated at the succedent phases, especially at the BD3 stage. 1986 (66.89%), 1920 (62.95), and 1662 (60.57%) genes were up-regulated at the BD3 phase compared with the CK3, BD1, and BD2 phase, respectively (Table 2).

[Figure omitted. See PDF.]

Table 2. Differentially expressed genes across all libraries.

https://doi.org/10.1371/journal.pone.0065737.t002

DEGs spanned over a wide variety of processes (pairwise comparison of BD3 and CK3 for example, Figure 3) and in functional groups (Table 3). All genes with altered expression in citrus vein at the four phases were mapped to GO terms (p-value≤0.05) in the database and were classified into several categories. The DEGs were involved in numbers of processes such as metabolic, cellular, primary metabolic, cellular metabolic and macromolecular metabolic processes (Figure 3). Some metabolic and signal transduction pathways were identified through the pathway enrichment analysis by comparing them to the whole genome background (Figure 4A). Comparing with BDs and CKs, an obvious decrease in the expression of genes was associated with photosynthesis (Table 3, Figure 4G), which is often a universal phenomenon of B deficiency treatment in plants. In addition, genes associated with the cell cycle (Table 3, Figure 4B), DNA replication (Table 3, Figure 4C), lignin biosynthesis (Table 3, Figure 4D), cytokinin signal transduction (Table 3, Figure 4E) and vascular development (Table 3, Figure 4F) were up-regulated in early or later period.

[Figure omitted. See PDF.]

Figure 3. Gene classification based on gene ontology (GO) for differentially expressed genes in BD3 and CK3 libraries of citrus vein.

Only the biological processes were used for GO analysis. The y-axis and x-axis indicated the names of clusters and the ratio of each cluster, respectively.

https://doi.org/10.1371/journal.pone.0065737.g003

[Figure omitted. See PDF.]

Figure 4. Overview of major pathways and their expression profiles during the corky split vein development.

[A] Overview of major pathways in corky split vein caused by B deficiency. [B–G] Expression profiles of the DEGs involved in cell cycle [B], DNA replication [C] lignin biosynthesis [D], cytokinin signal transduction [E], vascular development [F], and photosynthesis [G] pathways from BD1 phase to BD4 phase. Pathways with transcripts that were up-regulated or down-regulated were indicated with + or −. CTK, cytokinin. DEGs, differentially expressed genes.

https://doi.org/10.1371/journal.pone.0065737.g004

[Figure omitted. See PDF.]

Table 3. Selected genes with altered expression across all libraries.

https://doi.org/10.1371/journal.pone.0065737.t003

Genes associated with major functional group were affected in corky split veins

As the global gene expression results showed that B deficiency could affect cell cycle, DNA replication, cytoskeleton, cytokinin signal transduction, vascular development, lignin biosynthesis, and photosynthesis in citrus vein (Table 3). Genes associated with cytokinin signal transduction pathway were affected in corky split veins (Figure 4E). The cytokinin receptor authentic histidine kinase gene (HK/WOL) and type-B authentic response regulator (ARR) family member ARR12 accumulated to high levels at the BD1 and BD2 phases (Figure 4E). In addition, genes related to the three groups which were involved in cell division such as cell cycle (Figure 4B), DNA replication (Figure 4C) and cytoskeleton (data not shown in Figure 4, see Table 3) were up-regulated in the corky split veins of citrus, and these genes included CYCs, CDKs, MCMs, POLD2, PCNA2, and TUs.

Furthermore, genes involved in vascular development such as WOX4 and VND7 (Figure 4F) increased apparently in the corky split vein development. The dynamic accumulation of lignin biosynthesis genes (Figure 4D) during cellular differentiation also increased. These genes including MYB85, MYB63, MYB42, PAL, C4H, CCR, CAD, 4CL, CCoAOMT, and PER, which encode transcription factors or enzymes involved in lignin biosynthesis pathway, were up-regulated at the BD3 and BD4 phases. Moreover, several genes, PSAs, PSBs, and FNR2 involved in photosynthesis (Figure 4G), were remarkably down-regulated by B deficient treatment.

Differentially expressed genes were confirmed by qRT-PCR

To confirm the results obtained through DGE analysis, a total of nine DEGs were selected to analyze their expression profiles by qRT-PCR in the four stages. The representative genes selected for the analysis were those involved in DNA replication, as well as vascular development, lignin biosynthesis and photosynthesis pathways. The relative expression levels of CKs and BDs were compared with those of DGE data respectively (Figure 5). Despite some quantitative differences at expression level, qRT-PCR results revealed the same expression tendency as DGE results. Eight genes (MCM3, MCM5, POLD2, WOX4, MYB85, CAD, CCR, and PETE1) expressed consistently both in qRT-PCR and DGE data (Fig. 5A–D, F–I) and only one gene (VND7) by using qRT-PCR did not agree well with the DGE analysis pattern (Fig. 5E).

[Figure omitted. See PDF.]

Figure 5. qRT-PCR confirmation of differentially expressed genes from digital gene expression analysis.

[A–I] The transcript levels from qRT-PCR and digital gene expression analysis of differentially expressed genes: MCM3 [A], MCM5 [B], POLD2 [C], WOX4 [D], VND7 [E], MYB85 [F], CAD [G], CCR [H], and PETE1 [I]. Relative expression levels were calculated by qRT-PCR with actin as a standard. Columns and bars represented the means and standard error of relative expression levels (n = 3) respectively (Y axis at left). Lines represented the means of RPKM value calculated by digital gene expression analysis (Y axis at right). The x-axis indicated the phases of vein development. RPKM, Reads Per kb per Million reads.

https://doi.org/10.1371/journal.pone.0065737.g005

Discussion

In this study a symptom of corky split vein was observed in B-deficient citrus. In addition, the surveys on morphology, anatomy, and transcriptome were performed at the four stages of corky split veins to compare to the control. The histological micrographs indicated that vascular hypertrophy was related to the symptom occurrence. Through the use of high-throughput sequencing, we mapped in details the transcriptional changes that occurred during hypertrophic development. Genes associated with plant cytokinin signal transduction, cell division, vascular development were affected in corky split veins.

Corky split vein caused by B deficiency was the result of vascular hypertrophy

The deficiency of B in citrus leaf was defined on the condition that B concentration was no more than 20 mg kg−1DW [47]. B concentrations in mature leaves were determined (Figure 2). In the views of previous research, the situation of 8.86 mg kg−1DW measured in BD leaf belonged to the level of B deficiency. Meanwhile, the symptom of corky split vein (Figure 1D, E) was similar to vein of B-deficient mulberry leaves [9]. Based on their morphologic characters, the four stages were selected for the following analysis and named BD1, BD2, BD3 and BD4, respectively. Vein in the BD1 stage (Figure 1B) was similar with CK veins (Figure 1A). Vein at the BD2 phase (Figure 1C) slightly protruded from the blade compared with CK vein (Figure 1A). Corking vein was clearly observed in the BD3 stage and it protruded from the blade prominently (Figure 1D). In the BD4 stage, the leaf vein exhibited seriously corky split (Figure 1E). In addition, micrographs of light microscopy on LVs were observed through the four stages (Figure 1F–J). Compared with CK (Figure 1F), the area of vascular tissues at four phases (Figure 1G–J) became larger with treating time extending. Interestingly, due to the vascular hypertrophy, epidermis was broken at the BD3 phase and then expanded extremely at the BD4 phase. Thus cracked vein could be found through morphological detection when the epidermis was destroyed. Therefore, these dynamic observations suggested that corky split vein caused by B deficiency was the result of vascular hypertrophy.

Global gene expression was changed during the corky split vein development

DGE profiling could facilitate the identification of systemic gene expression. It reveals quantitative changes in transcript abundance on a genome-wide scale. Here, 8 transcriptomic profiles of vein were performed to identify differentially expressed genes in the four stages of corky split vein. A sequencing depth of more than 7.1 million reads per library (Table 1) was reached, and highly similar tendencies of randomness assessments (Figure S3) were analyzed in eight libraries, suggesting that the sequencing experiment was successful. To our konwledge, it is the first time that changes in transcriptomic profiles are shown in corky split vein when the naval oranges were subjected to long-term B deficiency.

Results of DGE analyses showed that we could effectively identify DEGs across a wide range of transcript abundances. The global gene transcription in the stressed leaf vein significantly altered 418–3050 genes assayed during the four developmental stages of corky split vein (Table 2). Many of the genes are known to be responding stress and stimuli (Figure 3). Decreased expression of photosynthesis was observed in this study (Table 3, Figure 4G), which is similar to the results of previous physiological experiments on other B-deficient plants [48], [49].

Cell division and vascular development were involved in the formation of corky split vein

Vascular hypertrophy caused the formation of the corky split vein by the intensifying of cell division and cell differentiation. Numerous studies on plants showed that the plant cell cycle is governed by cyclin-dependent kinases (CDKs) that are associated with their activator proteins called cyclins (CYCs), and the activity of CYC–CDK is modulated at both transcriptional and post-translational levels at the G1-to-S and G2-to-M transitions [50], [51]. CDKB2, one of the plant-specific B-type CDKs (CDKBs) genes, its expression was specific to the G2 and M phases [52], [53]. A recent study demonstrated that CDKB2;2 was required both for normal cell cycle progression and for meristem organization in Arabidopsis thaliana [54]. Genes involved in cell cycle in our experiment were differentially expressed (Table 3). The expressions of cell cycle genes were up-regulated at the BD2, BD3 and BD4 stages (Figure 4B). Thus up-regulated expressions of cell cycle genes may contribute to the proliferating vascular cylinder.

WOX4, a member of the WUSCHEL-related HOMEOBOX (WOX) gene family, is required for promoting the proliferation of procambial/cambial stem cells in A. thaliana. Compared with Col-0, the mutant wox4-1 displays less procambial cell number and smaller stele width [55]. The function of WOX4 on maintaining procambial/cambial stem cell activity has been also confirmed in Vitis vinifera L. [56] and Solanum lycopersicum [57]. VND7, a member of Vascular-Related NAC Domain gene family, acts as a positive regulator for protoxylem vessel formation in A. thaliana [58], [59]. In vascular vessels, VND7 controlls both secondary wall development and programmed cell death (PCD) of the vessels in both root and shoot tissues [60], [61]. Over-expression of VND7 could induce transdifferentiation of various cells into protoxylem-like vessel elements in both A. thaliana and Populus tomentosa. In addition, PCD-related gene transcripts are highly accumulated in VND7 over-expressing plants [59]. In this study, expressions of WOX4 and VND7, were highly increased along the development of corky split vein (Figure 4F), implying that the formation of corky split vein under B-deficient condition was related to the expression alteration of the two vascular development genes (Figure 4A).

Cytokinin performs a role in initiating the vascular hypertrophy

In previous physiological experiments, the plant phytohormone cytokinin promotes the differentiation of tracheary elements [62], [63] and contributes to the maintenance and proliferation of cambial cells [64], [65]. It is generally believed that cytokinin executes these physiological activities by regulating cell division and cell differentiation through cytokinin signal transduction pathway [66], [67], [68], [69]. In A. thaliana, WOL, one of primary cytokinin receptors, is necessary for early procambial cell division in embryogenesis [70]. The wol mutant exhibits a reduced cell number and a cell division process failing to take place in the root and lower hypocotyl region soon after the torpedo stage [71]. Type-B ARRs are implicated in the cytokinin signal transduction pathway and acts downstream of the cytokinin receptors [66], [67], [68], [69]. As transcriptional activators, the phosphorylated type-B ARRs could induce rapid transcription of cytokinin-associated target genes [72]. Elegant genetic studies in A. thaliana showed that ARR12, one of type-B ARRs, redundantly played an important (or essential) role in cytokinin signal transduction [73], [74]. An arr10 arr12 double mutant exhibited the inhibition of root elongation [73] and an arr1 arr10 arr12 triple mutant displayed stunted growth and abnormality in vascular development [74]. Our results showed that the two genes, WOL and ARR12, were up-regulated at the BD1 phase in corky split veins (Figure 4E) and these indicated that the cytokinin signal transduction pathway may play a role in initiating the vascular hypertrophy/observed phenotype.

Conclusions

In summary, our research described changes in the transcriptome of leaf veins during four stages of exposure to B deficiency using DGE analysis. In this study a symptom of corky split vein was observed in B-deficient naval orange. The histological micrographs indicated that vascular hypertrophy was related to the symptom occurrence. Through the use of high-throughput sequencing, we mapped in detail the transcriptional changes that occurred during hypertrophic development. The changes in the pathways were usually consistent with the observed symptoms. Genes associated with plant cytokinin signal transduction, cell division, vascular development were affected in corky split veins. The cytokinin signal transduction pathway may play a role in initiating the observed phenotype. This research provides a new insights into the molecular mechanisms of corky split vein caused by B deficiency.

Supporting Information

[Figure omitted. See PDF.]

Figure S1.

Examples of collecting samples for total RNA extraction (whole) and light microscope observation of vein (asterisked in white). Whole leaf vein was used for total RNA extraction experiment. The position of lateral vein, at the 3rd item near the petiole, for light microscope observation was marked using a white asterisk.

https://doi.org/10.1371/journal.pone.0065737.s001

(TIF)

Figure S2.

Classification of raw reads of eight libraries.

https://doi.org/10.1371/journal.pone.0065737.s002

(TIF)

Figure S3.

Randomness assessments of the eight libraries. [A-H] Randomness assessments of library CK1 [A], BD1 [B], CK2 [C], BD2 [D], CK3 [E], BD3 [F], CK4 [G] and BD4 [H].

https://doi.org/10.1371/journal.pone.0065737.s003

(TIF)

Figure S4.

Sequencing saturation analysis of the eight libraries. [A–H] Sequencing saturation analysis of library CK1 [A], BD1 [B], CK2 [C], BD2 [D], CK3 [E], BD3 [F], CK4 [G] and BD4 [H]. With the number of reads increasing, the number of detected genes was increasing in the eight libraries. However, when the number of reads reached 3 million reads or higher, the growth rate of detected genes became flatten.

https://doi.org/10.1371/journal.pone.0065737.s004

(TIF)

Table S1.

Primer sequences for qRT-PCR.

https://doi.org/10.1371/journal.pone.0065737.s005

(XLSX)

Acknowledgments

The authors thank the International Citrus Genome Consortium (ICGC) for free use of Citrus sinensis genome sequence (http://www.phytozome.net/search.php?method=Org_Csinensis).

Author Contributions

Conceived and designed the experiments: SAP CQY. Performed the experiments: CQY JCA SL LFJ. Analyzed the data: CQY LFJ GFZ QJW. Contributed reagents/materials/analysis tools: CQY NNW LNF XL XMH TSY. Wrote the paper: SAP CQY YZL HQY. Revised the manuscript: SAP CQY.

Citation: Yang C-Q, Liu Y-Z, An J-C, Li S, Jin L-F, Zhou G-F, et al. (2013) Digital Gene Expression Analysis of Corky Split Vein Caused by Boron Deficiency in ‘Newhall’ Navel Orange (Citrus sinensis Osbeck) for Selecting Differentially Expressed Genes Related to Vascular Hypertrophy. PLoS ONE 8(6): e65737. https://doi.org/10.1371/journal.pone.0065737

References

1. Warington K (1923) The effect of boric acid and borax on the broad bean and certain other plants. Annals of Botany 37: 629–672.

2. Tanaka M, Fujiwara T (2008) Physiological roles and transport mechanisms of boron: perspectives from plants. Pflügers Archiv European Journal of Physiology 456: 671–677.

3. Shorrocks VM (1997) The occurrence and correction of boron deficiency. Plant and Soil 193: 121–148.

4. Rerkasem B, Jamjod S (2004) Boron deficiency in wheat: a review. Field Crops Research 89: 173–186.

5. Garg O, Sharma A, Kona GRSS (1979) Effect of boron on the pollen vitality and yield of rice plants (Oryza sativa L. var. Jaya). Plant and Soil 52: 591–594.

6. Alves M, Moes S, Jenö P, Pinheiro C, Passarinho J, et al. (2011) The analysis of Lupinus albus root proteome revealed cytoskeleton altered features due to long-term boron deficiency. Journal of Proteomics 74: 1351–1363.

7. Bolanos L, Brewin NJ, Bonilla I (1996) Effects of boron on rhizobium-legume cell-surface interactions and nodule development. Plant Physiology 110: 1249–1256.

8. Silva DH, Rossi ML, Boaretto AE, Nogueira NL, Muraoka T (2008) Boron affects the growth and ultrastructure of castor bean plants. Scientia Agricola 65: 659–664.

9. Tewari RK, Kumar P, Sharma PN (2010) Morphology and oxidative physiology of boron-deficient mulberry plants. Tree Physiology 30: 68–77.

10. Liu Z, Zhu QQ, Tang LH (1989) Regularities of content and distribuition of boron in soils. Acta Pedologica Sinica 4: 006.

11. Liu Z, Zhu QQ, Tong LH (1980) Boron-deficient soils and their distribution in China. Acta Pedologica Sinica 17: 228–239.

12. Jiang CC, Wang YH, Liu GD, Xia Y, Peng SA, et al. (2009) Effect of boron on the leaf etiolation and fruit drop of Newhall Navel orange in southern Jiangxi. Plant Nutrition and Fertilizer Science 15: 656–611.

13. Xiao JX, Yan X, Peng SA, Fang YW (2007) Seasonal changes of mineral nutrients in fruit and leaves of ‘Newhall’ and ‘Skagg's Bonanza’ navel oranges. Journal of Plant Nutrition 30: 671–690.

14. Kobayashi M, Matoh T, Azuma J (1996) Two chains of rhamnogalacturonan II are cross-linked by borate-diol ester bonds in higher plant cell walls. Plant Physiology 110: 1017–1020.

15. O'Neill MA, Ishii T, Albersheim P, Darvill AG (2004) Rhamnogalacturonan II: structure and function of a borate cross-linked cell wall pectic polysaccharide. Plant Biology 55: 109–139.

16. Fleischer A, O'Neill MA, Ehwald R (1999) The pore size of non-graminaceous plant cell walls is rapidly decreased by borate ester cross-linking of the pectic polysaccharide rhamnogalacturonan II. Plant Physiology 121: 829–838.

17. Ryden P, Sugimoto-Shirasu K, Smith AC, Findlay K, Reiter WD, et al. (2003) Tensile properties of Arabidopsis cell walls depend on both a xyloglucan cross-linked microfibrillar network and rhamnogalacturonan II-borate complexes. Plant Physiology 132: 1033–1040.

18. Brown PH, Bellaloui N, Wimmer MA, Bassil ES, Ruiz J, et al. (2002) Boron in plant biology. Plant Biology 4: 205–223.

19. Cakmak I, Kurz H, Marschner H (1995) Short-term effects of boron, germanium and high light intensity on membrane permeability in boron deficient leaves of sunflower. Physiologia Plantarum 95: 11–18.

20. Blaser-Grill J, Knoppik D, Amberger A, Goldbach H (1989) Influence of boron on the membrane potential in Elodea densa and Helianthus annuus roots and H+ extrusion of suspension cultured Daucus carota cells. Plant Physiology 90: 280–284.

21. Ruiz JM, Bretones G, Baghour M, Ragala L, Belakbir A, et al. (1998) Relationship between boron and phenolic metabolism in tobacco leaves. Phytochemistry 48: 269–272.

22. Cara FA, Sánchez E, Ruiz JM, Romero L (2002) Is phenol oxidation responsible for the short-term effects of boron deficiency on plasma-membrane permeability and function in squash roots? Plant Physiology and Biochemistry 40: 853–858.

23. Alves M, Chicau P, Matias H, Passarinho J, Pinheiro C, et al. (2011) Metabolic analysis revealed altered amino acid profiles in Lupinus albus organs as a result of boron deficiency. Physiologia Plantarum 142: 224–232.

24. Takano J, Wada M, Ludewig U, Schaaf G, Von Wirén N, et al. (2006) The Arabidopsis major intrinsic protein NIP5; 1 is essential for efficient boron uptake and plant development under boron limitation. The Plant Cell 18: 1498–1509.

25. Miwa K, Tanaka M, Kamiya T, Fujiwara T (2010) Molecular mechanisms of boron transport in plants: Involvement of Arabidopsis NIP5; 1 and NIP6; 1. Mips and their Role in the Exchange of Metalloids 83–96.

26. Tanaka M, Takano J, Chiba Y, Lombardo F, Ogasawara Y, et al. (2011) Boron-dependent degradation of NIP5; 1 mRNA for acclimation to excess boron conditions in Arabidopsis. The Plant Cell 23: 3547–3559.

27. Tanaka M, Wallace IS, Takano J, Roberts DM, Fujiwara T (2008) NIP6; 1 is a boric acid channel for preferential transport of boron to growing shoot tissues in Arabidopsis. The Plant Cell 20: 2860–2875.

28. An JC, Liu YZ, Yang CQ, Zhou GF, Wei QJ, et al. (2012) Isolation and expression analysis of CiNIP5 , a citrus boron transport gene involved in tolerance to boron deficiency. Scientia Horticulturae 142: 149–154.

29. Takano J, Miwa K, Yuan L, von Wirén N, Fujiwara T (2005) Endocytosis and degradation of BOR1, a boron transporter of Arabidopsis thaliana, regulated by boron availability. PNAS 102: 12276–12281.

30. Miwa K, Takano J, Fujiwara T (2006) Improvement of seed yields under boron-limiting conditions through overexpression of BOR1, a boron transporter for xylem loading, in Arabidopsis thaliana. The Plant Journal 46: 1084–1091.

31. Takano J, Yamagami M, Noguchi K, Hayashi H, Fujiwara T (2001) Preferential translocation of boron to young leaves in Arabidopsis thaliana regulated by the BOR1 gene. Soil Science and Plant Nutrition 47: 345–357.

32. Nakagawa Y, Hanaoka H, Kobayashi M, Miyoshi K, Miwa K, et al. (2007) Cell-type specificity of the expression of Os BOR1, a rice efflux boron transporter gene, is regulated in response to boron availability for efficient boron uptake and xylem loading. The Plant Cell 19: 2624–2635.

33. Kasajima I, Ide Y, Yokota Hirai M, Fujiwara T (2010) WRKY6 is involved in the response to boron deficiency in Arabidopsis thaliana. Physiologia Plantarum 139: 80–92.

34. Reguera M, Bonilla I, Bolaños L (2010) Boron deficiency results in induction of pathogenesis-related proteins from the PR-10 family during the legume–rhizobia interaction. Journal of plant physiology 167: 625–632.

35. Zhao ZK, Wu LK, Nian FZ, Ding GD, Shi TX, et al. (2012) Dissecting quantitative trait loci for boron efficiency across multiple environments in Brassica napus. PloS ONE 7: e45215.

36. Papadakis I, Dimassi K, Bosabalidis A, Therios I, Patakas A, et al. (2004) Effects of B excess on some physiological and anatomical parameters of ‘Navelina’ orange plants grafted on two rootstocks. Environmental and Experimental Botany 51: 247–257.

37. Papadakis IE, Dimassi KN, Bosabalidis AM, Therios IN, Patakas A, et al. (2004) Boron toxicity in ‘Clementine’ mandarin plants grafted on two rootstocks. Plant Science 166: 539–547.

38. Sheng O, Song SW, Peng SA, Deng XX (2009) The effects of low boron on growth, gas exchange, boron concentration and distribution of ‘Newhall’ navel orange (Citrus sinensis Osb.) plants grafted on two rootstocks. Scientia Horticulturae 121: 278–283.

39. Hoagland DR, Arnon DI (1950) The water-culture method for growing plants without soil. Circular California Agricultural Experiment Station 347.

40. Li HP (2009). Plant microscopy technique (in Chinese). Beijing, China: Science Press.

41. Ambion website. Available: http://www.phytozome.net/searchphp?method=Org_Csinensis. Accessed 2011 Jan 15.

42. Gmitter FG (2012) New Citrus Genome Sequence Resources; 2012. Plant and Animal Genome

43. Main D The Citrus Genome Database: A Community Resource for Citrus Research; 2012. Plant and Animal Genome

44. Li RQ, Yu C, Li YR, Lam TW, Yiu SM, et al. (2009) SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics 25: 1966–1967.

45. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B (2008) Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nature Methods 5: 621–628.

46. Audic S, Claverie JM (1997) The significance of digital gene expression profiles. Genome Research 7: 986–995.

47. Obreza TA, Morgan K (2008) Nutrition of Florida citrus trees. Florida: University of Florida.Press. 27 p.

48. Kastori R, Plesničar M, Panković D, Sakač Z (1995) Photosynthesis, chlorophyll fluorescence and soluble carbohydrates in sunflower leaves as affected by boron deficiency. Journal of Plant Nutrition 18: 1751–1763.

49. Han S, Chen LS, Jiang HX, Smith BR, Yang LT, et al. (2008) Boron deficiency decreases growth and photosynthesis, and increases starch and hexoses in leaves of citrus seedlings. Journal of Plant Physiology 165: 1331–1341.

50. Inzé D, De Veylder L (2006) Cell cycle regulation in plant development1. Annual Review of Genetics 40: 77–105.

51. Komaki S, Sugimoto K (2012) Control of the plant cell cycle by developmental and environmental cues. Plant and Cell Physiology 53: 953–964.

52. Magyar Z, Mészáros T, Miskolczi P, Deák M, Fehér A, et al. (1997) Cell cycle phase specificity of putative cyclin-dependent kinase variants in synchronized alfalfa cells. The Plant Cell 9: 223–235.

53. Porceddu A, Stals H, Reichheld JP, Segers G, De Veylder L, et al. (2001) A plant-specific cyclin-dependent kinase is involved in the control of G2/M progression in plants. Journal of Biological Chemistry 276: 36354–36360.

54. Andersen SU, Buechel S, Zhao Z, Ljung K, Novák O, et al. (2008) Requirement of B2-type cyclin-dependent kinases for meristem integrity in Arabidopsis thaliana. The Plant Cell 20: 88–100.

55. Hirakawa Y, Kondo Y, Fukuda H (2010) TDIF peptide signaling regulates vascular stem cell proliferation via the WOX4 homeobox gene in Arabidopsis. The Plant Cell 22: 2618–2629.

56. Dai R, Jin HP, Wang Z, Avihai P, Xu HY, et al. (2011) Cloning and characterization of WOX4 Gene from Vitis vinifera L. involved in stem cell regulation. Agricultural Sciences in China 10: 1861–1871.

57. Ji J, Strable J, Shimizu R, Koenig D, Sinha N, et al. (2010) WOX4 promotes procambial development. Plant Physiology 152: 1346–1356.

58. Kubo M, Yamaguchi M, Demura T (2004) VND7 Regulates xylem vessel formation in Arabidopsis. Journal of Plant Research 117: 121.

59. Kubo M, Udagawa M, Nishikubo N, Horiguchi G, Yamaguchi M, et al. (2005) Transcription switches for protoxylem and metaxylem vessel formation. Genes and Development 19: 1855–1860.

60. Wang HZ, Dixon RA (2012) On–off switches for secondary cell wall biosynthesis. Molecular plant 5: 297–303.

61. Yamaguchi M, Kubo M, Fukuda H, Demura T (2008) Vascular-Related NAC-Domain7 is involved in the differentiation of all types of xylem vessels in Arabidopsis roots and shoots. The Plant Journal 55: 652–664.

62. Aloni R (1987) Differentiation of vascular tissues. Annual Review of Plant Physiology 38: 179–204.

63. Fukuda H (1997) Tracheary element differentiation. The Plant Cell 9: 1147.

64. Matsumoto-Kitano M, Kusumoto T, Tarkowski P, Kinoshita-Tsujimura K, Václavíková K, et al. (2008) Cytokinins are central regulators of cambial activity. PNAS 105: 20027–20031.

65. Nieminen K, Immanen J, Laxell M, Kauppinen L, Tarkowski P, et al. (2008) Cytokinin signaling regulates cambial development in poplar. PNAS 105: 20032–20037.

66. Heyl A, Schmülling T (2003) Cytokinin signal perception and transduction. Current Opinion in Plant Biology 6: 480–488.

67. Kakimoto T (2003) Perception and signal transduction of cytokinins. Annual Review of Plant Biology 54: 605–627.

68. Ferreira FJ, Kieber JJ (2005) Cytokinin signaling. Current Opinion in Plant Biology 8: 518–525.

69. Muller B, Sheen J (2007) Advances in cytokinin signaling. Science Signalling 318: 68.

70. Mähönen AP, Bonke M, Kauppinen L, Riikonen M, Benfey PN, et al. (2000) A novel two-component hybrid molecule regulates vascular morphogenesis of the Arabidopsis root. Genes and Development 14: 2938–2943.

71. Scheres B, Di Laurenzio L, Willemsen V, Hauser MT, Janmaat K, et al. (1995) Mutations affecting the radial organisation of the Arabidopsis root display specific defects throughout the embryonic axis. Development 121: 53–62.

72. Sakai H, Aoyama T, Oka A (2000) Arabidopsis ARR1 and ARR2 response regulators operate as transcriptional activators. The Plant Journal 24: 703–711.

73. Yokoyama A, Yamashino T, Amano YI, Tajima Y, Imamura A, et al. (2007) Type-B ARR transcription factors, ARR10 and ARR12, are implicated in cytokinin-mediated regulation of protoxylem differentiation in roots of Arabidopsis thaliana. Plant and Cell Physiology 48: 84–96.

74. Ishida K, Yamashino T, Yokoyama A, Mizuno T (2008) Three type-B response regulators, ARR1, ARR10 and ARR12, play essential but redundant roles in cytokinin signal transduction throughout the life cycle of Arabidopsis thaliana. Plant and Cell Physiology 49: 47–57.