Comparative Transcriptome Profiling Indicated that Leaf Mesophyll and Leaf Vasculature have Different Drought Response Mechanisms in Cassava

Drought stress is one of the major environmental factors that limited crop’s growth and production. Cassava known as a tropical crop that is widely distributed in Sub-Saharan Africa. It has a strong drought tolerance and can grow well under tough environmental conditions. Therefore, understanding how cassava responds to drought stress and coordinates survival and accumulation has great theoretical significance for improving crop drought resistance breeding. Many studies on cassava drought responses mainly focused on the leaf and whole seedling. Nevertheless, how the vasculature plays an important role in plant response to water deficiency remains to be fully elucidated. Here, comparative transcriptome analysis was performed on isolated mesophyll tissue and leaf vein vascular tissue of cassava variety KU50 after mild drought treatment to determine the molecular mechanism behind drought resistance in cassava vasculature. Our results showed that KU50 leaves had increased leaf temperature, with characters of rapidly decreased net photosynthetic rate, stomatal conductance, and transpiration rate in leaves, and the intercellular CO2 concentration accumulated under drought stress. Comparative transcriptome profiling revealed that under drought stress, leaf mesophyll tissue mainly stimulated the biosynthesis of amino acids, glutamic acid metabolism, and starch and sucrose metabolism. In particular, the arginine biosynthesis pathway was significantly enhanced to adapt to the water deficiency in leaf mesophyll tissue. However, in vascular tissue, the response to drought mainly involved ion transmembrane transport, hormone signal transduction, and depolymerization of proteasome. Concretely, ABA signaling and proteasome metabolism, which are involved in ubiquitin regulation, were changed under drought stress in KU50 leaf vascular tissue. Our work highlights that the leaf vasculature and mesophyll in cassava have completely different drought response mechanisms.


Introduction
Drought is a common abiotic stress that can curtail crop productivity (Tuberosa and Salvi 2006). It has many negative effects on plant growth and production. To overcome these negative effects, plants have evolved several bioprocesses to respond to drought stress. These bioprocesses cover many layers of plants' activities. At the biological level, plants adjust the stomata status to reduce water loss under abiotic stress (Steiner et al. 2014). Plants can also optimize carbohydrate concentrations to cope with drought stress (Gong et al. 2013). At the biochemical level, plants can induce antioxidants to scavenge reactive oxygen species under drought stress (Valliyodan and Nguyen 2006). At the molecular level, plants can manipulate many genes to respond to drought stress (Gong et al. 2015).
In the past decades, several studies on drought stress have been conducted. These studies involved different plant species and tissues. The plant species included crop plants, grass plants, and woody plants (Anjum et al. 2017;Cavin and Jump 2017;Li et al. 2017). Moreover, the plant organs that were used for studying drought stress included shoots, roots, and whole plant tissue (Bowne et al. 2012;Chen et al. 2017;Guo et al. 2018). Guo et al. used drought-tolerant and drought-sensitive wheat seedlings to investigate the various physiological processes related to drought tolerance. They found that the wheat metabolome, which could enhance drought tolerance, was dominated by sugars, organic acids, and amino acids (Guo et al. 2018). In addition, seedlings, shoots, and roots had been used as materials to study water deficit stress.
Several studies have been conducted to elucidate how the shoots cope with water deficit situation. For example, Li et al. found that various processes in leaf had been changed when drought stress was coming, such as ribosomes, biosynthesis of secondary metabolites, and carbon metabolism (Li et al. 2016). In another study, the overexpression of the ribosomal protein family member RPL23A in rice could increase water use efficiency under drought conditions (Moin et al. 2017). This study indicated that the shoots could adjust ribosome proteins to enhance drought tolerance. Moreover, Arabidopsis shoots could alter the flavonoid content to protect themselves from abiotic stress (Nakabayashi et al. 2014). This finding indicated that the shoots could change the secondary metabolite to cope with abiotic stress. In addition, carbon metabolism plays an important role in plant leaf abiotic resistance. Durand et al. found that carbon (C) export from leaves to the roots had been enhanced to deal with water deficit situation (Durand et al. 2016).
The shoots will change the metabolism to optimize the consumption of water and nutrients under drought stress, whereas the roots will enhance the uptake of water and nutrients (Gargallo-Garriga et al. 2014). Therefore, the shoots and roots have different responses to drought stress. For example, Wu et al. found that the primary roots of Alhagi sparsifolia had enriched the genes of multiple metabolic pathways under drought stress, especially the glutathione metabolism pathway (Wu et al. 2015). Moreover, the plant root tips have primary signal receptors for water deficiency as they play a role to uptake and transport water from soil to the shoot.
Given the different drought responses of the shoots and roots, the drought resistance mechanism of plant vasculature has become a research hotspot. However, due to the lack of efficient methods to study the vasculature, related studies are few. Compositional analysis of phloem and xylem sap from plants subjected to abiotic stress is a popular method to study how the vasculature is involved in the stress condition (Hu et al. 2016). Because collecting vasculature sap is easy by using EDTA facilitated, insect stylet, and cucurbit exudation. However, concerns have been raised in recent years. For example, contamination is an unavoidable issue related to EDTA-facilitated method and cucurbit exudation (Zhang et al. 2012). Although insect stylet could provide relative pure phloem sap, it is not compatible with many studies due to its trivial amount of sap exudates from the stylet (Reidel et al. 2009). In addition, these methods are associated with a critical issue in that only mobile molecules can be identified. Immobile molecules synthesized in the phloem can play important roles in response to stress (Sevanto 2014). How the vasculature protects against water deficit situation and its molecular mechanism remain unknown.
Cassava (Manihot esculenta Crantz) is an important tropical crop that has high fresh root productivity and extraordinary drought tolerance. Because of its metabolic efficiency under adverse environmental conditions, cassava can produce more carbohydrates per unit area than other crops under conditions of water stress and in poor soils (El-Sharkawy 1993). Therefore, understanding how cassava responds to water deficiency in vasculature and mesophyll tissue of leaf is valuable for improving drought tolerance breeding of cassava and other crops. Physiologically, cassava leaves have decreased water content, ratio of free water content to bound water content (FW/BW), net photosynthetic rate (Pn), intercellular CO 2 concentration (Ci), stomatal conductance (Gs), and transpiration rate (Tr) under drought stress (Shan et al. 2018). Consequently, most photosynthesis-related proteins in cassava leaves are downregulated during drought stress (Chang et al. 2019). However, little is known about vascular tissue response to drought stress in cassava. Here, we performed drought treatment on cassava variety KU50 and sampled the vasculature and mesophyll cells of leaf. By comparatively analyzing the transcriptome data, we found significantly different metabolic responses between the mid-vein vascular bundle and mesophyll cells under slight drought stress. This work will expand our knowledge about how plants adapt to environmental stresses.

Phenotype and Physiological Changes of Cassava Plants Under Drought Stress
Forty-five days after cassava variety KU50 (Fig. 1A) was planted in a greenhouse, water was withdrawn for 3 days, and leaf temperature increased (Fig. 1C). The leaf mid-vein stained with toluidine blue mainly consisted of vasculature tissue, including xylem (which was stained in deep green), phloem (which was around the xylem), and ground tissue (which was in the middle) (Fig. 1B). This structure could guide us to dissect the vascular bundle from leaf by artificial separation. To verify the drought stress intensity, we used the "Flir-One" system to check the leaf temperature (Liu et al. 2011). Thermal imaging showed that the leaf color was red under drought treatment and blue under control.
The red color means that the leaf temperature is high, and blue means that the leaf temperature is low. This result suggested that water stress took place in the leaf of KU50.

Transcripts Annotated Response to Drought Stress
Three repeat leaf mesophyll samples and mid-vein samples under drought stress (DL and DV) and their well-watered  Table 1. We generated a total of 88.47G bases and 589.74 million clean reads.
Because we used a scalpel to dissect leaf mid-vein vasculature to obtain small mesophyll tissues, we used different comparison sets to enrich the specific differentially expressed genes (DEGs, log2 fold change > 2, p < 0.01) in the vasculature. A total of 116 DEGs were commonly involved in the different leaf mesophyll comparison sets: DL compared with CL, CL compared with CV, and DL compared with DV ( Fig. 3A). Simultaneously, 108 DEGs were commonly involved in the different mid-vein vasculature comparison sets: DV compared with CV, CV compared with CL, and DV compared with DL.
We further analyzed the enriched DEGs in mesophyll tissue and mid-vein, found some DEGs significantly enriched in leaf vasculature than leaf mesophyll. They were XETs (Xyloglucan endo-transglycosylase) and SEOs (Sieve element occlusion). XETs are involved in metabolism of xyloglucan and have a function during the formation of secondary cell walls of vascular tissues (Bourquin et al. 2002;Nishikubo et al. 2007). SEOs encode structural phloem proteins (Ernst et al. 2012;Knoblauch et al. 2014).
And we also found the significantly enriched DEGs in leaf mesophyll than leaf vasculature, they were Late embryogenesis abundant (LEA) proteins, which are hydrophilic, mostly intrinsically disordered proteins, and play major roles in desiccation tolerance and expressed at the transcript level in vegetative tissues (Candat et al. 2014). These tissue specific genes' expression levels in our data might be used for quality control to estimate the rate of contamination during sampling (Table 2).

Gene Ontology (GO) Enrichment of DEGs in Mid-Vein and Mesophyll Cells Under Drought Stress
GO annotation of all DEGs specific in the mid-vein and mesophyll was used to generally describe the classification characteristics of responsive genes under drought stress (Fig. 4). Different DEGs were enriched in leaf mesophyll and mid-vein vasculature in response to drought stress. In leaf mesophyll, DEGs were mainly involved in hydrolase activity, nucleoside-triphosphatase activity, cytoskeletal protein binding, pyrophosphatase activity, and threoninetype endopeptidase activity. In the mid-vein, DEGs were mainly involved in the movement of cell or subcellular component, DNA replication, inorganic cation transmembrane transport, proton transmembrane transport, and cellular carbohydrate metabolic process. To cope with drought stress, hydrolase activity was enhanced in leaf mesophyll, and ion transmembrane transport activity was increased in mid-vein.

Main Biological Processes Involved in Leaf Mesophyll and Mid-Vein Vasculature Under Drought Indicated by Kyoto Encyclopedia of Genes and Genomes (KEGG)
KEGG enrichment further enriched DEGs involved in specific biological processes in leaf mesophyll and mid-vein vasculature (Fig. 5). In leaf mesophyll, DEGs related to the biosynthesis of amino acids, such as alanine, aspartic acid, and glutamic acid; pyruvate metabolism; and starch and sucrose metabolism had been significantly enriched under drought stress (Fig. 5A). In mid-vein vasculature, DEGs related to carbon metabolism, plant hormone signal transduction, and proteasome and oxidative phosphorylation had been significantly enriched (Fig. 5B). These results showed that the effects of drought signaling in cassava leaf mesophyll and mid-vein vasculature were significantly different.

KU50 Leaf Mesophyll Significantly Enhanced the Arginine Biosynthesis Pathway to Adapt to Drought
We further investigated the expression changes of genes in the arginine biosynthesis pathway in KU50 leaf mesophyll under drought (Fig. 6). The results showed that most genes were  (Nasibi et al. 2011). These results indicated that enhanced arginine biosynthesis could be the specific adaptive response in leaf mesophyll against drought stress.

Hormone Signal Transduction and Depolymerization of Proteasome Occurred in Vasculature of Mid-Vein in Response to Drought Stress
Two pathways were induced in vascular cells under drought stress, namely, plant hormone biosynthesis and signal transduction, and proteasome-related pathway. The ABA biosynthesis and signaling pathway were highly induced in the mid-vein of KU50 leaves under drought stress (Fig. 7). The ABA biosynthesis genes (NCED and ABA1) and ABA signal elements (PP2C, SnRK2, and ABF) were upregulated. ABA2, NCED, and AAO3 were predominantly observed in vascular bundles (Endo 2008), while PP2C and SnRK2 were found in the plant vasculature (Fujii et al. 2007;Sugimoto et al. 2014). Therefore, our results implied that drought initiated ABA biosynthesis and signaling in vasculature predominantly than in leaf mesophyll. According to proteasome analysis in the mid-vein of KU50 under drought stress, the expression of several genes for the 20S and 19S proteasome was downregulated to adapt to drought stress. Although some genes decreased in the leaf mesophyll, several genes related to the proteasome were induced in the mid-vein. Plant cells contain a mixture of 26S and 20S proteasomes that mediate ubiquitin-dependent and ubiquitin-independent proteolysis, respectively (Kurepa et al. 2009). Ubiquitin is well established as a major modifier of signaling in eukaryotes (Sadanandom et al. 2012). These genes covalently modify target biopathways to regulate drought resistance.

Discussion
Drought is one critical abiotic stress that negatively affects crop development and production. Cassava is a tropical crop that has notable tolerance to abiotic stress and barren soil (Angelov et al. 1993). Cassava can be produced adequately in drought conditions, making it the ideal food security crop in marginal environments (Okogbenin et al. 2013). Thus, several studies on cassava drought resistance have been conducted. However, few have investigated the cassava vasculature response to drought stress. To determine the mechanism behind cassava vasculature's resistance to drought stress, we comparatively analyzed the RNA-Seq data of cassava leaf mesophyll and mid-vein.
To verify the physiological structure of cassava leaf midvein, we stained the leaf mid-vein transection using toluidine blue. The result showed that cassava mid-vein consisted of phloem, xylem, and ground tissue. This result indicated that the cassava mid-vein that we harvested is a completely vasculature tissue.
To determine the compatible drought stress to cassava, we used a thermograph to measure the leaf temperature because  (Liu et al. 2011). The thermograph revealed that drought stress increased the cassava leaf temperature. In addition, the water deficit stress was evaluated by photosynthesis-related parameters, including net photosynthesis, transpiration rate, and stomatal conductance. These parameters all decreased in stressed cassava leaves compared with control leaves. Taken together, these results demonstrated that cassava KU50 is under drought stress (Zargar et al. 2017). The mid-vein from the leaf of stressed KU50 and control sets was dissected, and the tissues were used for comparative transcriptome analysis.
Profiling experiments revealed 8,052 stress-responsive transcripts in cassava leaf and 8,946 stress-responsive transcripts in cassava mid-vein. This result indicated that cassava leaf vascular bundle is more actively responsive to drought than cassava mesophyll cells. According to gene expression pattern analysis for leaf mesophyll and mid-vein vascular bundle, the results implied that leaf mesophyll and mid-vein vascular bundle have completely different drought response patterns.
To further investigate the drought stress response mechanism in leaf mesophyll and mid-vein vascular bundle, GO and KEGG enrichments were performed. Drought stressed KU50 leaf mesophyll had enriched GO pathways related to the following categories: hydrolase activity, cytoskeletal protein binding, and pyrophosphatase activity. These pathways are all involved in plant responsive bioprocesses under abiotic stress. Studies found that numerous hydrolases could be induced by drought, and some of them could improve plant drought resistance. A patatin-like gene encoding a galactolipid acyl hydrolase was stimulated by drought in the leaves of tropical legume and Arabidopsis (Matos et al. 2008). Meanwhile, water deficit triggered extensive AF cytoskeleton reorganization within different types of leaf-blade cells in barley (Śniegowska-Świerk et al. 2015). Moreover, pyrophosphatase activity could increase drought resistance, and the overexpression of a vacuolar pyrophosphatase gene in cotton could enhance drought tolerance (Zhang et al. 2011).
The DEGs in drought stressed KU50 mid-vein vascular bundle were enriched in GO pathways related to the following categories: movement of cell or subcellular component, DNA replication, and proton transmembrane transport. Water deficit stress could stunt growth which consisted with the DNA replication pathway had been decreased in maize (Danilevskaya et al. 2019). Plant membrane transport systems play a significant role when adjusting to water scarcity (Jarzyniak and Jasinski 2014). To sum up, these results suggested that cassava leaf mesophyll has a different molecular pathway compared with mid-vein to respond to the water deficit situation.
To comprehensibly understand the drought response molecular mechanism in tissues of KU50, KEGG analysis was performed. The leaf mesophyll and mid-vein vascular bundle had common molecular pathways, including ATP activity and pyrophosphatase activity. This might be because pyrophosphatase is involved in many bioprocesses during water deficiency (Zhang et al. 2011;Pizzio et al. 2015). KEGG enrichment analysis of leaf mesophyll (Fig. 4A) showed that alanine, aspartic acid, and glutamic acid metabolism and pyruvate metabolism are involved in cassava leaf mesophyll drought response. Asparagine is a notable amino acid in plants for nitrogen storage and transport. It is synthesized by asparagine synthetase via the ATP-dependent reaction (Canales et al. 2012). Proline is a key protective amino acid. Arginine, lysine, and aromatic and branched chain amino acids are all increased under drought stress in different plants to increase drought resistance (Bowne et al. 2012;You et al. 2019). Notably, the arginine biosynthesis pathway was dramatically increased in KU50 leaf mesophyll, which was not significantly changed in mid-vein vascular bundle. Arginine could decrease the oxidant damage caused by water deficit in tomato (Nasibi et al. 2011). On the basis of these results, we hypothesized that cassava leaf mesophyll improved the content of protective components to protect themselves against drought stress.
By contrast, DEGs were mainly sorted in the proteasomerelated pathway, oxidative phosphorylation, and hormone signal transduction pathway in KU50 mid-vein vascular bundle. ABA is a notable phytohormone involved in drought response (Davies et al. 2005). Several ABA responsive genes are involved in cassava vascular drought resistance mechanism. For example, we found that the ABA signal transduction gene PP2C was changed. Studies showed that ABA biosynthesis is initiated in vascular parenchyma and activates a signaling network in neighboring bundle sheath cells (Galvez-Valdivieso et al. 2009). Taken together, we supposed that cassava vascular bundle might optimize the ABA content to induce the expression of downstream genes to resist drought stress. Among these genes, a cluster gene (proteasome-related pathway) was mapped to one pathway, which was significantly decreased in the mid-vein vascular bundle than in leaf mesophyll under drought. The ubiquitinproteasome pathway in plants regulates several cellular signaling processes, such as those elicited by hormones, sucrose, and light, as well as development and responses to pathogen invasion (Smalle and Vierstra 2004). Abiotic stress inhibits 26S proteasome activity, either by directly affecting 26S proteasome functions or by decelerating the turnover rate of other 26S proteasome targets (Ali and Baek 2020).
On the basis of these results, we speculated that cassava vascular bundle has complex drought resistance mechanisms and positive adaptation mechanisms. This may explain the huge difference between cassava vascular bundle drought resistance and cassava leaf mesophyll cell drought resistance (Fig. 8). Further studies are needed to elucidate these mechanisms.

Growth Condition and Drought Treatment
Cassava variety KU50 was used in this study. When cassava seedlings emerged, only one seedling was kept in each pot. The plants were continuously grown in a greenhouse with day/night temperature of 35 °C/28 °C. Drought stress treatment was initiated at approximately 45 days after planting by withdrawing water for 3 days. The control set was watered normally.

Physiological Measurements
Before the treatment, the cassava leaf vein was stained with toluidine blue O and observed using a light microscope. Then, water was withdrawn for 3 days. During the drought stress, the "Flir-One" system (FLIR, Nashua, NH) was used to measure leaf temperature on three leaves localized on top of the cassava seedlings. During the treatment, Li-Cor 6800 (LI-COR, Lincoln, NE) photosynthesis apparatus was used to measure the net photosynthetic rate, transpiration, and stomatal conductance on the top three fully expanded leaves following the manufacturer's instructions. All indexes were collected with three biological replicates and three reads for technical replicates. The results were presented as mean ± SE of independent replicates. All measurements were started at 9 am and ended at 11 am.

RNA Extraction and RNA-Seq Library Construction
We collected three top cassava leaves from each of the control and drought stressed KU50 and then used a scalpel to carefully dissect the mid-vein from the leaf. We grouped mid-veins together as one sample and leaves without midvein together as one sample; each of them had three biological replicates. Three plants were sampled for each biological replicate. All samples were collected in the morning and placed in liquid nitrogen. Total RNA was extracted by using the E.Z.N.A. Plant RNA Kit (Omega Bio-Tek Company, Norcross, GA). RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was checked using a NanoDrop spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE). RNA integrity was assessed using the RNA 6000 Nano Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). About 1 μg of RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer's recommendations, and index codes were added to attribute sequences to each sample.

Transcriptomic Data Analysis
The library preparations were sequenced on an Illumina NovaSeq platform, and 150 bp paired-end reads were generated. Raw data (raw reads) of fastq format were processed through in-house perl scripts. Then, clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N (paired reads with N proportion > 10%), and low-quality reads (the number of bases whose Qphred ≤ 20 accounted for more than 50% of the reads) from raw data. At the same time, the Q20, Q30, and GC contents of clean data were calculated. Manihot esculenta v7.1 was used as the reference genome, and gene model annotation files were downloaded from the genome website directly (https:// data. jgi. doe. gov/ refinedownl oad/ phyto zome? organ ism= Mescu lenta & expan ded= 520). Hisat2 v2.0.5 was used to build the index of the reference genome and to align paired-end clean reads with the reference genome. Then, the FPKM of each gene was calculated on the basis of the length of the gene and read counts mapped to this gene. Differential expression analysis of drought treatment and control KU50 was performed using the DESeq2 R package. Genes with an adjusted p-value less than 0.05 found by DESeq2 were considered differentially expressed.
GO enrichment analysis of DEGs was implemented by the clusterProfiler R package, in which gene length bias was corrected. GO terms with corrected p-value less than 0.05 were considered significantly enriched by differential expressed genes. KEGG is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, organism, and ecosystem, from molecular level information, especially large-scale molecular datasets generated by genome sequencing and other highthroughput experimental technologies (http:// www. genome. jp/ kegg/). We used clusterProfiler R package to test the statistical enrichment of DEGs in KEGG pathways. Availability of Data and Materials All data generated or analyzed during the study are included in this published article and its supplementary information files.

Declarations
Ethics Approval and Consent to Participate Not applicable.

Conflicts of Interest
The authors declare that there is no conflict of interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.