text:
Sequencing and assembly of the Egyptian buffalo genome

Water buffalo (Bubalus bubalis) is an important source of meat and milk in countries with relatively warm weather. Compared to the cattle genome, a little has been done to reveal its genome structure and genomic traits. This is due to the complications stemming from the large genome size, the complexity of the genome, and the high repetitive content. In this paper, we introduce a high-quality draft assembly of the Egyptian water buffalo genome. The Egyptian breed is used as a dual purpose animal (milk/meat). It is distinguished by its adaptability to the local environment, quality of feed changes, as well as its high resistance to diseases. The genome assembly of the Egyptian water buffalo has been achieved using a reference-based assembly workflow. Our workflow significantly reduced the computational complexity of the assembly process, and improved the assembly quality by integrating different public resources. We also compared our assembly to the currently available draft assemblies of water buffalo breeds. A total of 21,128 genes were identified in the produced assembly. A list of milk virgin-related genes; milk pregnancy-related genes; milk lactation-related genes; milk involution-related genes; and milk mastitis-related genes were identified in the assembly. Our results will significantly contribute to a better understanding of the genetics of the Egyptian water buffalo which will eventually support the ongoing breeding efforts and facilitate the future discovery of genes responsible for complex processes of dairy, meat production and disease resistance among other significant traits.

Introduction
Water buffalo (Bubalus bubalis) is an important source of red meat and milk in countries with relatively warm weather, like India, Egypt, Brazil, Italy, and Turkey. The water buffalo population is increasing in these countries, because of its superiority to cattle in terms of higher meat and milk production, more adaptability to harsh environmental conditions, and resistance to a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 diseases. The water buffalo (Bubalus bubalis) belongs to the Bovidae family, Bovinae subfamily, and Bovini tribe. Among others, the tribe Bovini contains the geneses bos, Syncerus, and Bison. The species Syncerus includes the African buffalo (Syncerus caffer), and the species Bison includes the American bison (Bison bison). The species Bubalus bubalis is subdivided into two major subspecies: The river buffalo that has 50 chromosomes and swamp buffalo that has 48 chromosomes [1][2][3].
The water (river) buffalo has been first domesticated in the Indian subcontinent. Then it spread to northern Africa (mainly in Egypt) and to southern Europe. With human migration water buffalo has been domesticated in Italy, Balkan and Turkey, giving them unique characteristics that differs from Egyptian water buffalos [4][5][6][7][8]. This suggests that the Egyptian breed is an intermediate breed between the Indian (Eastern) and European (Italian) ones [2]. The Egyptian water buffaloes are dual purpose animals used for (milk and meat). Local Egyptian breeds were shown to be more efficient when crossed with either the Italian or Pakistani buffaloes under Egyptian local conditions, proving that Egyptian breeds are more adaptive to our managerial practices, climatic and feed changes. Also, they are more resistance to diseases [9,10].
Compared to the cattle genome that is sequenced and well-characterized [11][12][13], little has been done to reveal the genome structure of the water buffalo. This is in spite of its importance as a major member of the Bovidae family and in spite of its interesting physiological traits related to immunology, lactation, and ruminant physiology [14]. The lack of genomic sequences limited further advanced studies addressing functional and evolutionary aspects as well as breeding optimization.
Like any other mammalian genome, the sequencing of the buffalo is complicated by the large genome size, the complexity of the genome, and its high repetitive content. To overcome these obstacles, Next Generation Sequencing (NGS) technology was used to yield a large number of DNA fragments covering the genome with a reasonable depth suitable for the assembly. A number of projects have been recently launched worldwide to sequence the genome of different buffalo breeds. In addition to the one described here, there are large-scale sequencing efforts to sequence the breeds of Italian, Indian, Bangladeshy as well as the African buffalo [15][16][17][18][19].
This manuscript describes the first draft assembly of Egyptian water buffalo genome. The sequencing endeavor was based on multiple runs and libraries of different insert sizes using Sequencing by Oligo Ligation and Detection (SOLiD) technology (Thermo Fisher). The assembly was achieved using a reference-based assembly workflow, highly tuned for SOLiD data. This workflow dramatically reduced the computational complexity of the assembly process and improved the assembly quality by integrating different public resources. The workflow produced an assembly with a scaffold N50 of 3.579 Mbp. Additional analyses, including the genetic content and the identification of milk genes, were also conducted. Our results, together with the other water buffalo assemblies, will significantly help in understanding the genetics of the water buffalo and opens the door to a better understanding of disease susceptibility and to further applications related to the increase of milk and meat production.

Chosen animal
An inbred, female Egyptian water buffalo (Hathour) from the Beheiry region/type was selected for sequencing. To confirm the homozygosity level across the genome and to simplify the genome assembly, the animal was chosen after confirming half-brother and half-sister mating for three generations according to the records kept at Animal Production Research Institute (Kafr El-Shiek, Egypt). Blood samples were collected by a competent veterinary surgeon at her home farm. The experiments in this study were performed in compliance with the official decree of the Ministry of Agriculture in Egypt relevant to animal welfare No. 27/1967 regarding the humane treatment of animals [20]. Ethics approval was obtained from the Ethics Committees of The Agricultural Genetic Engineering Research Institute and Animal Production Research Institute.

Library preparation
DNA was extracted from whole blood using the QIAamp DNA Blood Mini Kit (Qiagen, Hilden, Germany). The buffalo genomic DNA was used to prepare short fragment libraries and mate-paired libraries. Fragment libraries were prepared by shearing two micrograms of DNA to generate fragments of approximately 100-250 bp using the Covaris S2 system (Covaris Inc., Woburn, MA, USA). The amount of the sheared DNA was confirmed by Qubit™ dsDNA HS Assay Kit (Thermo Fisher Scientific, Waltham, MA). Fragment library was prepared with the SOLiD TM fragment library construction kit (Thermo Fisher Scientific), and the SOLiD TM fragment library standard adaptors kit (Thermo Fisher Scientific) following the manufacturer's instructions. Mate-paired libraries were prepared by shearing 5-10 μg of DNA to 1 Kb and to 3 Kb fragments by the Covaris S2 system. Two mate-paired libraries, 1 Kb and 3 Kb, were generated using the 5500 SOLiD TM Mate-Paired library kit (Thermo Fisher Scientific).

Library amplification and sequencing
Each library was amplified using emulsion PCR. 35 femto moles or 20.4 billion molecules (~5-6 ng) of a fragment library were emulsified with mineral oil in a SOLiD EZ Bead Emulsifier according to the E120 scale protocol (Thermo Fisher Scientific). The emulsion PCR was carried out in a SOLiD EZ Bead Amplifier (Thermo Fisher Scientific) using the E120 setting. To enrich for the beads that carried amplified template DNA, the beads were purified on a SOLiD EZ Bead Enricher using the recommended chemicals and software (Thermo Fisher Scientific). The beads were then purified, quantified, deposited, and attached to a sequencing flowchip. The flowchip was mounted on the 5500xl sequencer and the library was sequenced. The short fragment library was sequenced using the settings and recommended chemicals for sequencing 75 nucleotides in the forward direction and 35 nucleotides in the reverse direction, with Exact Call Chemistry (ECC) module (Thermo Fisher Scientific). Mate-Paired libraries were sequenced using, 2 x 60-bp reads were on a SOLiD 5500xl instrument, with Exact Call Chemistry (ECC) module (Thermo Fisher Scientific). Fig 1 shows the basic steps of the analysis workflow (pipeline). Both Single-end reads and mate-paired reads were mapped to the cattle reference genome (Bos_taurus_UMD_3.1.1). The bowtie software [21] was used to map the highly similar reads to the genome. The SHRiMP software [22], which takes color space into account, was used to map the less similar reads to the reference genome. The contigs were constructed by stitching the reads mapped to the same cattle locus and overlapping with at least 50 bps. The reads that could not be mapped with Bowtie or SHRiMP were collected and assembled separately using the de novo program velvet [23].

Genome assembly
The mate-paired reads were then used to confirm the assembled contigs, and to assemble the contigs into larger scaffolds. A contig is confirmed if the two parts of the paired-end reads are mapped to it. Paired-end reads mapped to two different contigs are used in scaffolding the two contigs. The greedy path-merging algorithm, like the one implemented in Opera, Sopra, and Grass [24][25][26][27] was used in the scaffolding process.
The S1 File includes more details about the implementation of the assembly workflow, which is based on the Tavaxy Workflow Management System [28] and the distributable version available as Docker image. Web links to the workflow implementation is also available at the project website https://egyptianbuffalo.org/ and https://buffalo.tavaxy.org.

Genome annotation
Repeat analysis. RepeatMasker [29], with default parameters, was used to locate the repeats in the assembly. The resulting output was analyzed using in-house scripts to compute coverage of repeats, classify them into families, and produce final statistics and charts. Gene analysis. A pipeline similar to the NCBI Eukaryotic Genome Annotation Pipeline [30], was used to annotate the assembly. The coding genes in the Egyptian buffalo genome were identified over three phases. During the first phase, the mammalian mRNA sequences available in the NCBI were aligned to the draft genome. The mammalian mRNA sequences included RefSeq mRNA data from cattle, human, mouse, rat, and pig. The alignment of mRNA to the Egyptian water buffalo genome was achieved using BLAT [31], with a cross-species setting of 80% identity and a minimum hit length of 50b. The second phase was performed using homology-based prediction by transferring the cattle's annotation to our assembly. The program GASS [32] was then used to map the gene sequences to the draft genome. In the final stage, extra genes were identified using the de novo prediction Augustus program [33].
Comparison to other assemblies. The assembly of each draft genome was collected from the NCBI Assembly database. Each assembly was then processed using RepeatMasker [29] to mask out the repeated regions. The alignment between the different assemblies was achieved using Minimap2 [34], CoCoNUT [35] and Mauve [36], where the minimum length of a reported alignment is 100bp. Colinear alignments were coalesced into larger syntenic blocks. Water buffalo genome assembly workflow. Both mate-paired and short fragment reads, after trimming at Q15, were mapped to the cattle genome. The mapped reads were used to construct contigs. The un-mapped reads were collected and assembled separately using the de novo program velvet. The mate-paired reads were used to confirm assembled contigs, and to assemble the contigs into larger scaffolds. The greedy path-merging algorithm, Opera, Sopra, and Grass algorithms were used in the scaffolding process. https://doi.org/10.1371/journal.pone.0237087.g001 To facilitate the plotting of the synteny blocks, we used the chromosome-level assembly of the cattle and UOA buffalo genomes as references. The plots were generated using the Rideogram package in R.

Results and discussion
In this paper, the Egyptian water buffalo genome was sequenced using the SOLiD technology. To the best of our knowledge, this is the first mammalian genome to be sequenced using this technology. Two-base encoding offered by the SOLiD system is critical for robust detection of sequencing errors to produce highly accurate base calls necessary for genome sequencing [37]. Reference-based assembly strategy was used to produce the assembly, this strategy can overcome the short read length produced by this sequencing platform, speed up the computation, and produce a high-quality assembly. Specifically, the mapping step including the color-space based SHRiMP program [22] could map most of the data to the reference genome and facilitated the assembly of the contigs. The use of de novo assembler was then limited to a muchreduced set of unmapped reads. The reported assembly has been verified using different information sources, including mate-paired reads and gene sequences from public databases.
The incorporation of sequences from both mate-paired and single-read libraries significantly improved the production of long scaffolds. In addition to their importance in genome annotation, the use of publicly available gene sequences from related organisms also helped to verify the genome assembly and improve its quality by correcting some scaffolds.

Sequencing
For sequencing, a female Egyptian water buffalo was selected from the Beheiry region/type. DNA was extracted from peripheral venous blood, and two short fragment libraries, one 1 Kb mate-paired, and one 3 Kb mate-paired library were constructed. The genome was sequenced using the SOLiD DNA sequencing platform (SOLiD 5500XL sequencer). In total, 3.4 billion reads with 229.24 Gb were generated, achieving about 80X-fold coverage of the water buffalo genome. After quality trimming at Q15, 2028 million reads (~60%) of total length 100.85 Gb passed the filtering, achieving about 36-fold coverage of the water buffalo genome. Table 1 summarizes the resulting sequencing data and related quality metrics. The reads of different runs were submitted to NCBI under project number PRJNA267486, and a copy of it is accessible from our website http://egyptianbuffalo.org.

Short-read de novo sequencing and assembly
The river buffalo genome has 5 bi-armed and 19 acrocentric chromosome pairs in addition to the XY sex chromosomes [14,38]. The water buffalo genome was assembled using a referencebased assembly workflow using the "Tavaxy" workflow management system [28]. Tavaxy is an extended version of Galaxy enabling the execution of the workflow using cloud computing. The use of reference-based assembly workflow was preferred over the de novo assembly workflow due to the short read length produced by the SOLiD platform; which will not be easily assembled into longer contigs using de novo assembler. Furthermore, the abundance of repeats in the genome leads to an ambiguity that cannot be solved in the assembled regions using short reads/contigs. The cattle genome (Bos_taurus_UMD_3.1.1) was selected as the reference genome. At the chromosomal level, the reads of the Egyptian buffalo covered all chromosomes with an average of 98.1%. Also, the assembled contigs from the mapped reads covered 97.5% of the cattle genome (S2 File) these primary results justified the use of the cattle genome as a reference for the assembly.

Assembly results
Raw reads were successfully assembled into 447,289 contigs with a total length of 3.005 Gb ( Table 2). The maximum contig length was 148 Kb, and the N50 of contigs was 14.568 Kb. The number of contig lengths larger than 1 Kb was 298,457 with a total length of about 2.475 Gb. Mate-paired read libraries of different lengths (1 Kb and 3 Kb) have been used to confirm the contigs. A contig is considered verified if both ends of at least three mate-paired reads are mapped to it and there are no mate-paired reads connecting this contig to another one in a non-suffix prefix manner.
The scaffolds have been confirmed using mate-paired read libraries of different lengths (1 Kb and 3 Kb), and also by aligning the publically available cattle mRNA genes to the draft assembly. Out of the tested 11653 cattle genes, 11138 genes were mapped correctly with complete conservation of the gene structure within one scaffold.

Comparison to other buffalo assemblies
A comparison of the Egyptian buffalo to all available buffalo assemblies deposited in the NCBI database until June 2020 was performed. These assemblies included the Mediterranean UOA_WB assembly [17], the Bangladeshy assembly [18], the Mediterranean UMD_Caspur [15,39], the Indian [19], and the African buffalo assembly [16] (Table 3). All new assemblies appeared in 2019 and 2020, and the UMD water buffalo and the African buffalo assemblies were updated within the last year using the recent long read sequencing methods.
Comparisons of the Egyptian buffalo to these water and African buffalo assemblies as well as to the cattle genome were performed (Table 3B). At the contig level, the Egyptian assembly has a moderate N50 compared to other Bubalus bubalis draft assemblies, which is expected due to the short-read output of the SOLiD technology. Comparing the scaffolds, the assembly of the Egyptian buffalo had N50 of 3.57 MB which is the second largest N50 among the Bubalus bubalis draft assemblies based on short-read sequencing. The usage of mate-paired read libraries, aligning the publically available cattle mRNA genes to draft assembly, and the usage of UMD assembly are the main reasons for refining our assembly.
Irrespective of the differences in the number of scaffolds, the Egyptian draft assembly is in good agreement with the Mediterranean buffalo (UOA and UMD assemblies). At the read level, the reads of the draft Egyptian assembly covered 99.59% of the UMD/UOA assembly. At the scaffolding/contig level, 99.3% of the Egyptian buffalo assembly is mapped to the UMD/ UOA assembly.
The synteny blocks among the draft assemblies were computed, taking the cattle and UOA assemblies as references. The cattle and UOA genomes are currently the reference genome assemblies: They were assembled up to the level of chromosome and they were intensively verified using different support methods like Hi-C, RH and optical mapping. Fig 2 shows different plots of the synteny blocks with the UOA and Cattle genomes as references. The plots are based on plotting syntenic blocks of at least 1Mbp, except for the Indian buffalo where we had to reduce the threshold to 500Kbp. As expected, the new UOA Mediterranean assembly [17] has better placement of scaffolds and chromosomal structure. It is in high synteny with the cattle genome. It can also be observed that the assemblies based on short read sequencing, notably the Bangaladeshi and ours, could achieve good level of synteny as did the new UOA assembly.  Future work to enhance these assemblies using other runs of long read assemblies would definitely lead to significant improvement in scaffolding.

Repeat analysis
Repeats in the Egyptian water buffalo genome assembly were identified using the program RepeatMasker [29]. Analysis of the repeats revealed that a 1.2 Gb, 39.3% based upon genome size 3.05 Gb, of the Egyptian water buffalo genome is composed of repeats. The repeat content in Egyptian water buffalo assembly is less than that of the cattle genome (1.27 Gb). The LINE family comprises around 49% of the repeats (46% in the cattle), followed by the SINE family representing 36.8% (35.6% in the cattle). LTR makes up 10% in both water buffalo and cattle genomes (Table 4).

Gene content
Annotation of genes, transcripts, was performed using the NCBI Eukaryotic Genome Annotation Pipeline [30]. The annotation of the Egyptian water buffalo genome yielded 21128 genes (S3 File). Of these, 11,702 ones came directly from the mammalian mRNAs. A total of 9426 genes that have no mRNA hit were identified using Augustus annotations. A total of 3753 pseudogene loci from 603 genes were identified (S4 File). In addition, there were 5393 non-coding RNAs (1202 rRNA, 2242 tRNA, and 1949 other small RNAs), identified by RepeatMasker. Milk genes. Buffalo milk accounts for 13% of the total world milk production (http:// faostat.fao.org/). Although the buffalo provides much lower average milk yield than that of Holstein cows [40], its milk is one of the best raw materials for making dairy products [41]. Buffalo milk has higher fat content, higher unsaturated fatty acids, higher milk protein and lower levels of cholesterol compared to that's of the cow's milk [42].
Although milk composition is variable across species, the milk and mammary genes are more conserved in most mammals than are other genes in the genome. Therefore, milk genes were identified in our assembly using the annotation of milk lactation gene sets [13]. Collectively, 3889 e milk virgin-related genes; 1,383 milk pregnancy-related genes; 3,111 milk lactation-related genes; 867 milk involution-related genes; and 840 milk mastitis-related genes were identified inside the Egyptian draft assembly (Fig 3 and S5 File).

Pathway and gene ontology analysis of milk genes
The molecular functions of milk virgin-related genes; milk pregnancy-related genes; milk lactation-related genes; milk involution-related genes; and milk mastitis-related genes were mapped using PANTHER pathway analysis kit [43]. Each milk lactation gene-set was classified according to their possible involvement in molecular functions, biological process, cellular component, and different pathways (Fig 3 and S6 File). The identified genes encompass a wide spectrum of molecular functions, biological processes, pathways, cellular components and Panther protein classes. The genes were successfully annotated and classified into the following categories: binding, proteins with different catalytic activity, regulators of molecular function, molecular transducer activity related proteins, structural molecule activity, transcription regulator activity, translation regulator activity, and transporter activity.

Conclusions
The availability of the Egyptian buffalo genome is an addition to the global efforts for unravelling the genomic basis of the buffalo biology. The paper at hand introduces a high quality draft assembly based on short-read sequencing using the SOLiD platform. To the best of our knowledge, this is the first large mammalian genome to be sequenced using this technology. Combining libraries of single and paired-end reads were crucial to refine the assembly and produce long scaffolds. The raw sequences and the assembled contigs and scaffolds have been made public and deposited at the NCBI repository. Genome annotation was conducted and gene lists have been produced. We focused our analysis on milk genes and analysed the respective pathways and gene ontology. The usage of modern long-read sequencing has shown great advantage in producing better scaffolds, as demonstrated in the recent Italian and African buffalo assemblies and the African buffalo. We plan to use similar technology in the future to improve the scaffolding and provide well-phased reference assembly for the Egyptian breed.
We expect that the next logical step in buffalo genome research is to run large scale genome variation projects on a large number of animals and to conduct association studies to predict genes and variations contributing to the unique characteristics of the buffalo breeds. This step will make use of the available well-phased reference assemblies to efficiently compute these variations.
The Egyptian water buffalo genome offers a unique insight for understanding the variations in sequence between buffalo around the world and cattle imporatant genes. This variation could help in the future discovery of genes behind complex dairy and meat production traits and to understand the species-specificity of milk composition. Moreover, the availability of various buffalo assemblies with gene annotation will offer a good opportunity for understanding the disease resistance/susceptibility in buffaloes. This work will eventually help in identifying and validating specific genes related to the quantitative and qualitative aspects of milk production. Such tools will enable the farmer to predict the characteristics of the desired animal and be able to invest in appropriate animals with higher milk or meat production without wasting effort, time and money.
Supporting information S1 File. Description of implementing the assembly workflow on the Tavaxy workflow management system.


extfieldsofstudy: Medicine Chemistry

provenance: AgriculturalAndFoodSciences,2019-2020-Agricultural And Food Sciences-2020.jsonl.gz.gz

text:
Bioactive Compounds from Leaves and Twigs of Guayule Grown in a Mediterranean Environment.

Guayule (Parthenium argentatum) is one of the most promising alternatives to produce natural rubber. As the guayule leaves represent a waste byproduct in the rubber extraction process, their exploitation might contribute to the valorization of the whole plant biomass. The specific aim of our study was to determine the antioxidant capacity and the content of phenolic compounds in leaves and twigs of different guayule lines cultivated in a Mediterranean environment. The antioxidant capacity and the contents of phenolic were affected significantly by guayule accession and harvest time. Overall means in twigs represented about 60% of the corresponding leaf values for antioxidant capacity as DPPH (1,1-diphenyl-2-picrylhydrazyl) and not tannic phenolic and about 55% as ABTS (2,2'-azinobis (3-ethylbenzothiazoline-6-sulphonic acid) diammonium salt), total phenolic, tannic phenolic and total flavonoid contents. Five individual phenolic compounds were identified in methanolic extracts of leaves. Neochlorogenic acid (62.5-174.8 mg g-1 DW) and chlorogenic acid (28.6-41.1 mg g-1 DW) were the most abundant phenolic acids. In addition to the compounds quantified in leaves, cynarin was identified only in twigs and for the first time in guayule biomass. Our results indicate that guayule leaves and twigs represent a rich source of antioxidants for potential applications in fodder, nutraceutical and pharmacological sectors.

Introduction
Natural rubber, for its excellent and unique properties, represents a strategic raw material whose demand is increasing worldwide [1]. Traditionally, it is extracted as latex from the plantations of Hevea brasiliensis (Willd. ex A. Juss.) Müll. Arg, which are mainly grown in tropical areas of Asian-Pacific countries and supply 90% of the global production [2].
However, potential shortages in natural rubber, caused by several factors (fungi attacks, land use change in production areas, vulnerability to climate change, etc.) and the onset of allergies due the presence of some proteins in its latex affect availability and use of Hevea rubber [3]. Since 1997, European Union has listed natural rubber as a critical material and imports 1.2 million tons each year [4].
These shortcomings and the rising global demand in natural rubber have cyclically re-addressed interest in alternative sources, of which the most important are represented by guayule (Parthenium argentatum A.Gray) and Kazakh dandelion (Taraxacum kok-sagyz L.E.Rodin) [3].
Guayule is a shrub from the Asteraceae family native to Northern Mexico and Southwest regions of the USA. As it grows within a range of temperature from −15 to +40 • C and an annual rainfall requirement of 350 to 800 mm [5], guayule is well suited to semi-arid and Mediterranean regions.
Currently, guayule is considered as an alternative source, having the potential to become commercial crop, to produce natural rubber in Europe. Sfeir et al. [6] evaluated technical and economic feasibility of a guayule commodity chain in Mediterranean Europe, indicating that this crop might be a sector commercially profitable only in the framework of a biorefinery approach. This requires the valorization of the whole plant biomass and that of the production of bioenergy and non-rubber co-products, to be converted into products and/or raw materials with high added value for different sectors [7]. Therefore, it is justified an overall interest in developing new agricultural commodities such as guayule for Mediterranean or semi-arid climates over the next years [8].
As guayule leaves contain negligible amounts of rubber compared to branches, they represent a waste by-product [9] in the rubber extraction processes from its biomass. In fact, the complete removal of leaves before homogenization is a necessary step to produce high-quality latex, avoiding fine solids contamination mostly attributable to the co-homogenized leaves [10]. As regards the presence of bioactive compounds in leaves of guayule, Mears [11] reported the distribution of 34 flavonoids detected in leaves of North American Parthenium spp. A recent paper focused on procedures for the extraction of bioactive compounds from leaves of guayule [12].
Despite the abundance of literature dealing with studies on guayule for rubber production and the shared approach for the exploitation of all biomass fractions, very little information is available regarding secondary plant metabolites and, in particular, the phenolic fraction contained in its leaves. Concurrently, no information dealing with the phenolic fraction in twigs (i.e., young thin stems) of guayule, as well as the possible influence of genotype and harvest time seems obtainable.
Plant bioactive compounds have been widely studied as natural phytochemicals due to superior antioxidant properties compared with synthetic antioxidants. Indeed, bioactive compounds from leaf extracts of tree species have been investigated for the beneficial effects (antioxidant, anticancer, anti-inflammatory, antidiabetic, etc.) on human health [13]. Therefore, the investigation, identification and quantification of the bioactive compounds of leaf extracts are highly important for its antioxidant properties and related health effects, of which phenolic compounds are primarily responsible for [14]. Additionally, phenolic compounds play a role in plant defense reactions against biotic and abiotic stresses, as well as in animal nutrition and welfare. [15,16].
Within a research framework started in Sardinia (Italy) to evaluate adaptation and performances of guayule lines cultivated under Mediterranean conditions, our study focused on the exploitation of guayule leaves and twigs, as source of phenolic compounds and antioxidants for potential applications in fodder, nutraceutical and pharmacological sectors. Therefore, the specific aims of our study were to determine, in the above-mentioned plant fractions, the antioxidant capacity and content of phenolic compounds as well as their variations due to different genotypes and harvest times.

Antioxidant Capacity and Phenolic Content
Phenolic compounds are secondary plant metabolites, which have been classified into several classes, namely, flavonoids, phenolic acids, stilbenes and lignans, according to their chemical structure. They are of pivotal importance in regulating plant defenses mechanisms against several stresses and in healing wounds of rubber plants. Additionally, many studies have explored the pharmacological benefits of phenolic compounds, including the high antioxidant potential and other health-promoting properties [14]. In leaves, the antioxidant capacity and the contents of total phenolics (TotP), non-tannic phenolics (NTP) and total flavonoids (TotF) were affected significantly by both guayule accession and harvest time (Table 1). Line P803 showed the highest antioxidant capacity values in winter. The content of TotP ranged from 147 to about 204 g gallic acid equivalent (GAE) kg −1 in P803, the only accession where the effect of harvest time was significant. The same trend was observed for NTP values.
Overall average for tannic phenolics (TP) content was unaffected by line and season. AZ-1, P803 and 11591 exhibited the highest TotF values in winter and spring, whereas a marked seasonal increase was noticed for AZ-2.
In twigs, the harvest time did not affect ABTS (2,2 -azinobis (3-ethylbenzothiazoline-6-sulphonic acid) diammonium salt) and NTP values, whereas TP values were irrespective of guayule line (Table 2). At both harvest times, the highest values for the antioxidant capacity and the contents of TotP were found in AZ-1 and AZ-2. The peak values of NTP were reached in AZ-2 and AZ-1, in winter and spring, respectively. A remarkable seasonal increase in TotF content was recorded for P803. Overall means in twigs represented about 60% of the corresponding leaf values for antioxidant capacity (DPPH) and NTP and about 55% for ABTS, TotP, TP and TotF contents (Tables 1 and 2).
For guayule leaves, significant correlations were found between the antioxidant capacity by means of ABTS and DPPH (1,1-diphenyl-2-picrylhydrazyl) methods and the phenolic contents (Table 3). DPPH and total phenolics showed a highly significant correlation (R 2 , 0.6774). Significant correlations were also established between antioxidant capacity and NTP contents. Regarding twigs, closer correlations were found than in leaves. In fact, significant correlations were found between the antioxidant capacity by means of the ABTS and DPPH methods and the phenolic content reaching a significant correlation of R 2 , 0.8581 and 0.9174, respectively (Table 3).
Our results agree with other studies that reported the relationship between antioxidant capacity and total phenolic compounds in species belonging to the Asteraceae family such as Kazakh dandelion and Taraxacum officinale Weber ex Wiggers, as well as in many representative herbs within Labiatae and Asteraceae [17][18][19][20]. Unfortunately, very few papers deal with chemical composition and antioxidant capacity in leaves of guayule. Banigan et al. [21] investigated crude protein content of deresinated and derubberized leaf tissue and the amino acid composition. Mouse feeding studies showed deresinated leaf meal to be nontoxic but deficient as a sole source protein. Leaf bromatological composition of guayule is being investigated by Sulas et al. [22]. Mears [11] determined the distribution of flavonoids detected in the North American species belonging to Parthenium genus. P. argentatum showed a flavonoid pattern rich in aglycones and glycosides: derivatives of quercetagetin, 6-hydroxykaempferol, quercetin, kaempferol and kaempferol 3-methyl ether. Table 1. Trolox equivalent antioxidant capacity (TEAC) by 2,2 -azinobis (3-ethylbenzothiazoline-6-sulphonic acid) diammonium salt (ABTS) and 1,1-diphenyl-2picrylhydrazyl (DPPH) methods, total phenolics (TotP), non-tannic phenolics (NTP), tannic phenolics (TP), total flavonoids (TotF) in leaves of the examined accessions of Parthenium argentatum harvested at 30 (winter) or 33 months (spring) after transplanting.   It is noteworthy that our values for total phenolics reported in Table 1 were up to 7-fold or even 11-fold higher than leaf peak values indicated by Piana et al. [12]. According to that paper, dry guayule leaf material from Mexico was tested for both solid-liquid and ultrasounds-assisted extractions according to a range of conditions for green solvent composition and extraction time. Even if the total phenolic content of extracts was quantified by the Folin-Ciocalteu's assay as in our study, the magnitude of differences might be, presumably, explained by differences in the tested genotype, growth environment, sampling and storage procedures. Conversely, our samples were collected on site and immediately frozen in liquid nitrogen, then stored at −80 • C, and afterwards freeze-dried. An important implication arising from this direct comparison is that the potential of guayule leaves as source of total phenolic and antioxidants, based on our results, is markedly higher than that reported so far. For flavonoids quantification, a different method and standard were used, so a quantitative comparison of data is difficult. Regarding comparison with other species belonging to same genus, Panwar et al. [23] found in the weed Parthenium hysterophorus L. a content of phenolics of 20.8 mg GAE g −1 , but they performed the quantification with an alkaline hydrolysis method. Within the same botanical family, Molinu et al. [20] found in Taraxacum kok-saghys leaves a content of antioxidant capacity, total phenolics, total flavonoids of 19.6 mmol Trolox equivalent antioxidant capacity (TEAC) 100 g −1 DW, 106.4 g GAE kg −1 DW, 22.9 g catechin equivalent (CE) kg −1 DW respectively, lower than our results. A study dealing with bioactive compounds and antioxidant capacity from Mediterranean garland (Chrysanthemum coronarium L.) showed in leaves a content of antioxidant capacity (14.4 mmol TEAC 100 g −1 ) lower than our results, and total flavonoids (32 g CE kg −1 DW) very close to our results [24]. Sardinian milk thistle (Silybum marianum (L.) Gaertn.) populations were studied for chemical and productive properties by Sulas et al. [19] reporting, in acetone/water leaves extract, a content of antioxidant capacity (7.6 mmol TEAC 100 g −1 DW) and total flavonoids (11.5 g CE kg −1 DW) lower than our results. A content of phenolic compounds lower (18.7 mg GAE g −1 DW) than our results was also reported by Benchaachoua et al. [25] in a aqueous-methanol extract of leaves of S. marianum collected in Algeria. Ben Salem et al. [26] found phenolic contents of artichoke (Cynara scolymus L.) leaves extracts from 39.9 to 54.5 mg GAE g −1 DW lower than our results. Outside the Asteraceae family, a comparison with the legume shrub tagasaste (Chamaecytisus proliferus var. palmensis) grown in the same site than guayule, revealed that, at late autumn cutting, the antioxidant capacity and total flavonoids were quite similar to guayule results for leaves, but different for twigs [27]. Overall, comparisons with other study regarding the similar species are very difficult due to variations in methods, procedures and standard used for the analysis quantification According to the literature, seasonal changes in the phenolic content of plant leaves might depends on species, site conditions (altitude, solar radiation), the position of the leaves and variety, etc. [28][29][30].

RP-HPLC Analysis of Phenolic Compounds
Twenty-eight individual phenolic compounds were screened by RP-HPLC based on the literature concerning phenolic composition in Parthenium genus and other members of the Asteraceae family. Five individual phenolic compounds were identified in leaf methanolic extracts of guayule. Identification of compounds was carried out by comparison with their corresponding standard and with data reported by other authors. The five individual compounds identified in leaves were neochlorogenic (5-caffeoylquinic acid), chlorogenic (3-caffeoylquinic acid), cryptochlorogenic (4-caffeoylquinic acid), isovanillic, and 3,5-DCQ (3,5-Di-O-caffeoylquinic acid) acids (Table 4). Statistical analysis evidenced variations due to accession and harvest time. The neochlorogenic (62.5-174.8 mg g −1 DW) and chlorogenic acids (28.6-41.1 mg g −1 DW) were the most abundant phenolic acids, which taken together contributed on average to about 92% of the quantified total cinnamic acid derivatives. The neochlorogenic, cryptochlorogenic and 3,5-DCQ acids were significantly affected by harvest time in AZ-1, P803 and AZ-2, respectively, whereas chlorogenic and isovanillic acids were irrespective of harvest time. In particular, the spring content of neochlorogenic acid only represented 38% of the corresponding winter value in P803, whereas the spring content of 3,5-DCQ showed a 72% increase compared to the winter value in AZ-2. In guayule leaves, Piana et al. [12] identified chlorogenic and isovanillic acids as our results, and others caffeoylquinic acid isomers, identified in our study as criptochlorogenic and neochlorogenic acids. In addition, the same authors detected different dicaffeoylquinic acid isomers, identified as 3,5-DCQ and cynarin (1,3-Di-caffeoylquinic acid) in our study. Chlorogenic, caffeic and vanillic acids were also detected in leaves and flowers of P. hysterophorus grown in different regions of Ethiopia [31]. On the contrary, caffeic acid was not detected in our study. Moreover, a review of Khan and Ahnmad [32] reported the presence of chlorogenic, neochlorogenic, vanillic acids in leaves of P. hysterophorus. Chlorogenic acid (lower than 0.2 µg g −1 ) was also detected in the ethanolic extract of aerial parts of Chrysanthemum parthenium Bernh, very low value, compared to our results, probably due to a different solvent extraction [33]. A study regarding Bellis perennis L. revealed that a greater amount of chlorogenic and neochlorogenic acids were detected in leaves than in root; these compounds were studied for their allelopathic properties [34]. Table 4. HPLC analysis of phenolic acids (mg g −1 DM) in leaves of the examined accessions of Parthenium argentatum harvested at 30 (winter) or 33 months (spring) after transplanting. (Means ± SD, n = 3, Fisher's test).

Neochlorogenic Acid
Chlorogenic Acid Cryptochlorogenic Acid Isovanillic Acid 3,5-DCQ [ Derivatives of hydroxycinnamic acid, and particularly chlorogenic acid, are considered a characteristic trait of the Asteraceae plants, although they are ubiquitous of the plant kingdom [34,35]. Moreover, chlorogenic acid is an important and biologically active dietary polyphenol, playing several important and therapeutic roles in human health, such as antioxidant, antibacterial, hepatoprotective, cardioprotective, anti-inflammatory, antipyretic, neuroprotective, anti-obesity, antiviral, anti-microbial, anti-hypertension, free radical scavenging and central nervous system stimulator activity [36]. The same authors referred that chlorogenic acid has shown hepato-protective effects in animals; its positive effect was also revealed on performance and antioxidant function in piglets [37]. Cynarin was identified only in guayule twigs, in addition to the same five individual compounds quantified in leaves (Table 5), and for the first time in guayule biomass. Again, statistical analysis evidenced variations due to accession and harvest time. Overall, a dramatic decrease in neochlorogenic acid content in twigs compared to leaves was detected. Taken together chlorogenic and 3,5-DCQ acids contributed to about 95% of the total cinnamic acid derivatives that were quantified in guayule twigs. Spring reduction were recorded in all accessions for neochlorogenic acid, in AZ-2 and 11591 for chlorogenic acid, in AZ-1 and 11591 for isovanillic acid and in 11591 for 3,5-DCQ acid, respectively.
Conversely, within the Asteraceae family, cynarin was detected in methanolic extract of artichoke leaves as reported by Farhan et al. [38]. This compound was also identified in methanolic extract of sunflower (Helianthus annuus L.) sprouts, showing antioxidant and antidiabetic effects in vitro models [39]. It is noteworthy that cynarin has even been reported as anti-HIV active [40].
In leaves, the antioxidant capacity with both ABTS and DPPH assays was positively correlated with neochlorogenic and isovanillic acids, whereas chlorogenic acid was correlated with Ptot and NTP ( Table 6). The relatively higher concentration of neochlorogenic acid was the most probable determinant for observed antioxidant activities in leaves. In twigs, all the six identified compounds were weakly correlated with ABTS and four of them with DPPH. Only chlorogenic and 3,5-DCQ acids were also correlated with TotP, NTP, TP and TotF. According to Kozyra et al. [41], chlorogenic acid significantly influence the antioxidant activity of methanolic extracts from inflorescences of Carduus sp. In Loquat (Eriobotrya japonica Thunb. Lindl.) fruits, contents of chlorogenic and neochlorogenic acids were significantly correlated with antioxidant capacity, assessed by DPPH and ABTS radical scavenging activity [42].

Materials and Methods
The field experiments were carried out in Southern Sardinia (39 • 31 N, 8 • 51 E, Italy), where the climate is Mediterranean with mild winter. The area has a long-term average annual rainfall of 446 mm received mainly in autumn and winter months, and a mean annual air temperature of 17.6 • C. The soil, classified as Typic Fluvaquents, is sandy-clay-loam, with pH 7.8, average nitrogen content of 1.1% and phosphorous 16.2 ppm.
Four guayule accessions kindly provided by the Agricultural Research Service of the United States Department of Agriculture (ARS-USDA) were evaluated: the improved guayule lines AZ-1 and AZ-2, the accession P803, and the line 11591, which was developed during the Emergency Rubber Project and used as control. Additional details on the above-mentioned plant materials are reported by Ilut et al. [43]. Seeds were treated to break dormancy using the method of Naqvi and Hanson [44] and were sown in wet filter paper inside a growth chamber. Two weeks after germination, seedlings were moved into jiffy pots and afterwards put in plastic bags. Plantlets were then moved outdoor and transplanted in the field in May 2015. The experimental design was a randomized complete block with three replicates. Seedlings were transplanted in row spaced at 1 m and 0.5 m within the row. Each plot consisted of eight rows (two for each line). Additional outer rows were also transplanted as border. Seedlings were irrigated daily during the first week after transplanting and then at ten-day intervals in summer months from June to early October by using a drip irrigation system. Fertilizer was applied only at the time of field preparation when 30 kg N ha −1 was incorporated into the soil.

Leaves and Twigs Sampling
At 30 and 33 months after transplanting, corresponding to winter 2017 and spring 2018, respectively, three plants of each line per plot were cut at ground level and harvested. Aboveground plant dry matter was partitioned into leaves, twigs (i.e., the youngest apical stem portion, 6 to 8 mm in diameter, developed in last growing season), and remaining branches. Leaf and twigs samples were immediately frozen in liquid nitrogen, and stored at −80 • C, until lyophilization with Heto Lyolab 3000 (Heto-Holten A/S, Allerød, Denmark) for 48 h (−55 • C). Leaf and twig subsamples were ground in a mill to a fine powder and stored at −20 • C until analysis. Sample preparation procedures were performed according to Molinu et al. [20]. All samples were analyzed in triplicate.

Phenolic Content and Antioxidant Capacity
Total phenolics (TotP), non-tannic phenolics (NTP) and tannic phenolics (TP) of extracts were determined using the Folin-Ciocalteau reagent, according to procedures previously described [30]. Results were expressed as g of gallic acid equivalent (GAE) kg −1 dry matter of plant material (g GAE kg −1 DM) by means of a calibration curve of gallic acid. Total flavonoids (TotF) were quantified by colorimetric assay with the AlCl 3 method, following procedures previously reported [30]. TotF were quantified by catechin calibration curve and results were expressed as g of catechin equivalent (CE) kg −1 dry matter (g CE kg −1 DM). Antioxidant capacity was evaluated by ABTS (2,2 -azinobis (3-ethylbenzothiazoline-6-sulphonic acid) diammonium salt) and by DPPH (1,1-diphenyl-2-picrylhydrazyl) assays [30]. Trolox (6-hydroxy-2,5,7,8-tetramethylchroman-2carboxylic acid) was used as the reference standard. For each assay, 0.1 mL of appropriately diluted methanolic extracts was used and a calibration standard curve with Trolox (2-12 µM; R 2 = 0.997 for DPPH assay and R 2 = 0.998 for ABTS assay) was made. Briefly, for ABTS assay, 3.9 mL of the ABTS radical solution were mixed with the sample. The spectrophotometric readings were carried out after 6 min at 734 nm. For DPPH assay, sample was mixed with 3.9 mL of DPPH 60 µM and stored in the dark for 120 min, and the absorbance was recorded at 515 nm. The results were expressed in terms of Trolox equivalent antioxidant capacity (TEAC), as mmol Trolox equivalents 100 g −1 dry weight of leaves (mmol TEAC 100 g −1 DW).

RP-HPLC Analysis of Phenolic Compounds
RP-HPLC method was used to determinate phenolic compounds by an Agilent 1260 series HPLC instrument (Agilent Technologies, Palo Alto, CA, USA) equipped with a quaternary pump (G1311B), degasser, column thermostat (G1316A), auto-sampler (G1329B) and diode array detector (G1315 B, DAD). Chromatographic separation was carried out according to Molinu et al. [20]. The column was a Zorbax Eclipse plus C18 (250 × 4.6 mm, 5 µm; Agilent). Column temperature was set to 30 • C and the flow rate was 0.8 mL min −1 . The injection volume was 10 µL and the detection wavelengths were set to 280, 330 and 350 nm. Data were processed using the Agilent OpenLAB CDS ChemStation edition 2012. Peak assignment of phenolic compounds was carried out comparing their retention times and UV-vis spectra with those of analytically pure standard (Table 7), and as well as by adding the standard solution to the sample. Table 7. Phenolic compounds identified in Parthenium genus and members of the Asteraceae family screened in leaf and twig extracts of guayule by HPLC-DAD.

Data Analyses
Laboratory measurement values (http://dx.doi.org/10.17632/9zsgwd4gp6.1) were subjected to a two-way analysis of variance, using Statgraphics Centurion XVI version [55] to test the effects of guayule accession, harvest time and their interaction on the following variables: concentrations for antioxidant capacity, total phenolics, total flavonoids and individual phenolic compounds. Differences between means were assessed with the Fisher's least significant difference (LSD) procedure for means separation. The significance level was fixed at p ≤ 0.05 for all the statistical analyses. The regression analyses between polyphenols and antioxidant capacity were calculated using Microsoft Excel 2016.

Conclusions
The current study is the first report regarding phenolic content, antioxidant activity, and quantification of individual phenolic compounds from biomass fractions of guayule plants grown in a Mediterranean environment.
Our research evidenced that the antioxidant capacity and the content of total bioactive compounds in leaves and twigs were affected significantly by guayule accession and harvest time. However, the improved lines AZ-1 and AZ-2 showed higher contents of total phenolic in twigs, irrespective of harvest time. The cinnamic acid derivatives namely neochlorogenic and chlorogenic acids in leaves and chlorogenic and 3,5-Di-caffeoylquinic acids in twigs were the main antioxidants isolated from guayule extracts.
Therefore, leaves and twigs, the latter at a lower extent, represent a rich source of bioactive compounds in guayule for potential exploitation in fodder, nutraceuticals and pharmacological sectors. Additional studies are needed to characterize biological activities of these extracts.


extfieldsofstudy:
Materials Science
Medicine

text:
Optimization of pulsed electric field conditions for sugar extraction from carrots.

Abstract Physical destruction and thermal treatment are pretreatment methods used to destroy cell membranes and facilitate the release of solute extraction. In this paper, sugar extraction from carrots under different pulsed electric field conditions (field strengths of 250, 750, and 1,250 V/cm, pulse numbers of 10, 45, and 80, and pulse frequency of 1 Hz) and simultaneous thermal treatments (at 20, 45, and 70°C) were studied based on full factorial design experiments with 27 runs. Carrot slices treated with PEF were suspended in water at the desired temperature and liquid‐to‐solid weight ratio of L/S = 2. Immediately after the PEF treatment, a significant increase in solute extraction was observed due to the permeability of cell membrane that could lead to the enhancement of solute convection on the surface of the tissue. Optimum extraction parameters were obtained as follows: PEF with the field intensity of 750 V/cm, 10 pulses, and temperature of 45°C.

| INTRODUC TI ON
Sugars exist in most of the plant tissues, but only sugar beet and sugar cane are significant sources of sugar for optimal extraction in the industry. In general, carrot with 4.7% of total sugar is placed after sugar beet, sugar cane, and sweet corn, respectively, and could be considered as a source of sugar (Dubounet, 2011).
Extraction as the first step of plant tissue study plays an important role in the final result and the following processes. The extraction process is influenced by factors such as matrix properties of the plant, solvent, temperature, pressure, and time (Azmir et al., 2013). The complete extraction of intracellular components in a solid-liquid extraction process involves the denaturation of cell membranes in the solid phase. In this way, sugar is extracted by countercurrent water at the high temperature of liquid phase within 70-74°C during the extracting times of 1-1.5 hr, which leads to the disintegration of cell membranes. The disadvantage of these processes is the decomposition of the inner chemical structure of cell walls by hydrolytic reactions. At such elevated temperatures, the amount of colloidal substances such as proteins, pectin, and other polluting materials entering the raw juice would be increased, resulting in low-purity raw juice which is very harmful for evaporation units (Lopez, Puertolas, Condon, Raso, & Alvarez, 2009;Sotudeh-Gharebagh, Shamekhi, Mostoufi, & Norouzi, 2009;Van der Poel, 1998). Nonconventional methods, which are more environmentally friendly, reduce operational time, and have better yield and quality of extract, have been developed during the last 50 years (Azmir et al., 2013). To enhance the overall yield and selectivity of components from plant materials, ultrasound (Karki et al., 2010;Mason, Paniwnyk, & Lorimer, 1996), pulsed electric field (A B Jemai & Vorobiev, 2006;Toepfl, Mathys, Heinz, & Knorr, 2006), enzyme digestion (Wang, Yang, & Wei, 2010), microwave heating (Li, Jia, Wei, & Liu, 2012), supercritical fluids (Cardenas-Toro et al., 2014), and accelerated solvents (Nayak et al., 2015) have been studied as nonconventional methods. PEF 1 is a nonthermal extraction method to improve the membrane that allows cold or hot aqueous penetration of valuable materials from the plant matrix. During recent years, several studies have been conducted about the effects of PEF on the extraction of soluble materials from plant tissues (El Belghiti & Vorobiev, 2004;Jemai & Vorobiev, 2003;El-Belghiti, Rabhi, & Vorobiev, 2005;Fincan, DeVito, & Dejmek, 2004;Jemai & Vorobiev, 2002;Corrales, Toepfl, Butz, Knorr, & Tauscher, 2008). When the sample is exposed to an electric field of sufficiently high intensity, the cell membrane can be irreversibly degraded . The destruction of a biological cell (plant or animal) with a high-intensity electric field in the form of very short pulses causes the formation of pores on the cell membrane that can be either temporary or permanent. This phenomenon which is known as electroporation enhances the permeability of the treated cells, eventually leading to the breakdown of cell membranes (Mandal, Mandal, & Das, 2015).
Applying the PEF processing at moderate electric field intensity (E < 100 V/cm) and room temperature, or moderate thermal treatment (T < 50°C) without any electric field, requires a long treatment time (PEF or thermal) and, subsequently, consumes a great deal of energy. Applying PEF and thermal treatment simultaneously has a strong effect on tissue damage (Lebovka, Praporscic, Ghnimi, & Vorobiev, 2005;Lebovka, Shynkaryk, El-Belghiti, Benjelloun, & Vorobiev, 2007). This effect can lead to structural transitions occurring inside membranes at high temperatures (Zimmermann, 1986). The kinetics of these processes depends on the parameters of PEF (field intensity, number of pulses, and pulse duration), external solute transfer conditions (stirring of solution), and the ratio of liquid-to-solid, and effects of sample size (El Belghiti & Vorobiev, 2004).
This work aims to apply full factorial design (FFD) methodology to study the simultaneous effect of PEF parameters (field intensity and the number of pulses) and temperature on sugar extraction yield and the disintegration index of carrot slices to obtain the optimal conditions of PEF and thermal treatments. Compared with other published findings, the data analysis and optimization process of these two factors resulted from the deficiency of a one-variable-at-a-time approach. The conventional method of maintaining other factors involved in a study at a constant level and studying each factor separately does not include the interaction effects of all the factors and is time-consuming as well. Therefore, a statistical experimental design minimizes the above-mentioned difficulties by studying all the significant factors simultaneously (Elibol, 2002). The factorial design methodology is an important tool for data analysis and finding optimal process conditions. Using this method leads to a more efficient experimental design that avoids misleading conclusions (Gheshlaghi, Scharer, Moo-Young, & Douglas, 2008).

| Sample preparation
Carrot samples used in this study were supplied from the special carrot farm with known species, in Mashhad, wrapped in sealed bags, and stored for 2-3 days in the refrigerator at 4 ± 2°C until experimental testing. The carrots were washed, and all the wastes were cleaned from the surface. Upper and lower parts of each carrot (approximately 2 cm) were removed. Then, the carrots were cut into disk pieces of 1 ± 0.1 cm in thickness, 3 ± 0.2 cm in diameter, and 10 ± 0.5 g in weight for the experiments. For each experiment, 10 pieces of carrot with the total weight of 100 ± 0.50 g were weighed by a balance (model FX3000, AND, Japan).

| PEF pretreatment
Carrot disk samples were placed inside a Plexiglas treatment chamber with dimensions of 4 × 10×10 cm 3 between two stainless steel electrodes with 4 cm of gap and area of 10 × 10 cm 2 . Then, tap water (0.05 S/m, 20°C) was added to the treatment chamber at the room temperature, so that the carrot slices would be immersed in water completely. The liquid-to-solid weight ratio was fixed at the value of 2 (L/S = 2). The PEF generator (1,750 V-20 A, Food Tech, Iran) could create electrical current of 1,750 V/cm with exponential decay wave pulses. Electrical current from an AC power supply (220-240 V 50 Hz) was transferred to a DC power supply. Then, electrical energy was transferred to a series of capacitances by a linear current, and the energy stored in the capacitors was discharged to the electrodes and treatment chamber by a pulse key. The energy of each pulse (Q) used in PEF treatment during exponential decay was given by Equation (1) (Maskooki & Eshtiaghi, 2012): where V is the peak value of decaying voltage (volt) in the sample and C is the capacity of the condenser (F). The temperature changes (ΔT) in the treatment chamber due to dissipation electrical energy in the chamber (nQ) during PEF treatment were calculated using the following equation: where c p and m are the specific heat and mass of samples, respectively.
The relation between cell disintegration index (Z) and energy input under different PEF treatment conditions was expressed as the equation below: where E eff is energy efficiency, Z cell disintegration index, and nQ is total energy input. This index showed effective cell disintegration per (2) mc p ΔT = nQ (3) E eff = Z∕nQ unit energy input by PEF treatment. This index showed effective cell disintegration per unit energy input by PEF treatment. PEF parameters included the electric field intensities of 250, 750, and 1,250 V/cm as well as different numbers of monopolar pulses of 10, 45, and 80. The pulse frequency was kept at constant 1 Hz. The juice elevated temperature was negligible (about 1°C).

| Thermal treatment
Immediately after PEF pretreatment, the samples were immersed in the juice that had previously reached the desired temperature (20-45-70°C) using a magnetic agitator (0-310°C, 0-1,500 rpm, IKA, RCT basic, Germany) according to the experimental design. In all the experiments, the solution was stirred at 250 rpm in a closed cylindrical device to prevent evaporation effects during the extraction process.

| Determination of extraction yield (c*)
To study the extraction kinetics, the solute concentration (Brix, g soluble 100 g −1 juice) was measured by a digital refractometer (RX-5000α, Atago, Japan) every 15 min up to 4 hr during the extraction process. The mathematical model used in this study was the one proposed by So and MacDonald (1986), which involved two first-order kinetics simultaneously. The first stage was washing which was related to the transfer of solute by convection and occured quickly; the second stage was diffusion related to the prolonged solute transfer by diffusion from the intracellular plant tissue matrix. This model can be described by Equation where c* can be described according to Equation (5): where c is the actual solute concentration in the solution, c ∞ is the equilibrium solute concentration, c ∞ ≅ c 0 ∕(n + 1), c 0 is the concentration of juice in fresh plant tissue in g[soluble]/100 g [juice], n is the liquid-tosolid ratio (L/S = 2), L is the extract volume in m 3 , S is the solid volume in m 3 , c * w is the c w /c ∞ ratio, c w is the equilibrium solute concentration due to the washing stage alone, c * d is the c d /c ∞ ratio, c d is the equilibrium solute concentration due to the diffusion stage alone, k w is the kinetic coefficient for the washing stage in min −1 , k d is the kinetic coefficient for the diffusion stage of extraction in min −1 , and t is time in min.

| Determination of disintegration index (Z)
At the end of each experiment, the electrical conductivity of each piece was determined using a conductivity meter (model IDSC004) before and after PEF and thermal treatment by measuring the electrical current passing through the sample connected to a 12 V DC source by the following equation (Sakr & Liu, 2014): where L is the gap between two electrodes in m, A is the cross-sectional area of carrot samples in the cell in m 2 , I is the alternative current passing through the material in A, and V is the voltage across the material in V that is constant. The degree of tissue damage was estimated from the electrical conductivity disintegration index Z (Angersbach, Heinz, & Knorr, 2002;Lebovka, Shynkaryk, & Vorobiev, 2007b, 2007c: where σ is the measured electrical conductivity in S/m, and the subscripts i and d refer to the conductivities of intact and completely damaged tissues, respectively. The maximally damaged carrot tissues σ d were obtained by measuring the electrical conductivity of freezethawed tissue (Lebovka, Shynkaryk, El-Belghiti, et al., 2007;Lebovka et al., 2007bLebovka et al., , 2007c).

| Experimental design and data analysis
To develop a mathematical model for extraction yield and disintegration index as dependent variables in terms of predictor variables (field intensity, number of pulses, and temperature), a sequential approach using the design of experiments was carried out.
Each variable was coded according to Equation (8): where x i is the coded value of variable X i , and X i,low and X i,high are the values of the variable at low and high levels, respectively (Gheshlaghi et al., 2008;Montgomery & Design, 2001).

| Modeling and optimization by full factorial design (FFD)
Full factorial design is a strong candidate in which all the possible replacements of all factors and possible interactions within factors are tested, and each test is replicated at least twice and the results are analyzed using ANOVA 2 (Gottipati & Mishra, 2010;Morrison, Block, Strickland, Collier, & Peterson, 2008). Thus, this method has greater precision in estimating the overall main factor effects and their interactions (Gottipati & Mishra, 2010).
If the combinations of k factors are applied at three levels, a factorial design will consist of 3 k experiments (Montgomery, 1982). In the present study, a three-factor three-level FFD with two replicates at each factor combination (3 3 runs) was carried out to develop a mathematical model for extraction yield and disintegration index.
The levels of the three variables for the FFD are indicated in Table 1.
The low and high levels of the factors were selected based on some preliminary experiments.
After modeling the responses by FFD, this method was developed for the optimization of variables affecting the two investigated responses.

| Statistical analysis
The data analysis and graphics as well as optimization process were performed using the experimental design software Design-Expert (version 6.0.4). The goodness of the fitted model was evaluated by the coefficient of determination, R 2 . Furthermore, an adequate precision statistic was used to measure the signal-to-noise ratio. A ratio of greater than 4 is desirable. The normality of the errors was checked by the normal probability plot of the residuals. To detect the correlation between the residuals, the plot of residuals in time sequence was checked so that the residuals are structureless in the correct model (Montgomery, 1982) (not shown). The significance of a fitted model and independent variables was determined using ANOVA with Fisher's test (P-value < 0.05).

| RE SULTS AND D ISCUSS I ON
Electric field intensity and the number of pulses are the main characteristics of a PEF protocol that influence extraction yield and disintegration index. The purpose of the first step of experiments (FFD) was to investigate the contribution of each factor in responses and determine their interaction effect as well. The second step was performed for finding the best conditions of the variables to achieve high performance of responses. The levels and ranges of the three factors for FFD are given in Table 1. This value was based on preliminary investigations. The experimental matrix for the 3 n factorial design (n factors, each run at three levels) for extraction yield and disintegration index is summarized in Table 2. The results were analyzed in Design-Expert software, and the main effects and interactions between the factors were determined.

| Solute extraction yields (c*)
The main effects and the interacting factors affecting the extraction yield were determined by ANOVA. For extraction yield, the main effects A and C and the interaction AC were of higher statistical significance, with a 95% confidence level. The model F-value of 29.10 indicated that the model was significant. Based on F-ratio and p-value, statistically insignificant factors were discarded. Afterward, the following 2FI 3 refined model was statistically considered in explaining the extraction yield: Comparing the relative weights of various terms in the model signified that temperature (C) had a very important effect on the extraction yield. To ensure an appropriate model, the R 2 value (0.9699) of the refined model was compared and deemed in agreement with the adjusted coefficient of determination, R 2 (adj) value (0.9565). Figure 1 shows the interaction effects of electric field intensity and temperature on the solute extraction yield of carrot after PEF treatments at 45 pulses. From Equation 9, it is evident that the number of pulses and its interactions with electric field intensity (AB), as well as temperature (BC), are not significant. According to Angersbach, Heinz, and Knorr (2000), when the pulse width which is in the microsecond range (for the cell in potato tissue, τ = 0.7 µs; in apple tissue, τ = 1.4 µs; and in plant cell culture suspension, τ = 0.45 µs) became several times larger than the specific disintegration time (τ), the critical breakdown transmembrane potential that is the main cause of the membrane destruction did not depend on pulse duration and number of pulses (Angersbach et al., 2000).
Therefore, there was no significant difference by varying the number of pulses between 10 and 80 at the 95% confidence level.
According to Figure 1, at 20°C and the electric field intensity of 250 V/cm, only 17% of solute materials were extracted. As the energy consumption increased during PEF pretreatment and electric field intensity increased from 250 V/cm to 1,250 V/cm, the permeable cells and, subsequently, yield extraction enhanced up to 750 V/cm. However, by increasing the field intensity to higher than 750 V/cm, no significant difference was observed in solute extraction yield at the 95% confidence level. It can be expected that this amount of energy is enough for the maximum permeability of the cells. Also, with increasing temperature in the range of 20 to 70°C, the solute extraction yield increased.

| Cell disintegration index
The cell disintegration index of each run was calculated, and the responding values of FFD were utilized for ANOVA. The F-value of 22.57 for the model implied that the model was significant. Moreover, the main effects of temperature and electric field, as well as the interaction effect, highly influenced the disintegration index. The p-value of <0.0001 showed that the model was significant and there was only 0.01% chance that this level of fit could occur due to noise. The normal probability plot of the studentized residuals validated the assumption of normally distributed errors with constant variance and zero mean.
The plots of residuals against predicted values, factor levels, and run orders showed a virtually constant variance over the variable ranges (not shown). Therefore, all the plots and statistics were acceptable, and it could be concluded that the model was adequate to explain the behavior of the system in terms of disintegration index, explaining approximately 98% (due to R 2 ) of the variability of Z.
The following 2FI model was accordingly chosen for the disintegration index: This function describes how the experimental variables and their interactions influence the disintegration index. All the main factors and (10)  then, decreased. For example, in the sugar beet tissue, the maximum breakdown was achieved after 60 pulses, and excess pulsations were unnecessary (Maskooki & Eshtiaghi, 2012). The mechanism of this phenomenon has been described in other studies (Sale & Hamilton, 1967;Zimmermann, 1986). Longer pulses are more effective, and their effect is more obvious at the room temperature and moderate electric fields (De Vito, Ferrari, Lebovka, Shynkaryk, & Vorobiev, 2008) By increasing the number of pulses in the range of 45-80 at 1,250 V/cm, the disintegration index Z first increased and, then, decreased, showing a limit for the increase in the values for electric fields. Bazhal, Lebovka, and Vorobiev (2003) investigated the effect of electric field intensity on the maximal disintegration index of vegetable tissues (Bazhal et al., 2003). With increasing the electric field and, subsequently, the energy input, the specific disintegration factor passed through a maximum point. Increasing the electric field strength was necessary for the membranes' damage because the secondary cell walls increased resistance in the membrane area. At 250 and 750 V/cm, no significant difference was observed in the disintegration index by increasing the number of pulses from 45 to 80. (250 V/cm), the disintegration index was decreased by increasing the temperature to higher than 45°C. According to Lebovka, Shynkaryk, El-Belghiti, et al. (2007), this effect can be related to the generation of small gas bubbles that screen the electrodes at high temperatures (Lebovka Shynkaryk, & Vorobiev, 2007b, 2007c.

F I G U R E 3
The interaction effect of electric field intensity and temperature on the disintegration index at (a) n = 10, (b) n = 45, and (c) n = 80 pulses after 255 min of thermal treatment Based on Lebovka et al. (2005), PEF treatment under moderate heating (<50°C) can help achieve a high degree of tissue damage at a moderate electric field (Lebovka et al., 2005). The enhancement of field intensity resulted in high degree of tissue damage. The complex dependences of disintegration index Z on electric field intensity E and temperature T are due to the existence of a wide range of cells with different geometries and sizes (Lebovka Shynkaryk, & Vorobiev, 2007b, 2007c. From Figure 3, by increasing the electric field intensity to over 750 V/cm at moderate temperatures in the range of 20-45°C, the increase in the disintegration index was not significant at the 95% confidence level. These results were in agreement with those of other studies conducted on the apple tissue (Bazhal et al., 2003).
With increasing the electric field of 1,200 to 1,500 V/cm in the apple tissue, no significant difference was observed in the disintegration index. According to Bazhal et al. (2003), disintegration index Z in small quantities of E (E < 300-500) has one or more saturation levels and may stabilize after an initial increase. If the treatment is continued, a further increase in σ and, subsequently, Z is observed after a long time. At higher values of E, Z has only one saturation level and reaches the saturation level after a short time. Figure 4(a-c) depicts the interaction effects of the number of pulses and temperature on the disintegration index. By increasing the temperature to over 45°C at 10 and 80 pulses, no significant difference was observed in disintegration index Z at the 95% confidence level. This effect, as already mentioned, was because Z reached the saturation level.
According to Figures 2 and 4, the disintegration index increased by increasing the number of pulses up to 45 pulses at the electric field intensity of 250 V/cm. In the higher electric field (>250 V/cm), there was no significant difference in the disintegration index by increasing the number of pulses. According to Bazhal et al. (2003), the optimal electric field E opt for cells with a secondary cell wall in carrot, potato, cucumber, and apple is equal to 200-400 V/cm (Bazhal et al., 2003). Consequently, the electric field intensity of 750 V/cm and above as well as increasing the number of pulses and, subsequently, the treatment time did not create any significant effects on disintegration index Z.

| Process optimization
After analyzing the effects of different variables and their interactions on the responses, it is important to find the best conditions to achieve an optimal response in the range of experiments. The experimental design software has the ability to find the best response. As mentioned before, the FFD model was used for optimization. Table 2 shows the experiments used in this method. The optimization of the responses was performed by a multiple-response method, known F I G U R E 4 The interaction effect of the number of pulses and temperature on the disintegration index at (a) 250 V/cm, (b) 750 V/ cm, and (c) 1,250 V/cm after 255 min of thermal treatment as desirability (D) function, to optimize different combinations of process variables such as field intensity (E), the number of pulses (n), and temperature (T). This process aimed to improve the conditions to increase the yield of solute (c*) and disintegration index (Z) by considering minimum energy consumption that was hidden in all of the three factors. To achieve maximum desirability, the three factors were set within the range. Then, c* and Z were set to maximum.
As already mentioned, in the electric fields of higher than 300-400 V/cm in carrot tissue, 10 pulses (treatment time = 10 s) were enough to achieve maximum yield extraction. The yield was increased by increasing the electric field up to 750 V/cm. Afterward, no significant difference was found. However, by increasing the number of pulses to more than 10, no significant difference was seen in c*. Increasing the field intensity alternatively resulted in larger amounts of Z at 10 pulses. However, when the temperature was higher than 45°C, no significant difference was observed in Z.
According to the two observed situations, the optimal conditions for both of the responses simultaneously included the electric field intensity of 750 V/cm, pulse number of 10, and temperature of 45°C.
In this case, the solute yield extraction and the disintegration index were 64% and 0.43, respectively.

| CON CLUS ION
In this study, the method of pulsed electric field was utilized for sugar extraction from carrot. Design-Expert software and FFD were employed for obtaining the interaction of factors (electric field intensity, number of pulses, and temperature) and to find optimal conditions. The sugar extraction yield and disintegration index were measured as responses. The importance levels of electric field strength and temperature were almost identical based on the responses obtained. As a result, by the simultaneous application of a moderate electric field to the tissue at an average temperature, the same yield that was obtained at high temperatures can be obtained without any electric field, thereby saving energy consumption. The optimal conditions for the extraction process including the highest extraction yield and disintegration index can also be achieved under these situations.

ACK N OWLED G M ENT
Authors would like to thank Dr Abdolmajid Maskooki although he is not with us now, but his invaluable research projects in PEF technology were our main reference in Iran. We would also like to thank Dr.
Morteza Mohammadi for his assistance with PEF processing of the samples and comments that greatly improved the manuscript.

CO N FLI C T O F I NTE R E S T
None declared.