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 1^st 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 2^nd and 3^rd 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 [49][1], and soluble B content in soil
   and irrigation water is a crucial for determining the yield of
   agricultural products [50][2]. B deficiency (BD) has been reported in
   132 plants in more than 80 countries [51][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 [52][4], [53][5], the
   reduction of root biomass and nodules in legume [54][6], [55][7], the
   deformation of young leaves and the hypertrophy of petioles in castor
   bean [56][8], and the cracked bend veins of leaves in mulberry [57][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 [58][3].
   In China, large south and east areas contain extremely low level of
   soluble B (<0.25 mg kg^−1) [59][10], [60][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.)
   [61][12], which declines the vigour of tree rapidly after the fruit
   set, and eventually affects fruit yield and quality in the coming years
   [62][13].

   The primary function of B in higher plants is to form borate esters
   with apiose residues of rhamnogalacturonan II (RG-II) [63][14]. The
   formation of the complex is essential for cell wall structure and
   function [64][15] since it contributes significantly to the control of
   cell wall porosity [65][16] and tensile strength [66][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 [67][18], [68][19] and
   the activity of proton-pumping ATPase in sunflower [69][20], causes an
   accumulation of phenolics through the stimulation of the enzyme
   phenylalanine-ammonium lyase (PAL) in tobacco and squash [70][21],
   [71][22] and alters amino acid profiles in white lupin [72][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 [73][2], [74][24], [75][25], [76][26], [77][27], [78][28],
   [79][29], [80][30], [81][31], [82][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 [83][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 [84][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 [85][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 [86][36],
   [87][37], [88][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 [89][38],
   [90][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
   [91][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 [92][40]. The 3^rd LV with
   the lenghth of 0.5 centmeter (cm) near the petiole were collected
   ([93]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 [94]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
   ([95]http://www.ncbi.nlm.nih.gov/Traces/sra_sub/sub.cgi?view=submission
   s) 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
   [96][41], [97][42], [98][43]
   ([99]http://www.phytozome.net/search.php?method=Org_Csinensis) using
   SOAPaligner/soap2 [100][44] and the NCBI database
   ([101]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 [102][45] method
   (Reads Per kb per Million reads), and the formula is shown as follows:
   graphic file with name pone.0065737.e001.jpg

   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 [103][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
   ([104]http://www.geneontology.org/) and pathway terms in KEGG database
   ([105]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 [106][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
   [107]Table S1.

Results

B deficiency caused leaf vein tissue hypertrophy

   Corky split veins were observed in mature leaves of the B-deficient
   plants ([108]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) ([109]Figure 2).

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

   Figure 1
   [110]Open in a new tab

   [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].

Figure 2. B concentration of whole leaf.

   Figure 2
   [111]Open in a new tab

   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.

   Morphological ([112]Figure 1A–E) and anatomical ([113]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
   ([114]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, [115]Figure 1B) of being
   similar with CK ([116]Figure 1A), a phase (BD2, [117]Figure 1C) of
   veins protruding from the blade, a phase (BD3, [118]Figure 1D) of
   slight corking veins and a phase (BD4, [119]Figure 1E) of seriously
   corky split veins.

   Light micrographs of lateral veins exhibited a continuous variation of
   hypertrophy in vascular tissue compared with CK ([120]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
   ([121]Figure 1I).

Analysis of DGE libraries and reads mapping

   Being operated in ice, veins were separated from leaves ([122]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 ([123]Figure S2)
   before mapping them to the reference database ([124]Table 1).

Table 1. Major characteristics of eight libraries.

   Summary CK1 BD1 CK2 BD2 CK3 BD3 CK4 BD4
   Raw reads TN 7179660 7623813 7259087 7178343 7607537 7596798 7579221
   7143559
   Clean reads TN 7132497 7565451 7209850 7132317 7554064 7536975 7528213
   7091768
   TP 99.34% 99.23% 99.32% 99.36% 99.30% 99.21% 99.33% 99.27%
   Mapped reads TN 5415328 5834339 5488209 5488088 5720491 5618233 5609421
   5356143
   TP 75.92% 77.12% 76.12% 76.95% 75.73% 74.54% 74.51% 75.53%
   Mapped genes TN 23290 22897 23433 23054 22374 23792 23422 23228
   TP 50.47% 49.62% 50.78% 49.96% 48.48% 51.56% 50.76% 50.33%
   Perfect match TN 4079188 4393468 4123087 4110504 4250953 4184929
   4198715 4011645
   TP 57.19% 58.07% 57.19% 57.63% 56.27% 55.53% 55.77% 56.57%
   ≤2 bp mismatch TN 1336140 1440871 1365122 1377584 1469538 1433304
   1410706 1344498
   TP 18.73% 19.05% 18.93% 19.31% 19.45% 19.02% 18.74% 18.96%
   Unique match TN 2076283 2235674 2090868 2109145 2152356 2348318 2187769
   2261917
   TP 29.11% 29.55% 29.00% 29.57% 28.49% 31.16% 29.06% 31.89%
   [125]Open in a new tab

   TN, Total Number. TP, Total Percentage.

   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 ([126]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
   ([127]Table 1).

   To estimate whether the sequencing depth was sufficient for the
   transcriptome coverage, the sequencing saturation was analyzed in the
   eight libraries ([128]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 ([129]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 ([130]Table 2).

Table 2. Differentially expressed genes across all libraries.

   CK1 vs. BD1 CK2 vs. BD2 CK3 vs. BD3 CK4 vs. BD4 BD1 vs. BD2 BD1 vs. BD3
   BD1 vs. BD4 BD2 vs. BD3 BD2 vs. BD4 BD3 vs. BD4
   Total 1387 418 2969 2168 798 3050 2350 2744 2150 710
   Up-regulated 683 173 1986 1195 348 1920 1367 1662 1232 295
   Down-regulated 704 245 983 973 450 1130 983 1082 918 415
   [131]Open in a new tab

   All the genes mapped to the reference sequence were examined for their
   expression differences across the eight libraries. CK1, CK2, CK3, CK4,
   BD1, BD2, BD3 and BD4 represent eight libraries, respectively. Numbers
   of differentially expressed genes represent across sense transcripts,
   using the threshold with a FDR (False Discovery Rate) ≤0.001 and the
   absolute value of log2 Ratio ≥1 for controlling false discovery rates.

   DEGs spanned over a wide variety of processes (pairwise comparison of
   BD3 and CK3 for example, [132]Figure 3) and in functional groups
   ([133]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 ([134]Figure
   3). Some metabolic and signal transduction pathways were identified
   through the pathway enrichment analysis by comparing them to the whole
   genome background ([135]Figure 4A). Comparing with BDs and CKs, an
   obvious decrease in the expression of genes was associated with
   photosynthesis ([136]Table 3, [137]Figure 4G), which is often a
   universal phenomenon of B deficiency treatment in plants. In addition,
   genes associated with the cell cycle ([138]Table 3, [139]Figure 4B),
   DNA replication ([140]Table 3, [141]Figure 4C), lignin biosynthesis
   ([142]Table 3, [143]Figure 4D), cytokinin signal transduction
   ([144]Table 3, [145]Figure 4E) and vascular development ([146]Table 3,
   [147]Figure 4F) were up-regulated in early or later period.

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

   [148]Figure 3
   [149]Open in a new tab

   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.

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

   Function-al group Database ID Symbol RPKM value
   BD1 BD2 BD3 BD4 CK1 CK2 CK3 CK4
   Cell cycle
   orange1.1g015387m CYCB2;3 1.3 3.6 14.8 6.7 4.7 1.4 1.8 3.1
   orange1.1g018657m CYCD1;1 4.8 3.1 21.1 15 8 7.1 3.1 2.6
   orange1.1g021982m CDKB2;2 2.8 4 17.9 8.4 3 1.5 1 3.8
   orange1.1g016078m CYCB1;4 3.3 8.3 35.2 14.5 12 6.6 3.4 12.4
   orange1.1g011434m CYCA2;1 1.7 1.8 9.3 2.6 0.7 1.9 1.1 1.1
   orange1.1g016163m CYCD3;1 10.2 20.1 132.3 79.8 33.1 15.7 13.3 27
   DNA replication
   orange1.1g002353m MCM2 7.4 19.1 37.1 16.4 12.6 13.3 5.5 12.4
   orange1.1g003637m MCM4 3.5 7.5 15 7.5 6.1 3.6 1.8 4
   orange1.1g003809m MCM6 2.8 7.3 9.4 4.1 6.9 5.4 1.6 4
   orange1.1g004502m MCM3 1.6 8.5 18.1 6.8 7.6 4.6 4.5 7
   orange1.1g004862m MCM5 3.2 9.1 19.2 8.5 7.5 6.4 3.1 4.7
   orange1.1g013632m POLD2 13.8 18.7 24.1 29.4 16.3 15.6 11.1 16.6
   orange1.1g024527m PCNA2 20.7 27.7 57 38.4 25.8 19.8 21.5 25.9
   Cytoskeleton
   orange1.1g013066m TUA2 9.4 26 80.9 33.4 12.2 12.2 5.3 16.6
   orange1.1g013291m TUB8 19.7 33.7 82.2 55.2 26.8 26.3 26.4 21.5
   orange1.1g013439m TUB6 3.9 5.6 21.7 12.6 4.8 1.8 3.2 3.7
   Cytokinin signal transduction
   orange1.1g001846m WOL 93.4 105.3 22 46.4 44.2 52.9 10.4 38.6
   orange1.1g005719m ARR12 53.4 32.5 20.9 14.8 13.8 27.8 34.9 14.3
   Vascular development
   orange1.1g027928m WOX4 67.6 57.9 304.5 218 73.4 95.1 103.7 87.2
   orange1.1g019821m VND7 1.3 - 8.4 7.9 3.3 2.8 1.4 4.1
   Lignin biosynthesis
   orange1.1g024849m MYB85 13.4 6.5 30.3 33 11.1 16.5 16.5 12.6
   orange1.1g045411m MYB63 0.4 - 23.3 34.4 - - - -
   orange1.1g024441m MYB42 0.3 - 9.6 4.1 0.8 0.7 5.4 0.4
   orange1.1g037382m PAL 5.6 4 59.7 38.7 1.4 2.9 2.8 3.8
   orange1.1g012770m C4H 0.4 0.1 14.5 13.4 0.2 0.1 - 0.1
   orange1.1g042783m 4CL2 33.4 19.8 348.9 339.6 29 31.8 38.3 22.3
   orange1.1g026418m CCR 38.9 42.9 176.1 258.9 34.8 30.9 50.4 38.2
   orange1.1g025824m CCoAOMT 90.2 59.2 432.4 377 79.2 83.6 60.8 78.7
   orange1.1g020266m CAD 29.3 12 83.2 69.4 21.3 14.5 15.6 36.6
   orange1.1g018871m PER 10.7 10.3 388.2 186.6 15.3 12.5 4 16.9
   Photosynthesis
   orange1.1g017613m ATPC1 632 844.5 230.9 371.9 942 873.3 952.6 1085.9
   orange1.1g019870m PSBO2 74.9 121.3 38.8 58.6 134.4 108.8 156 131.9
   orange1.1g024892m ATPDZ 379.2 526.9 107 180.9 718.5 522.4 364.7 651.5
   orange1.1g026780m PSBQ2 62 90.6 43.4 56.9 111.7 89.2 189.6 140.1
   orange1.1g026929m PSAF 8.1 8.3 5.6 8.3 11.5 15 20.7 12.9
   orange1.1g027843m PSAL 1690.7 2199.2 1122.7 1613.9 2453.2 2297 5119.5
   3232.8
   orange1.1g028581m PSAD2 1078.4 1490.4 452.4 661.4 1570.3 1446.6 2151.9
   1650
   orange1.1g029024m PSBY 1155.4 1273.5 428.7 601.6 1180.8 1196.8 2557.4
   1333.6
   orange1.1g030950m PETE1 700.9 1207.5 475.4 793.7 1122.8 1036.6 2660.4
   1525.1
   orange1.1g031508m PSAG 2578.3 2862.3 1040.8 1661.6 2931.6 3269.3 4827.5
   3419.4
   orange1.1g032245m PSAH2 935.4 1221.4 678 737.1 1151.6 1164.6 2804.9
   1688.6
   orange1.1g032691m PSBW 523.1 870.4 175.4 223.1 903.9 748.1 1388.9
   1196.9
   orange1.1g033067m PSAK 855.8 1312 594.4 674.2 952.4 933.6 2363.3 1288.5
   orange1.1g040346m ATPD 191.1 260.3 78.1 120.7 226 184.9 227.1 258
   orange1.1g043022m FNR2 1112.5 2032.7 444.3 420.5 2539.4 1983.5 1883.2
   2855.3
   [150]Open in a new tab

   RPKM value of which gene expression was not detected in a library was
   indicated with -. RPKM, Reads Per kb per Million reads.

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

   [151]Figure 4
   [152]Open in a new tab

   [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.

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 ([153]Table 3). Genes associated with
   cytokinin signal transduction pathway were affected in corky split
   veins ([154]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 ([155]Figure 4E). In addition, genes related to the three groups
   which were involved in cell division such as cell cycle ([156]Figure
   4B), DNA replication ([157]Figure 4C) and cytoskeleton (data not shown
   in [158]Figure 4, see [159]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 ([160]Figure 4F) increased apparently in the corky split vein
   development. The dynamic accumulation of lignin biosynthesis genes
   ([161]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 ([162]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 ([163]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 ([164]Fig. 5A–D, F–I) and only one gene (VND7) by
   using qRT-PCR did not agree well with the DGE analysis pattern
   ([165]Fig. 5E).

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

   [166]Figure 5
   [167]Open in a new tab

   [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.

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 [168][47]. B
   concentrations in mature leaves were determined ([169]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 ([170]Figure 1D, E) was similar to vein of
   B-deficient mulberry leaves [171][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
   ([172]Figure 1B) was similar with CK veins ([173]Figure 1A). Vein at
   the BD2 phase ([174]Figure 1C) slightly protruded from the blade
   compared with CK vein ([175]Figure 1A). Corking vein was clearly
   observed in the BD3 stage and it protruded from the blade prominently
   ([176]Figure 1D). In the BD4 stage, the leaf vein exhibited seriously
   corky split ([177]Figure 1E). In addition, micrographs of light
   microscopy on LVs were observed through the four stages ([178]Figure
   1F–J). Compared with CK ([179]Figure 1F), the area of vascular tissues
   at four phases ([180]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 ([181]Table 1) was reached, and highly similar tendencies
   of randomness assessments ([182]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.

   [183]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
   ([184]Table 2). Many of the genes are known to be responding stress and
   stimuli ([185]Figure 3). Decreased expression of photosynthesis was
   observed in this study ([186]Table 3, [187]Figure 4G), which is similar
   to the results of previous physiological experiments on other
   B-deficient plants [188][48], [189][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 G[1]-to-S and G[2]-to-M transitions [190][50], [191][51]. CDKB2,
   one of the plant-specific B-type CDKs (CDKBs) genes, its expression was
   specific to the G[2] and M phases [192][52], [193][53]. A recent study
   demonstrated that CDKB2;2 was required both for normal cell cycle
   progression and for meristem organization in Arabidopsis thaliana
   [194][54]. Genes involved in cell cycle in our experiment were
   differentially expressed ([195]Table 3). The expressions of cell cycle
   genes were up-regulated at the BD2, BD3 and BD4 stages ([196]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 [197][55]. The
   function of WOX4 on maintaining procambial/cambial stem cell activity
   has been also confirmed in Vitis vinifera L. [198][56] and Solanum
   lycopersicum [199][57]. VND7, a member of Vascular-Related NAC Domain
   gene family, acts as a positive regulator for protoxylem vessel
   formation in A. thaliana [200][58], [201][59]. In vascular vessels,
   VND7 controlls both secondary wall development and programmed cell
   death (PCD) of the vessels in both root and shoot tissues [202][60],
   [203][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 [204][59]. In this
   study, expressions of WOX4 and VND7, were highly increased along the
   development of corky split vein ([205]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
   ([206]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 [207][62], [208][63]
   and contributes to the maintenance and proliferation of cambial cells
   [209][64], [210][65]. It is generally believed that cytokinin executes
   these physiological activities by regulating cell division and cell
   differentiation through cytokinin signal transduction pathway
   [211][66], [212][67], [213][68], [214][69]. In A. thaliana, WOL, one of
   primary cytokinin receptors, is necessary for early procambial cell
   division in embryogenesis [215][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 [216][71].
   Type-B ARRs are implicated in the cytokinin signal transduction pathway
   and acts downstream of the cytokinin receptors [217][66], [218][67],
   [219][68], [220][69]. As transcriptional activators, the phosphorylated
   type-B ARRs could induce rapid transcription of cytokinin-associated
   target genes [221][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 [222][73], [223][74].
   An arr10 arr12 double mutant exhibited the inhibition of root
   elongation [224][73] and an arr1 arr10 arr12 triple mutant displayed
   stunted growth and abnormality in vascular development [225][74]. Our
   results showed that the two genes, WOL and ARR12, were up-regulated at
   the BD1 phase in corky split veins ([226]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.

   (TIF)
   [227]Click here for additional data file.^ (2MB, tif)
   Figure S2

   Classification of raw reads of eight libraries.

   (TIF)
   [228]Click here for additional data file.^ (227.7KB, 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].

   (TIF)
   [229]Click here for additional data file.^ (80.3KB, 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.

   (TIF)
   [230]Click here for additional data file.^ (89.5KB, tif)
   Table S1

   Primer sequences for qRT-PCR.

   (XLSX)
   [231]Click here for additional data file.^ (11.5KB, xlsx)

Acknowledgments