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

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.


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 21 ) [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 Brelated proteins, for example, WRKY6, a low-B-induced transcrip-tion 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.

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 21 (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 28 th March, 2011 (CK1 and BD1), 7 th April, 2011 (CK2 and BD2), 16 th April, 2011 (CK3 and BD3) and 2 nd 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 2 nd 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 500uC 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 3 rd 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 mm by using microtome (Lecia RM-2255, Germany). The sections were transferred onto Superfrost Plus slides and dried for five days at 40uC-45uC, 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 TM 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 log 2 (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 2DDCT (gene/actin). Primers of actin and DEGs for qRT-PCR were shown in Table  S1.

B deficiency caused leaf vein tissue hypertrophy
Corky split veins were observed in mature leaves of the Bdeficient 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 21 DW) than that in CK treatment (125.98 mg kg 21 DW) ( Figure 2).
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 TM 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).  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).
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.
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   (Fig. 5E).

Discussion
In this study a symptom of corky split vein was observed in Bdeficient 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 21 DW [47]. B concentrations in mature leaves were determined (Figure 2). In the views of previous research, the situation of 8.86 mg kg 21 DW 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   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 posttranslational levels at the G 1 -to-S and G 2 -to-M transitions [50,51]. CDKB2, one of the plant-specific B-type CDKs (CDKBs) genes, its expression was specific to the G 2 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 cytokininassociated 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 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 3 rd item near the petiole, for light microscope observation was marked using a white asterisk.