extfieldsofstudy: Medicine Biology

s2fieldsofstudy:
Agricultural And Food Sciences
Biology

text:
Comparative Root Transcriptomics Provide Insights into Drought Adaptation Strategies in Chickpea (Cicer arietinum L.)

Drought adversely affects crop production across the globe. The root system immensely contributes to water management and the adaptability of plants to drought stress. In this study, drought-induced phenotypic and transcriptomic responses of two contrasting chickpea (Cicer arietinum L.) genotypes were compared at the vegetative, reproductive transition, and reproductive stages. At the vegetative stage, drought-tolerant genotype maintained higher root biomass, length, and surface area under drought stress as compared to sensitive genotype. However, at the reproductive stage, root length and surface area of tolerant genotype was lower but displayed higher root diameter than sensitive genotype. The shoot biomass of tolerant genotype was overall higher than the sensitive genotype under drought stress. RNA-seq analysis identified genotype- and developmental-stage specific differentially expressed genes (DEGs) in response to drought stress. At the vegetative stage, a total of 2161 and 1873 DEGs, and at reproductive stage 4109 and 3772 DEGs, were identified in the tolerant and sensitive genotypes, respectively. Gene ontology (GO) analysis revealed enrichment of biological categories related to cellular process, metabolic process, response to stimulus, response to abiotic stress, and response to hormones. Interestingly, the expression of stress-responsive transcription factors, kinases, ROS signaling and scavenging, transporters, root nodulation, and oxylipin biosynthesis genes were robustly upregulated in the tolerant genotype, possibly contributing to drought adaptation. Furthermore, activation/repression of hormone signaling and biosynthesis genes was observed. Overall, this study sheds new insights on drought tolerance mechanisms operating in roots with broader implications for chickpea improvement.

Introduction
Drought imposes a serious threat to crop growth and productivity [1]. Moreover, the intensity and frequency of drought periods are predicted to increase in the future climate [2]. Feeding the world's growing population under water limiting conditions will be a big challenge. To sustain agricultural patterns related to different functional categories were observed. We discuss the possible role of these genes in drought tolerance that can facilitate the development of drought-tolerant chickpea varieties in the future.

Phenotypic Responses
When exposed to drought stress, root-related traits such as dry weight, length, and surface area of the tolerant chickpea genotype (ICC8261) were significantly higher than the sensitive genotype (ICC283) at vegetative stage (VS) ( Figure 1A-C). In contrast, at the reproductive stage (RS), root length and surface area of the sensitive genotype were higher than the tolerant genotype ( Figure 1B,C). The average root diameter of the tolerant genotype was higher than the sensitive genotype at RS ( Figure 1D). Root volume and relative water content (RWC) were similar in both genotypes at all growth stages ( Figure 1E,F).
Int. J. Mol. Sci. 2020, 21, x FOR PEER REVIEW 3 of 20 genes in drought tolerance that can facilitate the development of drought-tolerant chickpea varieties in the future.

Phenotypic Responses
When exposed to drought stress, root-related traits such as dry weight, length, and surface area of the tolerant chickpea genotype (ICC8261) were significantly higher than the sensitive genotype (ICC283) at vegetative stage (VS) ( Figure 1A-C). In contrast, at the reproductive stage (RS), root length and surface area of the sensitive genotype were higher than the tolerant genotype ( Figure  1B,C). The average root diameter of the tolerant genotype was higher than the sensitive genotype at RS ( Figure 1D). Root volume and relative water content (RWC) were similar in both genotypes at all growth stages ( Figure 1E,F).  (SE). Physiological data were obtained from three independent biological replicates. Statistically significant differences between sensitive and tolerant genotypes obtained by one-way ANOVA at p < 0.05 are depicted by an asterisk (*). Furthermore, shoot dry biomass was significantly higher in the tolerant genotype under drought stress at all developmental stages ( Figure 1G). The chlorophyll content was significantly higher in the tolerant genotype during drought stress at RTS ( Figure 1H). The specific leaf area (SLA) was lower in the drought-tolerant genotype at RS ( Figure 1I). Overall, chickpea genotypes ICC8261 and ICC283 showed contrasting responses to drought stress, which is consistent with the previous findings suggesting that ICC8261 shows features of drought tolerance [14,16].

Transcriptional Responses and GO Analysis
To gain mechanistic insights into regulatory mechanisms underlying the drought stress response, we did RNA sequencing. The tolerant genotype showed a total of 2161 differentially expressed genes (DEGs) in response to drought stress at the vegetative stage, out of which 1214 DEGs Error bars represent standard errors (SE). Physiological data were obtained from three independent biological replicates. Statistically significant differences between sensitive and tolerant genotypes obtained by one-way ANOVA at p < 0.05 are depicted by an asterisk (*). Furthermore, shoot dry biomass was significantly higher in the tolerant genotype under drought stress at all developmental stages ( Figure 1G). The chlorophyll content was significantly higher in the tolerant genotype during drought stress at RTS ( Figure 1H). The specific leaf area (SLA) was lower in the drought-tolerant genotype at RS ( Figure 1I). Overall, chickpea genotypes ICC8261 and ICC283 showed contrasting responses to drought stress, which is consistent with the previous findings suggesting that ICC8261 shows features of drought tolerance [14,16].

Transcriptional Responses and GO Analysis
To gain mechanistic insights into regulatory mechanisms underlying the drought stress response, we did RNA sequencing. The tolerant genotype showed a total of 2161 differentially expressed genes (DEGs) in response to drought stress at the vegetative stage, out of which 1214 DEGs were upregulated  Figure 2A). Similarly, in the sensitive genotype, 1873 DEGs were obtained, out of which 821 DEGs were upregulated and 1052 DEGs were downregulated (Figure 2A). Venn diagram analysis showed that 261 DEGs were commonly upregulated and 262 DEGs were commonly downregulated in these genotypes ( Figure 2B) were upregulated and 947 DEGs were downregulated (Figure 2A). Similarly, in the sensitive genotype, 1873 DEGs were obtained, out of which 821 DEGs were upregulated and 1052 DEGs were downregulated ( Figure 2A). Venn diagram analysis showed that 261 DEGs were commonly upregulated and 262 DEGs were commonly downregulated in these genotypes ( Figure 2B). Interestingly, 61 DEGs showing upregulation in the tolerant genotype were downregulated in the sensitive genotype, and 20 DEGs showing downregulation in the tolerant genotype were upregulated in the sensitive genotype. At the reproductive stage, a total of 4109 DEGs were identified in the tolerant genotype in response to drought stress, out of which 2139 DEGs were upregulated and 1970 DEGs were downregulated ( Figure 2C). In the sensitive genotype, 3772 DEGs were identified wherein 2038 were up-and 1734 were downregulated. Venn diagram analysis showed that 1163 and 1059 DEGs were commonly up-and downregulated in both the genotypes, respectively ( Figure 2D). Out of these, 19 DEGs that were upregulated in the tolerant genotype were downregulated in the sensitive genotype, while 3 DEGs that showed downregulation in the tolerant genotype were upregulated in the sensitive genotype ( Figure 2D). The reliability of RNA-Seq gene expression was confirmed by qRT-PCR. The fold changes obtained from RNA-seq and qRT-PCR showed good correlation with each other, thus validating the RNA-seq data ( Figure S1).
Next, the functional over-representation of drought-responsive genes was performed by gene ontology (GO) analysis. GO terms classified genes into biological processes, molecular function, and cellular components (Figure 3). The most dominant terms recurring in the biological process category at vegetative and reproductive stages were 'cellular process', 'metabolic process', 'response to stimulus', 'response to abiotic stress', and 'response to hormones'. In the molecular function category, 'catalytic activity', 'binding', 'hydrolase', 'transporter activity', and 'antioxidant activity' terms were At the reproductive stage, a total of 4109 DEGs were identified in the tolerant genotype in response to drought stress, out of which 2139 DEGs were upregulated and 1970 DEGs were downregulated ( Figure 2C). In the sensitive genotype, 3772 DEGs were identified wherein 2038 were up-and 1734 were downregulated. Venn diagram analysis showed that 1163 and 1059 DEGs were commonly up-and downregulated in both the genotypes, respectively ( Figure 2D). Out of these, 19 DEGs that were upregulated in the tolerant genotype were downregulated in the sensitive genotype, while 3 DEGs that showed downregulation in the tolerant genotype were upregulated in the sensitive genotype ( Figure 2D). The reliability of RNA-Seq gene expression was confirmed by qRT-PCR. The fold changes obtained from RNA-seq and qRT-PCR showed good correlation with each other, thus validating the RNA-seq data ( Figure S1).
Next, the functional over-representation of drought-responsive genes was performed by gene ontology (GO) analysis. GO terms classified genes into biological processes, molecular function, and cellular components ( Figure 3). The most dominant terms recurring in the biological process category at vegetative and reproductive stages were 'cellular process', 'metabolic process', 'response to stimulus', 'response to abiotic stress', and 'response to hormones'. In the molecular function category, 'catalytic activity', 'binding', 'hydrolase', 'transporter activity', and 'antioxidant activity' terms were present. The cellular component category contained terms such as 'cell part', 'cytoplasm', 'intracellular', 'chloroplast', and 'cell wall'. Overall, these results suggest that roots underwent global transcriptional reprogramming during drought stress.

Root-Nodule Development Genes
Chickpea roots contain nodules that help in the establishment of a symbiotic association with nitrogen-fixing rhizobia. The expression of putative nodulation receptor kinase (LOC101507037) and calcium/calmodulin-dependent serine/threonine-protein kinase DMI3 (LOC101513751) were downregulated in the sensitive genotype both at VS and RS; however, the tolerant genotype showed marginal downregulation at RS only ( Figure 6, Table S3). In contrast, receptor-like kinases 3 (LOC101496137) and CLAVATA1 (LOC101488348) were upregulated at VS and RS in the tolerant genotype ( Figure 6).

Root-Nodule Development Genes
Chickpea roots contain nodules that help in the establishment of a symbiotic association with nitrogen-fixing rhizobia. The expression of putative nodulation receptor kinase (LOC101507037) and calcium/calmodulin-dependent serine/threonine-protein kinase DMI3 (LOC101513751) were downregulated in the sensitive genotype both at VS and RS; however, the tolerant genotype showed marginal downregulation at RS only ( Figure 6, Table S3). In contrast, receptor-like kinases 3 (LOC101496137) and CLAVATA1 (LOC101488348) were upregulated at VS and RS in the tolerant genotype ( Figure 6).

Discussion
Drought imposes a serious threat to agriculture. The root system plays a crucial role in the uptake of water and nutrients from the soil, thus helps in coping with drought stress. Root development and architecture are determined by intrinsic genetic properties and/or modulated by various environmental factors such as water availability [32,33]. So, we compared phenotypic and root-specific transcriptional responses of two contrasting chickpea genotypes at vegetative (VS), reproductive transition (RTS), and reproductive (RS) stages.
The tolerant genotype had higher root and shoot dry biomass under drought stress, possibly contributed by better water foraging capacity of roots and/or onset of resilience mechanisms. Indeed root length and surface area of tolerant genotype were higher at the vegetative stage; however, conservative growth was observed at the reproductive stage. In contrast, the root diameter showed the opposite trend. Longer roots help plants to forage water from deeper soils, and the larger surface area provides better opportunities for mycorrhizal colonization that facilitate nutrient acquisition [34]. The opposite response of root diameter could be explained by the fact that diameter controls the length and surface area for given biomass allocated to the root system, thus summarize the overall effect [35]. Chlorophyll content of the drought tolerant genotype was higher under drought stress at RTS, suggesting improved light-harvesting and energy production capacity to support growth and reproduction [36]. Also, lower specific leaf area in the tolerant genotype can improve water-use efficiency, as reported previously [37,38].
RNA sequencing allows genome-scale quantification of transcriptome changes underlying developmental transitions and stress responses [39]. In the present work, genes related to various biological functions were altered in response to drought stress. Transcription factor (TF) genes related to ABA-dependent (ABI5) and ABA-independent (DREB1A and DREB1C) pathways were upregulated in the tolerant genotype at RS, indicating robust induction of multiple regulatory networks to impart drought tolerance. Also, genotype-specific alterations in other TF genes such as AP2/ERF, bHLH, bZIP, DOF, and WRKY were observed. The AP2/ERF forms a large group of plant-specific TFs induced by multiple stresses and phytohormones [40,41]. ERFs were previously shown to be induced by wounding and biotic stresses [42]. Likewise, WRKY family TFs have been suggested to play a role in adaptation to biotic and abiotic stresses, and modulate stress-responsive signaling and ROS production [43]. WRKY33 is involved in the activation of peroxidases and glutathione S-transferases under drought stress, playing a role in scavenging ROS [44]. We found that bHLH family TF genes (VIP1, bHLH18-like, and bHLH93-like) were upregulated in the tolerant genotype. These TFs may contribute to improved defense; for instance, VIP1 activates stress-responsive genes via the mitogen-activated protein kinase pathway [45]. Also, other TF genes such as C2C2-Dof, MYB, BTB/POZ, and TAZ were upregulated in the tolerant genotype. Overall, upregulation of TFs in the tolerant chickpea genotype may activate stress-responsive genes conferring drought tolerance.
Protein kinases are integral components of signal transduction pathways playing an important role in the perception and activation of stress-responsive pathways [46]. The MKK1 gene was upregulated in tolerant genotype both at VS and RS upon drought exposure. It may contribute to the maintenance of cellular homeostasis by lowering ROS levels under stress [47,48]. CDPK family members (CDPK-SK5, CDPK3, and CDPK4) were upregulated in both the genotypes. CDPKs are activated by binding of intracellular Ca 2+ with a calmodulin-like domain that eventually regulates the downstream targets [49]. For instance, CDPK4 phosphorylates NADPH oxidase to regulate ROS production during stress [50][51][52]. We found upregulation of RLKs (IMK2, CRK25, and PERK4) in the tolerant genotype at VS, while PERK14 and HSL2 showed up-regulation in both genotypes at RS. RLKs get rapidly activated during the early stages of drought stress, as reported previously [53,54]. PERK4 modulates root tip growth via ABA signalling, and HSL2 regulates lateral root emergence by controlling cell separation [53,54]. The WAK-like 20 was upregulated in the tolerant genotype at VS, which possibly facilitate communication between the cell wall and cytoplasm, and regulate stress responses [55][56][57][58][59]. Therefore, their enhanced expression in the tolerant genotype may contribute to cell wall loosening that minimizes root growth inhibition under drought. Also, CIPK1-like was upregulated in the tolerant genotype at RS. The CIPK1 functions as a regulator of ABA-dependent and independent signaling pathways by interacting with calcium sensors [60]. Taken together, early and robust upregulation of protein kinases possibly contributes to drought tolerance by activation of downstream stress signaling pathways.
ROS act as both signaling and damaging molecules whose levels are determined by antioxidant molecules and enzymes [47]. The expression of many antioxidant enzymes, such as glutathione transferases (GST) and peroxidases (PER), was upregulated in the tolerant genotype. Simultaneous upregulation of GST and PER indicate that multiple ROS detoxification components efficiently scavenge ROS under drought; however, their downregulation in the sensitive genotype may lead to excessive ROS accumulation, causing root growth inhibition. The expression of RBOHH and RBOHD were downregulated in the tolerant genotype at VS, whereas RBOHE and RBOHA were upregulated at RS. Previously, downregulation of RBOH family genes has been reported under drought stress [61,62]. The down-regulation of RBOHD and RBOHH may cause inhibition of ABA-dependent ROS signaling, which allows marginal root growth inhibition in the tolerant genotype. Thus, the synergistic action of a battery of antioxidant enzymes confers stress tolerance in plants [63].
Stress-induced signaling molecules are transported via active or passive transporters. Active transport is mediated by ATP dependent transporters such as ATP-binding cassette (ABC), whereas passive transport occurs through ion-channels and carriers such as aquaporins [64][65][66][67].
ABC transporters were upregulated in the tolerant genotype at VS (ABCB15-like, ABCB19, and ABCC12-like) and RS (ABCB15-like, ABCB29, and ABCG22-like). Upregulation of ABC transporters can be linked to multiple functions, which include auxin transport, suberin formation, ABA signaling and transport, pollen development, leaf water retention, and stress tolerance [64][65][66][67]. Furthermore, sugar (N3, SWEET1-like, SWEET4) and nitrate (NRT1/PTR FAMILY) transporters were upregulated in both genotypes. The SWEET11 facilitates sugar supply to roots in a cell-type-specific manner and participates in the establishment of symbiotic association [68]. Furthermore, SWEET12 and SUC2 functions are associated with the reallocation of carbohydrates to roots for efficient root development under drought [69]. Also, NRT transporters are involved in drought stress responses by controlling ABA transport [70]. Aquaporins passively transport water molecules, and we found upregulation of PIP2-1-like and TIP2-2 in the tolerant genotype. Consistent with this, overexpression of OsPIP1-3 and SlTIP2-2 led to the maintenance of leaf water potential under drought stress in rice and tomato, respectively [71,72]. Also, ABA regulates aquaporin activity and hydraulic conductivity to maintain favorable plant water status [73]. Therefore, aquaporins could serve as ideal targets for engineering drought tolerance in chickpea.
Leguminous plants harbor nitrogen-fixing rhizobia in root-nodules, which facilitates the uptake of nutrients and water. It has been observed that Medicago with nodulated roots shows delayed leaf senescence as compared to non-nodulated plants under drought stress [74]. We also found downregulation of nodulation receptor kinase and calcium-dependent kinase DMI3 in the sensitive genotype. DMI3 acts upstream of nodulation-signaling pathway proteins (NSP1/NSP2) and is a regulator of Nod-factor genes [75,76]. Thus, their down-regulation in the sensitive genotype suggests a decrease in root nodulation, resulting in low nitrogen assimilation. In contrast, LYK3 and CLAVATA1 were upregulated in the tolerant genotype. LYK3 regulates rhizobial infection via Nod-factor genes [77], and CLAVATA1 regulates root length and nodulation [78]. Therefore, their upregulation may contribute to root nodulation and nitrogen assimilation under drought stress.
Cross-regulatory hormonal network modulates plant growth and stress response in plants [79]. Auxin regulates root and shoot growth under optimal and stressful conditions [80]. In this study, auxin biosynthesis gene (YUCCA2) was downregulated in the sensitive genotype at VS, indicating attenuation of auxin biosynthesis. Consistent with this inhibition of root biomass, length, and surface area were observed in the sensitive genotype. In contrast, PIN genes (PIN1b, PIN1c-like, and PIN4) were downregulated in the tolerant genotype at RS. PIN proteins act as auxin efflux carriers for basipetal auxin transport, which in turn controls lateral root formation [81]. Their downregulation may be associated with the reduction of root length and root surface area. Likewise, cytokinin controls lateral root, apical meristem, and vascular system development [82,83]. Cytokinin biosynthesis (IPT5) and signaling (AHK5) genes were downregulated in the tolerant genotype. It has been shown that AHK5 inhibits root elongation via ETR1-dependent ABA and ethylene signaling pathways [84]. Also, earlier transcriptome studies highlight the dual role of ethylene to regulate growth and defense responses during stress [85]. We found upregulation of ethylene-responsive transcription factors (ERF) such as ERF1-like, ERF113-like, ERF13-like, and ERF2-like in the tolerant genotype. ERFs integrate ethylene and jasmonic acid signaling pathways [86] and participate in stress responses [87,88]. Ethylene is known to control lateral root formation via modulating auxin biosynthesis and transport [89], and root and shoot growth via cross-talk with ABA which is a central mediator of drought stress-induced signaling [90]. The expression of ABA hydroxylases (CYP707A3-like and CYP707A1-like) was upregulated in the tolerant genotype, which could catalyze oxidative degradation of ABA, maintaining optimal ABA levels during drought [91]. The PYL genes were downregulated in tolerant genotype, while phosphatases (PP2C) were upregulated in the sensitive genotype. ABA binds with PYL/PYR/RCAR receptors to inhibit the activity of PP2Cs that negatively regulate ABA signaling through repression of SnRK2. PP2C inhibition leads to SnRK2 de-repression, eventually phosphorylating and activating downstream transcription factors [92]. The PYL4 downregulation and CYP707A3 activation in drought-tolerant chickpea genotypes suggest inhibition of ABA signaling at VS, whereas the up-regulation of NCED1 and PYL4 indicate activation of ABA signaling at RS. Also, ABA accumulation is known to inhibit root growth by promoting ethylene biosynthesis [93]. Gibberellic acid induces cell proliferation and elongation growth in plants [94,95]. DELLA proteins are negative regulators of GA signaling that act downstream of GA receptors. The binding of DELLA to GA receptor GID1 leads to their degradation and activation of GA function [96]. We observed downregulation GA receptor (GID1B-like) and DELLA genes in both the genotypes. This indicates that repression of GID1 and DELLA might result in reduced GA responses, thus attenuating root growth under drought stress. The role of lipid-derived phytohormone Jasmonate in abiotic stress tolerance is emerging [97]. We found that expression of JA biosynthesis genes such as LOX, LOX homologs, and TIFY/JAZ were significantly upregulated in the tolerant genotype, suggesting enhanced production of JA in this genotype. Previously, it was reported that activation of allene oxide synthase (AOS), an enzyme involved in JA biosynthesis, enhances drought tolerance in chickpea [97]. In addition, JA and ABA pathways are involved in the regulation of stomatal opening and closing to improve drought tolerance [98]. Overall, complex hormonal-crosstalk controls root growth under drought stress in chickpea.

Plant Material and Experimental Set-Up
The seeds of two chickpea (Cicer arietinum L.) cultivars viz. ICC283 (Desi type) and ICC8261 (Kabuli type), exhibiting contrasting responses to drought stress, were procured from ICRISAT (International Crops Research Institute for the Semi-Arid Tropics) Hyderabad, India. The genotype ICC283 is drought-sensitive, while ICC8261 is a drought-tolerant genotype. This phenotypic divergence in response to drought is attributed to their root system architecture, where the susceptible ICC283 genotype possesses shallow root systems, while the tolerant ICC8261 genotype possesses prolific root systems [14,16]. The seeds of these genotypes were sown in 27-cm polypropylene pots containing 9.5 kg of soil. Pots were placed in controlled growth conditions of the day (28 • C) and night (23 • C) temperatures, and randomly assigned as well-watered (WW) and drought-stressed (DS). The optimal water level of 95% ASWF (available soil water fraction) was maintained in WW pots by irrigating pots on alternate days. DS pots were maintained at 95% ASWF at sowing and were brought to 70% ASWF by dry-down fifteen days before the sampling. The pots were weighed every day and water was added to compensate the water lost through transpiration. The experiment was performed in a 3 × 2 × 2 (three time points, two genotypes, and two treatment conditions) completely randomized block design. Root tissue was harvested for physiological and RNA-seq analyses from WW and DS treatments at 30, 50, and 70 days after sowing (DAS). These time points reflect the vegetative stage (VS; 30 DAS), the reproductive transition stage (RTS; 50 DAS), and the reproductive stage (RS; 70 DAS). Root tissue harvested for RNA-seq was immediately put in liquid nitrogen. Three biological replicates were used to carry out physiological and transcriptional analyses.

Phenotypic Analysis
Phenotypic root and shoot traits were quantified in the sensitive and tolerant chickpea genotypes at three developmental stages viz. VS (30 DAS), RTS (50 DAS), and RS (70 DAS). Dry root and shoot biomass were determined after drying plant material at 65 • C for 3 days in the oven. To determine root-related traits such as length, surface area, diameter and volume harvested roots were washed with water and dried between the folds of filter paper. Then images were captured to compute these parameters. Relative water content (RWC) was calculated as RWC = [(fresh weight − dry weight)/(turgid weight − dry weight)] × 100. Leaf-related traits such as chlorophyll content (SPAD chlorophyll meter) and specific leaf area (SLA) was calculated as SLA = leaf area (cm 2 )/leaf dry weight (g). All the phenotypic analyses were performed in at least three replications, and data were reported as mean and standard error (SE). Means were statistically compared by one-way ANOVA at the probability level of 5% (p < 0.05) by using SPSS software (IBM, New York, NY, USA).

RNA Extraction and Library Construction
Total RNA was extracted with an RNeasy Mini Kit (Qiagen, Hilden, Germany). RNA quantification was done by using a Nanodrop Lite spectrophotometer (Thermo Fisher Scientific, Wilmington, MA, USA) and the integrity of RNA was checked with Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). For library preparation and sequencing, RNA samples with a 260/280 absorbance ratio between 2.1 to 2.2, 260/230 absorbance ratio between 2.0 to 2.5, and RIN (RNA integrity number) of more than 7.0 were used.
The enrichment of mRNA samples was done using a MicroPoly(A)Purist Kit (Thermo Fisher Scientific, Waltham, MA, USA) by following the guidelines of the manufacturer. The RNA-seq libraries were prepared using Ion Total RNA-Seq Kit v2 (Thermo Fisher Scientific, Waltham, MA, USA) and sequenced by using an Ion Proton System (Thermo Fisher Scientific, Waltham, MA, USA). The low-quality reads and primer/adapter sequences (QUADTrim; https://bitbucket.org/arobinson/ quadtrim) were removed to obtain single-end reads of length varying from 50-265 bp with an average of~80 bp. We sequenced thirty-three libraries in total, each genotype represented by two treatments, three time-points, and three biological replicates. The three samples of tolerant genotype (ICC8261) at the reproductive transition stage (RTS) were not sequenced due to RNA degradation.

Transcriptome Analysis
Sequenced reads were mapped to the Chickpea Reference Genome (Cicer arietinum v1.0) [28] by using Bowtie v2.2.3 [99] and Tophat v2.1.1 [100] with default settings. Then mapped reads were counted with HTSeq-count v0.61 [101], and differential expression analysis was performed by the EdgeR package [102]. This package runs in the R statistical environment, uses a count-based approach, and employs the over-dispersed Poisson model to account for both biological and technical variability. Finally, differentially expressed genes (DEGs) were selected by applying a cut-off of p-value ≤ 0.001, FDR ≤ 0.05, log2 fold change ≥ +1.0 and ≤ −1.0.


extfieldsofstudy: Chemistry Medicine

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

text:
Enzymatic conversion reactions of 5-hydroxymethylfurfural (HMF) to bio-based 2,5-diformylfuran (DFF) and 2,5-furandicarboxylic acid (FDCA) with air: mechanisms, pathways and synthesis selectivity

Background 2,5-Furandicarboxylic acid (FDCA) is one of the top biomass-derived value-added chemicals. It can be produced from fructose and other C6 sugars via formation of 5-hydroxymethilfurfural (HMF) intermediate. Most of the chemical methods for FDCA production require harsh conditions, thus as an environmentally friendly alternative, an enzymatic conversion process can be applied. Results Commercially available horseradish peroxidase (HRP) and lignin peroxidase (LPO), alcohol (AO) and galactose oxidase (GO), catalase (CAT) and laccase (LAC) were tested against HMF, 2,5-diformylfuran (DFF), 5-hydroxymethyl-2-furoic acid (HMFA) and 5-formyl-2-furoic acid (FFA). Enzyme concentrations were determined based on the number of available active sites and reactions performed at atmospheric oxygen pressure. AO, GO, HRP and LPO were active against HMF, where LPO and HRP produced 0.6 and 0.7% of HMFA, and GO and AO produced 25.5 and 5.1% DFF, respectively. Most of the enzymes had only mild (3.2% yield or less) or no activity against DFF, HMFA and FFA, with only AO having a slightly higher activity against FFA with an FDCA yield of 11.6%. An effect of substrate concentration was measured only for AO, where 20 mM HMF resulted in 19.5% DFF and 5 mM HMF in 39.9% DFF, with a Km value of 14 mM. Some multi-enzyme reactions were also tested and the combination of AO and CAT proved most effective in converting over 97% HMF to DFF in 72 h. Conclusions Our study aimed at understanding the mechanism of conversion of bio-based HMF to FDCA by different selected enzymes. By understanding the reaction pathway, as well as substrate specificity and the effect of substrate concentration, we would be able to better optimize this process and obtain the best product yields in the future.

inexpensive potential source of new, greener chemicals [2]. Therefore, to switch from petroleum-based to biomass-based chemicals, new processes and technologies have to be developed [3].
One such important platform chemical is 2,5-furandicarboxylic acid (FDCA). FDCA and its derivates can be applied in many fields, but the most promising use of this chemical is in the replacement of tetraphthalate in oil-based plastics like polyethylene terephthalates (PET) [4,5]. FDCA is mostly produced from fructose and other sugars via 5-hydroxymethylfurfural (HMF) as an intermediate [6]. There are two routes of HMF oxidation to FDCA ( Fig. 1): (1) the aldehyde group of HMF is oxidized to form a carboxylic acid, thus yielding 5-hydroxymethyl-2-furoic acid (HMFA); then HMF acid is oxidized to yield 5-formyl-2-furoic acid (FFA) and FDCA, and (2) the alcohol group of HMF is oxidized to yield the 2,5-diformylfuran (DFF), followed by further oxidations to FFA and FDCA. Most chemical methods for production of FDCA from HMF require harsh conditions like organic solvents, high temperature and pressure and special salts. This requires high energy expenditure as well as causes environmental pollution [7][8][9][10].
An alternative to chemical synthesis of FDCA by harsh oxidation methods is biocatalysis using either whole-cell [15][16][17][18][19] or enzymatic [11,13,20] conversion processes. These are usually mild and environmentally friendly since they are carried out at lower temperatures and pressures while producing less toxic waste [21]. The whole-cell conversion process has proved to be efficient in production of FDCA [15][16][17][18][19] and does not require isolation of enzymes which presents an extra cost to the process. However, while these microbes produce and secrete enzymes that catalyse the oxidation of HMF to FDCA, other compounds can also be present in the reaction mixture. These are either secreted by the microbes or are a part of the growth medium for the microbe and can subsequently affect the activity of the enzyme and thus the outcome of the reaction. So in order to determine Fig. 1 Schematic representation of oxidation routes of HMF to FDCA by different enzymes (E). Examples of enzymes for E1, E2 and E3 were taken from literature [11][12][13][14]. AAO aryl-alcohol oxidase, HMFO HMF oxidase, UPO unspecific fungal peroxidase, Lac laccase the mechanism of action of a specific enzyme, isolated enzymes have to be used.
However, only a few enzymes have been shown to be active toward HMF and/or its oxidation products, DFF, HMFA and FFA. The conversions are mostly slow and incomplete and with low yields of FDCA. Also, most of these enzymes are not capable of full conversion of HMF to FDCA and thus need to be combined in multi-enzyme reactions [11,15,17,20,22,23]. For example, an arylalcohol oxidase was able to oxidase HMF to FFA with 98% yield after 4 h but only 6% FDCA after 24 h [11]. An unspecific fungal peroxygenase oxidized HMF to up to 97% HMFA after 24 h, but only 10% FDCA. 90% FDCA was obtained with oxidation of FFA after 96 h with the same enzyme [11]. A fungal laccase, with TEMPO as a mediator, was able to produce 90.2% FDCA and 8.2% FFA after 96 h [12]. The only enzyme that was capable of almost complete conversion of HMF to FDCA was HMF oxidase, an FAD-dependent oxidase that was able to oxidase HMF to 92% FFA and 8% FDCA in 5 h, and at higher enzyme concentrations to 95% FDCA after 24 h [13,14].
As can be seen from previous research, some fast and complete conversions have already been achieved and an insight into the mechanism of action of these selected enzymes has been provided. However, not all enzymes work the same way, even if they belong to the same group [20,22]. Thus, it is evident that in order to determine the most efficient enzyme or enzyme combinations for production of FDCA, understanding the mechanism of action for each individual enzyme is crucial. In our current work, we selected six different commercially available enzymes: alcohol oxidase, galactose oxidase, horseradish peroxidase, laccase, catalase and lignin peroxidase. Although HMF can be produced from fructose and other sugars [6], we used a commercially available analytical standard grade HMF to avoid any possible contaminants that could arise from its isolation from biomass. We determined their substrate specificity as well as the HMF oxidation route they employ and the effect of substrate concentration and/or cofactor on enzyme activity. Oxygen content in solution was measured during the reaction to determine the need for its supplementation. Finally, we designed some simple one-pot multi-enzyme reactions with the aim of increasing the productivity of selected enzymes. Some of the enzymes are composed of multiple units (multiple active sites), so in order to obtain and compare the activity per one active site, the enzyme concentrations were calculated based on the number of subunits of a specific enzyme. The aim of our study was not to improve the conversion process of HMF to FDCA, as that has already been achieved in previous studies, but to better understand the mechanism of conversion of HMF to FDCA by our selected enzymes: substrate selectivity, HMF oxidation route, effect of substrate concentration and effect of cofactor concentration.

Oxidation of HMF with alcohol oxidase
Alcohol oxidase (AO) with 1 µM flavin adenine dinucleotide (FAD) was first tested against 10 mM HMF. Samples were taken after 24, 48 and 72 h. The results in Fig. 2a show a steep rise in DFF content after 24 h with some FFA also being formed. There was an additional increase of DFF after 48 h and a slight decrease after 72 h. The FFA content slowly increased throughout the whole incubation period, but did not reach more than a few percent.
During AO substrate oxidation H 2 O 2 is formed which has an inhibitory effect on the enzyme [24]. Thus, the decrease in conversion of HMF to DFF after 48 h can be explained by decreased enzyme activity due to H 2 O 2 accumulation. A small part of the decrease in DFF content could also be assigned to its conversion to FFA.
Previous research has shown that H 2 O 2 produced by the action of an alcohol oxidase can also oxidase HMF and/or its oxidation products [11]. Therefore, we incubated HMF with H 2 O 2 and analysed the oxidation products of the reaction after 72 h. There was only minimal (around 0.1%) conversion of HMF to DFF (data not shown), and thus, we could conclude that the oxidation of HMF was executed by the enzyme and not H 2 O 2 .
According to literature, the enzymes use molecular oxygen as a terminal electron acceptor [13,24]. To ensure that the oxygen concentration in solution remains stable or does not drop too low to affect the activity of the enzymes, Quin et al. [22] used daily air bubbling. We measured oxygen concentration at the beginning of the experiment and then every 24 h. Results in Table 1 show that despite the enzyme activity, which was presented as conversion of HMF to DFF, the oxygen concentration remained stable during the whole process.
AO is an FAD-dependent enzyme and previous research has shown that FAD concentration can have a notable effect on the activity of this type of enzyme [13]. There they determined that the best enzyme-tocofactor molar ratio was 1:1 and a tenfold increase in FAD had a notable inhibitory effect. Therefore, we used different FAD concentrations with AO and tested them against HMF. The results in Fig. 3 show that in the case of AO, the concentration of FAD does not have an effect on the activity of the enzyme. Thus, all of the following experiments with AO were carried out with 1 µM FAD.

Enzymatic conversion of HMF, DFF, HMFA and FFA by other enzymes
To determine substrate specificity and productivity, as well as which HMF oxidation route was employed, we tested each enzyme against HMF and its oxidation products, DFF, HMFA and FFA. Enzyme concentrations used were based on the number of their subunits and thus the measured activity was the activity of one enzyme active site. The products and product yields after 72 h of incubation are presented in Table 2 and the time course of oxidation of HMF, DFF and FFA with AO is presented in Fig. 2a-c (oxidation of HMFA is not included as no product was formed). The highest activity against HMF was measured for AO where 25.5% of DFF and 3.1% of FFA was formed. The second highest activity was measured for GAO, but with only 5.1% of DFF and no FFA.  LPO and HRP showed some low activity, but instead of DFF, HMFA was formed. CAT and LAC were not active against HMF. The activity of all enzymes against DFF was low, with the highest of 3.2% of FFA being obtained with AO. The reaction was also reversible and some HMF was formed. GAO, LAC, LPO and HRP showed some activity against HMFA forming FFA and FDCA. The highest yield of FDCA was 4% and was obtained with HRP. All of the enzymes showed activity against FFA. All but one produced only small amounts of FDCA, between 0.6 and 3.2%, whereas AO produced 11.6%. The reaction was reversible in some cases. CAT, LAC and LPO produced only small amounts of HMFA and AO produced a relatively large amount of 18.2%. There was also some spontaneous conversion of substrates in the absence of an enzyme. HMF converted to 0.3% DFF, DFF to 1% FFA and 0.5% HMF and FFA to 0.5% FDCA. This spontaneous conversion shows that the already low activity of some enzymes is in fact even lower or non-existent.
In previous research only a few enzymes have been tested against HMF and even in those cases many were not tested against other oxidation products (DFF, HMFA and FFA) and/or they were tested in combination with others in multi-enzyme reactions. A similar enzyme to our AO, a fungal aryl-alcohol oxidase (AAO), was tested against HMF, DFF, HMFA and FFA. HMF and DFF both presented as good substrates for this enzyme yielding 90 or more percent of FFA after 4 h and also some FDCA after 24 h [11]. In our case, only 25.6% DFF, 3.1% FFA and no FDCA was produced after 72 h. DFF was shown to be a good substrate for AAO, because it readily forms a hydrate at neutral pH [11,25]. However, that has proven not to be the case with our AO since only low quantities of FFA and no FDCA have been formed. As in the case of AAO [11], AO also showed no activity against HMFA. Because FFA does not readily form a hydrate, AAO was not active against this substrate [11], however, our AO was. On the other hand, another study showed that certain AAOs did have some activity against FFA [20] which was similar to our AO.
A recombinant galactose oxidase was previously tested against HMF and produced DFF with no FFA. This GAO was also capable of oxidizing HMFA, but not DFF, to FFA [20]. The oxidation of HMF with our GAO also produced only DFF, however, it was also active against DFF, producing some FFA, and against HMFA, producing not only FFA, but also small amounts of FDCA. Another group used the same GAO used in our study and tested it against HMF with similar results, 2% of DFF after 72 h. Like AO, GAO also produces H 2 O 2 as a by-product which inhibits its action [22], thus removing this byproduct might increase the enzymes activity.
To our knowledge, only one group tested CAT and HRP in a single-enzyme reaction and even then they were only tested against FFA where they produced around 55% or 22% FDCA after 48 h [20]. CAT and HRP were previously also tested in tandem against HMF, but the reaction yielded no products [22]. Both of these enzymes were mainly used in multi-enzyme reactions with AAO, AO or GAO for removal of the accumulated H 2 O 2 [22,23].
In previous research laccases have been used either in free form [22] or immobilized [12]. In both cases, TEMPO was used as a mediator. In the case of free enzyme reactions they tested three different LACs against HMF and obtained DFF, FFA and FDCA in yields of 68-82%, 4-6% and 5-10% after 48 to 96 h, respectively [22]. In the case of immobilized enzyme, the conversion of HMF to FDCA was over 90% in 96 h, showing a stabilizing effect of immobilization on the enzyme. Our results show a very low activity of this enzyme against DFF, HMFA and FFA and no activity against HMF, thus, emphasizing the need for a mediator.
LPO was chosen for our research because, like the unspecific peroxygenases (UPOs) used in previous research [11,20], it is a heme-peroxidase [26]. These heme-peroxidases have been shown to catalyse a H 2 O 2 -dependent hydroxylation of alcohols into aldehydes and carboxylic acids [27,28]. Therefore, like in the case of CAT and HRP, UPO was also used in multienzyme reactions for removal of the produced H 2 O 2 [11,20]. When HMF was oxidized by UPO, HMFA instead of DFF was formed with a yield of 72% after 72 h [11].
Our results with LPO also show this enzyme employs the HMFA oxidation route of HMF, however, the yields were much lower. Also, the reaction with UPO continued on to form FFA and FDCA, whereas the reaction with LPO did not. Despite UPO being a H 2 O 2 -dependent enzyme, previous research has also shown that it is active even in the absence of H 2 O 2 and that the addition of this cofactor only mildly increases the production of FDCA from FFA. The increase of FDCA production seemed to be due to FFA oxidation with H 2 O 2 , which was further supported by results that showed that adding CAT, which degrades H 2 O 2 , also decreased the production of FDCA [20]. Thus, our reactions with LPO were performed in the absence of exogenous H 2 O 2 and showed mild activity against all of the tested substrates.
The differences in our results compared to those presented in literature can not only be explained by differences in enzymes used, but also in different enzyme and substrate concentrations and their ratios as well as reaction conditions. It has been shown that the same types of enzymes from different origins have different activities and different pH optimums [20,22]. The change in pH of the reaction solution can not only affect the activity of the enzyme, but also the HMF oxidation route it employs (DFF or HMFA) [20]. Also, some of the enzymes used in these studies are composed of multiple units (AOoctamer, CAT-tetramer) and thus have multiple active sites with which to oxidize a substrate. To our knowledge, previous research not only compared different enzymes at different concentrations, but it also did not take into account the different number of active sites on each enzyme. In our current study, we determined the enzyme concentrations based on the number of subunits (active sites) which enabled us to compare enzyme activities based on a single active site.

Effect of substrate concentration on enzyme activity
To determine if the substrate concentration affects enzymatic activity, we tested the enzymes against 20, 10 or Table 2 Oxidation Table 3 show that substrate concentration had the most notable effect in the case of AO and HMF. At 20 mM HMF, the conversion to DFF and FFA was 19.8% and 3.7%, respectively. At lower substrate concentrations the conversion to FFA decreased to 3%; however, the production of DFF markedly increased and reached almost 40% in the case of 5 mM HMF. In the case of other enzymes, the conversion rates remained similar and thus showed no effect of substrate concentration on enzyme activity. The effect of substrate concentration on enzymatic activity was previously tested on a two-enzyme reaction with an aldehyde oxidase PaoABC and a galactose oxidase GOase M3-5 [24]. At 10 mM HMF, an almost complete conversion to FDCA took place, with DFF and FFA intermediates being formed. By increasing the concentration of HMF to 20 mM, the production of FDCA was reduced by half along with formation of 50% of HMFA. Another multi-enzyme reaction, a combination of galactose oxidase, HRP and PaoABC, was also tested [23]. There they showed a similar trend where high HMF concentrations resulted in a reduced production of FDCA and an increased production of HMFA. This was countered by reducing the concentration of PaoABC which resulted in a complete conversion to FDCA without any HMFA being formed. This goes to show that not only pH [20], but also substrate and enzyme concentrations can affect conversion rates and the oxidation route of certain enzymes.
Based on the reactions presented in Table 3, approximate kinetic parameters were also calculated. The time course oxidation of HMF with AO (Fig. 2a) showed that most of the substrate was converted in the first 24 h, therefore the kinetic parameters were determined based on the amount of product formed in that time period. The results in Table 4 show that the highest V 0 Table 3     were calculated for AO with HMF, ranging from 0.08 to 0.2 mM/h, and the lowest for LAC with HMFA, with the highest being 0.0002 mM/h and the rest at 0 mM/h. AO also had the highest affinity for its substrate, with a K m of 14 mM. However, the highest V max was determined for GAO with DFF, despite the K m of 770 mM showing a very low affinity for this substrate. Also, the productivity of HRP with FFA (product yield: 0.6-0.9%) was lower than that of GAO with DFF (product yield: 1.3%), yet the K m value of HRP with FFA was calculated to be 260 mM, almost 3 times lower than that of GAO with DFF. This discrepancy could be due to: (1) the fact that we only tested 3 different substrate concentrations (enzyme saturation was not reached); (2) that the conversion yields of DFF and FFA were very low and thus less accurate, and (3) that the Lineweaver-Burk plot is prone to error [29]. Thus, the values of V max and K m are most likely an overestimation.
To the best of our knowledge, kinetic parameters for this reaction have been determined for only a few enzymes [11,13,14,20]. The highest affinities presented in literature were determined for UPO and DFF with K m of 0.8 mM [20] and HMFO and HMF with K m of 1.4 mM [13]. A K m of 1.6 mM was determined for one AAO with HMF [11], which is almost 10 times lower than that of our AO. However, in another study, two other AAOs were tested against HMF and their K m values were 36.3 and 7.2 mM, respectively [20]. The same group also tested GAO with HMF and UPO with HMF and DFF. The K m value for GAO with HMF was 142 mM, indicating a relatively low affinity for this substrate. The group did not test GAO with DFF and, thus a comparison to our results is not possible. Although our calculated K m value for GAO with DFF is most likely an overestimation, based on the data on GAO with HMF [20] and the fact that our results showed the conversion of DFF with GAO is lower than that of HMF, a high value of K m for DFF was to be expected. Also to be taken into account is the type of measurement used to obtain the raw data for the determination of kinetic parameters. In our study, we directly measured the amount of product formed, whereas some other groups used a more indirect method and measured the amount of H 2 O 2 produced during the enzymatic conversion of a substrate in a peroxidase-coupled assay [11,13,14]. Since different measurements yield different types of data, a comparison of our results to those presented in literature might not be very accurate.

One-pot multi-enzyme reactions
With the aim of achieving a full conversion of HMF to FDCA or at least increasing the productivity of an enzyme, we devised simple one-pot multi-enzyme reactions (Fig. 4a-c; values of single-enzyme reactions for AO and GAO were taken from Table 2). AO and GAO were chosen as the main enzymes because of their relatively high activity against HMF compared to other selected enzymes. Since these two enzymes produce H 2 O 2 which in turn inhibits their activity [22,24], CAT and HRP were added to the reaction to remove this inhibitor and in that way increase the activity of AO and GAO. Results in Fig. 4a show that adding CAT to AO increased the production of DFF from 25.5 to 97.5%, but decreased the production of FFA from 3.1 to 0%. Adding HRP to AO and CAT had a similar effect, but to a smaller degree and produced 76.6% of DFF. Figure 4b shows that adding CAT to GAO slightly increased the production of DFF, from 5.1 to 7.4% and adding HRP further increased it to 18.1%. Combining all four enzymes (Fig. 4c) resulted in production of 36.7% of DFF, but without any other products being formed.
The oxidation of HMF to FDCA is rarely performed by only one enzyme as three consecutive oxidation steps are needed and thus, an enzyme would also have to accept all of the intermediates (DFF, HMFA and FFA) as a substrate. One of the exceptions was an HMF oxidase enzyme (HMFO) from Methylovorous sp. which was capable of oxidizing not only HMF, but also DFF, HMFA and FFA and thus producing high yields of FDCA [13,14]. Since no other enzyme has been shown to have such activity and wide range of substrates, other groups had to rely on finding the most favourable combinations of different enzymes to achieve a full conversion. Most often, an oxidase (alcohol or galactose) was combined with CAT, HRP or an unspecific peroxygenase [11,20,[22][23][24]. The latter were added not only because of their ability to utilize/remove the inhibitory H 2 O 2 produced by the oxidases, but also because of their activity against other substrates like FFA. Similar to previous research, our results show that unlike HMFO [13,14], none of the selected enzymes were able to catalyse the whole conversion of HMF to FDCA but that the addition of an enzyme that utilizes H 2 O 2 can increase the activity of another enzyme (an oxidase). The most promising enzyme combination used in our work was shown to be AO and CAT. The fact that adding HRP to the reaction lowered the production of DFF indicates that some inhibitory enzyme interactions might be taking place. This is also supported by results obtained with all four enzymes, where compared to only AO and CAT and AO, CAT and HRP, the production of DFF was lower. It has previously been shown that adding CAT and HRP to GAO increases the production of DFF from 2 to 28% [22]. There they also tested GAO with only CAT and GAO with only HRP and the results showed that while the addition of CAT and HRP increased DFF production to 23% and 46%, respectively, the combination of all three enzymes did not give the highest yield. This again points to an inhibitory enzyme interaction. On the other hand, our results show a more pronounced rise in DFF production when HRP was added to the reaction. However, we did not test GAO with only HRP and thus could not determine its effect on GAO. We would expect that this combination would yield the best results as HRP has been shown to activate GAO [11,23,30,31].

Conclusions
This study looked at the mechanism of enzymatic conversion of HMF and its oxidation products to FDCA. The activity was determined based on a single-enzyme active site and at atmospheric oxygen. Different enzymes had different substrate specificities as well as employed different HMF oxidation routes, but most had only low activity against the selected substrates. The highest yields were obtained for both oxidases. Substrate concentration affected only alcohol oxidase activity. Combining different enzymes in multi-enzyme reactions increased the conversion of HMF to DFF, and the best results were obtained with AO and CAT.


text:
Isolation and Pathogenic Characterization of Pigeon Paramyxovirus Type 1 via Different Inoculation Routes in Pigeons

Pigeon paramyxovirus type I (PPMV-1) causes regular outbreaks in pigeons and even poses a pandemic threat among chickens and other birds. The birds infected with PPMV-1 mainly show a pathological damage in the respiratory system, digestive system, and nervous system. However, there were few reports on the efficiency of the virus entering the host via routes of different systems. In the present study, a PPMV-1 strain was obtained from a dead wild pigeon in 2016 in Beijing, China. The mean death time (MDT) and the intracerebral pathogenicity (ICPI) of our isolate showed medium virulence. Phylogenetic analysis based on F gene sequence showed that the isolate belonged to subgenotype VIb, class II, which dominated in China in recent years. Then, we evaluated the infection efficiency of different routes. Pigeons were randomly divided into five groups of six as follows: intracephalic (IC), intranasal (IN), and intraoral (IO) infection routes, cohabitation infection (CO), and negative control (N negative). All pigeons were inoculated with 100 μl·106 EID50 PPMV-1 virus. After infection, pathological lesions, virus shedding, body weight change, survival rate, and tissue tropism were tested to compare the efficiency of the different infected routes. The mortality of groups IC, IN, IO, and CO were 100, 66.7, 50, and 33.3%, respectively. Weight loss in group IC was higher than the other groups, followed by groups IN and IO. The lesions observed in PPMV-1-infected pigeons were severe, especially in the lung and intestine in group IC. Viral shedding was observed from 2 dpi in groups IC and IN, but the shedding rate was higher in group IN than group IC. The longest period was in group CO. Tissue tropism experiment showed that our isolate has a wide range of tissue distribution, and the virus titer in the heart and intestine of group IC and in the brain of group IN was higher. Our data may help us to evaluate the risk of transmission of PPMV-1.

INTRODUCTION
Newcastle disease, caused by Newcastle disease virus (NDV), is classified as a notifiable disease by the World Organization for Animal Health (OIE), because of its high morbidity and mortality in avian species. At least 250 avian species, including chickens, ducks, and pigeons were reported to be susceptible to NDV (1)(2)(3). NDV is a single-stranded, negative-sense, and non-segment RNA virus, which belongs to the genus Orthoavulavirus, subfamily Avulavirinae, family Paramyxoviridae of the Mononegavirales order (4)(5)(6)(7). NDV isolates have been classified into classes I and II. NDV of class I are mostly of low virulence and contain only one genotype (8), and class II virus consists of 18 genotypes (9-13) based on the phylogenetic analysis of nucleotide sequence of the F gene and genomic size (14,15). Pigeon paramyxovirus type 1 (PPMV-1) found in birds of Columbidae family, mainly turtle dove (Streptopelia turtur) and Eurasian collared dove (Stretopelia decaocto) is a variant of Newcastle disease virus, and it is almost genotype VI, class II. PPMV-1 first reported in England during the late 1970s, affected racing pigeons with outbreaks in domestic chickens (16), which caused the third worldwide pandemics of NDV (17). In recent years, a high mortality (from 40 to 80%, even to 100% in some cases) caused by PPMV-1 has been observed in pigeons (18)(19)(20)(21)(22). Clinical signs of the infected birds involve nervous, respiratory, and digestive system symptoms (23), consisting of moderate to severe depression with neck twisting, ataxia, crouch, paralysis, eyelid edema, diarrhea, and green loose feces. NDV has been documented to remain infectious in feces and carcasses for at least a couple of weeks, several months in feathers, and up to 90 days in soil or water (24). Quite a few times, outbreaks in poultry have been ascribed to PPMV-1 (25)(26)(27)(28). The virulence can be enhanced after serial passages in chickens (17,21,29,30), which makes these pigeon-originated viruses a tremendous and continuous threat to the poultry industry (31)(32)(33). For these reasons, the potential transmission of the virus and the effective route where the virus will enter the host have been considered points of concern. In the present study, the pathogenicity of a PPMV-1, obtained from a dead wild pigeon in 2016 in Beijing, China, was investigated via different inoculation routes. Findings from our study showed intracephalic, intranasal, and intraoral infection routes were effective, but intracephalic was the most.

Viral Isolation, Amplification, and Full-Length Genome Sequencing
A moribund pigeon with neck twisting, diarrhea, and leg paresis or paralysis was found in Beijing, 5 September 2016. We initially diagnosed that the pigeon was infected with PPMV or avian influenza virus based on clinical symptoms, then, the avian influenza virus was excluded and PPMV-1 infection was confirmed by RT-PCR. The identification of Newcastle disease and the separation of strains are as follows: viral RNA was extracted from the tissues of the pigeon (i.e., heart, liver, spleen, lung, kidney, stomach, brain, trachea, intestine, and pancreas) using Trizol reagent (ambition by Life Technologies, Beijing, China) according to the manufacturer's instructions. Reverse transcription was performed as described (34). The detection gene (a part of fusion gene, 486 bp) was amplified from the cDNA by PCR utilizing Taq DNA Polymerase (CWBIO 2× Taq MasterMix, Cat. CW0682M), and the primers were designed according to the conserved sequence (34) (Primer sequence: F:CAG CTGCGGCCCTAATACA; R:TGGATGCCCAAGAGTTGAG). The program was as follows: 95 • C for 5 min; 35 cycles of denaturation at 95 • C for 30 s, annealing at 55 • C for 30 s, and extension at 72 • C for 25 s; and a final extension at 72 • C for 10 min. The PCR products were visualized by 1% agarose gel electrophoresis. Other pathogens (circovirus, avian influenza, and pathogenic bacteria) were negative. The viruses from different organs were plaque-purified three times on primary chicken embryo fibroblasts and inoculated into the allantoic cavity of 9-day-old specific-pathogen-free (SPF) eggs. The virus was isolated and RNA was extracted from the allantoic cavity, detected using PCR as abovementioned. The strain was designated as PPMV-1/pigeon/Beijing/China/01/2016 and abbreviated as PPMV-1/BJ-01/CH. Complete genome of our strain was amplified (primers are shown in Supplementary Table 1), sequenced, and submitted to GenBank (GenBank ID: MH807446).

Virulence-Test
There are three virulence evaluation indexes on Newcastle disease virus, and virulence is usually determined by no <2 indexes. The mean death time (MDT) was determined in 9-day-old SPF chicken embryo eggs, and the intracerebral pathogenicity index (ICPI) was determined in 1-day-old chick as previously described (35). The least-fatal dose and egg-infectious dose of the virus were tested in 9-day-old SPF chicken embryo eggs by multiple proportion dilution (from 10 −3 to 10 −9 ) with five repeats per dilution and calculated using the Reed and Muench method.

Phylogenetic Analysis Based on Complete F Gene
To determine the genetic relationships of the PPMV-1/BJ-01/CH isolate to others, a phylogenetic tree was constructed using the MEGA 7.0 software with the maximum likelihood method via the Kimura two-parameter model based on the complete F gene from 21 different subgenotypes of the reference PPMV-1 isolates and 1 APMV-1 as outgroup (8,12,29). All complete sequences of the F gene were downloaded from GenBank and aligned utilizing ClustalW. Nucleotide sequences obtained in the present study were aligned with reference sequences available in the GenBank database to determine the subgenotypes.

Animal Experimental Infection
To further determine the pathogenicity of the virus, a total of 31-month-old pigeons, with approximately equal body weight of −5-5%, were used in our study. These pigeons were bought from a hatchery in Miyun District, Beijing, and certified by hemagglutination inhibition (HI) experiment to have no antibodies of NDV and AIV. These pigeons were randomly divided into five groups, and a marked group was placed in a separate cage in an animal room under biosafety conditions. Adequate food and drinking water were provided. Pigeons were inoculated with the virus by intranasal, intraoral, intracephalic, and cohabitation infections (signed as groups IN, IO, IC, and CO, placing groups CO and IC together), which represented through respiratory system, digestive system, nervous system, and natural infection, respectively. Additionally, the negative control group received phosphate-buffered saline (PBS) solution at pH 7.2. All infected groups were inoculated with a dose of 10 6 median embryo lethal dose (ELD50)/100 µl each, calculated using the Reed and Muench method. Subsequently, all pigeons were observed daily for clinical signs, and clinical symptoms, mortality, and morbidity were recorded. Oropharyngeal and cloacal swabs were taken every day, and body weight was determined every other day until day post-infection (dpi) 14 (0 day post-infection in group CO means the day the pigeons were infected), when there is only one group left to shed virus.

Virus Shedding
All of the oropharyngeal and cloacal swabs were collected, placed in tubes with phosphate-buffered saline solution and 2% fetal bovine serum and stored at −80 • C until RNA extraction.
RNA was extracted, and reverse transcription and RT-PCR were performed as above to test the virus shedding. In addition, cDNA of the isolate in this study and reagent-grade water were used as positive and negative control, respectively.

Gross Lesion
Lesion formation and lesion size are often used to quantify virulence (36). In our study, organs of all dead pigeons were collected on dpi 5 and 14. The collected tissues (including the heart, brain, lung, intestine, and liver) were pathologically lesion examined. Individuals that died ahead of time point were examined and recorded in advance. Positive rate was tallied by percentage.

Histopathology
Tissues from each group were collected and fixed with 10% neutral formalin. The sections were stained with hematoxylin and eosin (HE), and all HE stained sections were examined for the presence of microscopic lesions.

Tissue Distribution
To understand the distribution and virus load of PPMV-1 in organs from different groups at dpi 5 when clinical symptoms emerged, we randomly choose three from each group, and collected samples from the heart, brain, lung, intestine, and liver. The virus was isolated from tissues of the same weight, and TCID 50 in DF-1 cells (chicken fibroblast cell line) was used  to estimate viral loads of the five groups; 3 × 10 4 DF-1 cells were seeded in 96-well plate with five repetitions 1 day before infection. Twenty-four hours later, the cells were infected with different dilutions of the virus for 1 h at 37 • C with shaking every 12 h and confirmed by the hemagglutination assay. TCID 50 was calculated using the Reed-Muench method. Data were analyzed using Prism (v.5.01) software. Statistical significance was set at a P-value of <0.05.

Ethics Statement
These animal studies were performed in strict accordance with the Guidelines for the Care and Use of Animals in Research, which are issued by the Institute of Zoology, Chinese Academy of Sciences (Approval Number IOZ12017).

Phylogenetic Analysis and Genetic
Characteristics of PPMV-1/BJ-01/CH Isolate Phylogenetic analysis showed our isolate belonged to subgenotype VIb (Figure 1). Our isolate possessed an 11 2R-R-Q-K-RF 117 segment at the cleavage site of F gene. This is characteristic of virulent NDVs. This motif is commonly found in NDV strains that are highly or moderately virulent in chickens, especially in genotypes VII and IX viruses and some pigeon paramyxovirus strains. Amino acid residues important for receptor recognition of HN is 174 (37), which is the same as most PPMV-1 strains.

Virulence of the PPMV-1/BJ-01/CH Isolate and Sequence Analysis
The MDT and ICPI were 72 h and 0.76, respectively. The results indicated that the PPMV-1 strain was medium according to the criteria (38). Moreover, the minimum lethal dose and 50% egg infection dose of the virus in chick embryo eggs were 10 −3 and 10 −5.37 , respectively.

Pathogenicity of the PPMV-1/BJ-01/CH in Different Routes of the Infection
Clinical signs were observed in all pigeons from 5 dpi. The pigeons lost their weight sharply once clinical symptoms emerged, then soon died. The body weight change (Figure 2A) and survival curve (Figure 2B) showed the damage of the infected pigeons from group IC were much severer than other groups. A steady increase of body weight was found in negative control group. Post-mortem examination showed the congestion or hemorrhages on the meninx and in the brain. The tracheal The lesions were shown in Supplementary Figure 1. Frontiers in Veterinary Science | www.frontiersin.org mucosa was congested or hemorrhagic. Extensive hemorrhages were observed in the mucosa of the small intestine. The spleens were atrophic, friable, and hemorrhagic. Congestion and hemorrhage were found in the pancreas. Our data showed the most serious tissue lesion was discovered in group IC. Groups IO and CO caused less lesions on dpi 5 ( Table 1).
Histopathological examination showed that the most severe lesions are found in group IC, including the desquamation of small intestine villi, lymphocyte infiltration, and severe epithelial cell loss in the intestine. The alveoli have been destroyed, and a small amount of lymphocyte infiltration and congestion were observed in the lungs. Hepatomegaly and the disappearance of hepatic cord structures were observed in the liver. A small amount of lymphocyte infiltration was also observed in the cerebrum. Milder lesions were observed in other groups (Figure 3).

Virus Shedding
The results and mortality of each group are compiled in Table 2. Pigeons shed virus from dpi 2 to dpi 14. The virus was detected from oropharyngeal swabs earlier than cloacal swabs in all groups. The earliest start time observed was in groups IC and IN in dpi 2, but the shedding rate was higher in group IN than group IC. The highest shedding rates were observed in group IN from the throat at 2 dpi and cloaca at 7 dpi. Viral shedding ceased at 14 dpi in all groups ( Table 2).

Tissue Tropism of the PPMV-1/BJ-01/CH in Different Routes of Infection
In this study, we found that the virus can effectively replicate in different tissues, no matter which route of infection, and virus titer in the brain, lung, and intestine were higher than in other tissues. The highest titer of the virus was found in the intestine in group IC and brain in group in (Figure 4).

DISCUSSION
In the present study, a PPMV-1 was isolated from a dead pigeon. In order to better understand the genetic characteristics of the strain, the full-length genome was amplified and sequenced. Phylogenetic analysis showed that the isolate belonged to VIb, class II, which was consistent with most of the genotypes found in China (20,21,(39)(40)(41). Amino acid sequence of F gene 112-117 is a typical virulent motif. But the MDT and ICPI showed medium virulence of our isolate. This contradiction often occurs on PPMV-1. Amino acid residues important for receptor recognition of HN are the same as most PPMV-1. All the characteristics of our strain consisted with current epidemic isolates in China (20,42,43).
The pathogenicity of the virus may relate to different inoculation routes. For example, the ducks infected intramuscularly with a virulent NDV strain showed the most severe clinical signs, while ducks infected intranasally and intraocularly sometimes also exhibited clinical signs but seldom died (44). To further understand the target system of the virus obtained in this study, several experiments were conducted to evaluate the pathogenicity via different system routes. Cohabitation infection was intended to mimic a naturalacquired infection. Intranasal infection was often used by labs as a substitute for aerosol infection of the virus. Oral infection mimic was used to simulate contaminated food and water. Inoculate intracranially was regularly used to evaluate NDV virulence, so it may be the most effective route of infection.
Viral shedding evaluates the effective replication and transmission of the virus in different groups. The results of our study showed infected pigeons of different groups except group CO which shed virus through the larynx from 2 dpi and through cloaca from 4 dpi. Pigeons shed virus through cloaca from 2 dpi after isolating NDV in previous studies (45,46). The earliest start time observed was in groups IC and IN, and the shedding rate was higher in group IN than group IC; the highest shedding rates were observed in group IN from the throat at 2 dpi and cloaca at 7 dpi. The result indicated that the respiratory tract is the fastest way for PPMV-1 to spread. The shedding rate in group IC remained high, which may confirm that virus can efficiently enter into the hosts through the nervous system. All infection groups shed virus, and the longest persisted until 14 days. Continuous virus shedding may contribute to circulatory infection. Additionally, we used body weight and survival rate to test the pathogenicity. Death occurred in dpi 6, but previous study of a PPMV-1 was inoculated in pigeons through the intranasal route did not result in mortality up to 31 dpi. (47), it may be caused by the differences in virulence of the isolates used (48). The rate of body weight decline of group IC was higher than in the other groups, followed by groups IN and IO. According to the result of TCID 50 in DF-1 cells, PPMV-1 was able to cause systemic infection in a relatively short time, and virus titer in groups IC and IN was higher than in other groups, which was broadly consistent with gross lesions at dpi 5, resulting in an effective infection process. Groups IO and CO caused less lesions which might be due to lower viral load in the organs at these days. By comparison of the viral load in different tissues, the virus had a wide range of tissue distribution, especially in the lung, brain, and intestine. It can be inferred that the virus could replicate well in these three organs during early infection. At 14 dpi, the gross lesions were more severe in these groups, but the rate of virus shedding was reduced which may be due to higher mortality rates. The results indicated that the virus can effectively infect pigeon via different routes, and the most pathogenic was infection through the nervous system and respiratory system, but infection through the nervous showed stronger pathogenicity. Previous studies described neurological lesions in NDV-infected birds, and virulent virus was capable of replication in the brain, but not the avirulent virus (49)(50)(51), thus successful replication in the nervous system determines its pathogenicity.
Overall, pigeons play an important role in the epidemiology of PPMV-1. The routes of inoculation greatly influenced pathogenicity of the pathogenic strain isolated in the present study. Considering the growing number of PPMV-1 cases in recent years, it is necessary to develop effective vaccines or other prevention and control methods.

DATA AVAILABILITY STATEMENT
Complete genome of our strain was amplified (primers were showed in Supplementary Table 1), sequenced and submitted to GenBank (GenBank ID: MH807446).

ETHICS STATEMENT
The animal study was reviewed and approved by Institute of Zoology, Chinese Academy of Sciences (Approval Number IOZ12017).


extfieldsofstudy: Biology Medicine

s2fieldsofstudy:
Agricultural And Food Sciences
Biology
Medicine

text:
Antimicrobial Resistance Profiles and Characterization of Escherichia coli Strains from Cases of Neonatal Diarrhea in Spanish Pig Farms

Escherichia coli is considered one of the most common agents associated with neonatal diarrhea in piglets. The aim of this work was to characterize the pathogenic and antimicrobial resistance (AMR) profiles of 122 E. coli strains isolated from pigs suffering diarrhea (n = 94) and pigs without diarrhea (n = 28) of 24 farms in Spain. Virulence factors, toxins and AMR (ESBL and colistin) genes and AMR phenotypes of E. coli isolates were analyzed. Low prevalence of pathogenic E. coli strains (26%) was found in both groups. However, ETEC and VTEC strains were more frequently isolated from diarrheic piglets. Irrespectively of diarrhea occurrence, 97.5% of the strains showed a multidrug-resistance (MDR) profile to aminopenicillins, sulfonamides and tetracyclines. It was found that 22% of E. coli was CTX-M+, with CTX-M-14 being the principal allelic variant. Remarkably, 81.5% of CTX-M+ strains were isolated from diarrheic animals and presented an extended MDR profile to aminopenicillins, quinolones and aminoglycosides. Finally, low frequencies of colistin resistance genes mcr-1 (4/122) and mcr-4 (1/122) were found. MDR E. coli strains are circulating in pig farms of Spain, representing a serious threat to animal and public health. More appropriate diagnostic approaches (genetic and AMR phenotypic analysis) should be implemented in animal health to optimize antibiotic treatments.

Introduction
Escherichia coli has historically been considered as one of the most common agents associated with diarrhea in suckling and post-weaned piglets [1]. There is high diversity and variants of E. coli strains integrating the normal gut microbiota, with most of them being considered not pathogenic [2]. The characterization of pathogenic E. coli strains is usually based on the presence of virulence factors [3]. In piglets, E. coli pathogenic strains can be classified into different pathotypes: enterotoxigenic (ETEC) strains releasing heat-labile (LT) and heat-stable Sta and Stb exotoxins, intimin (eae)-producing enteropathogenic (EPEC) strains and verotoxigenic (VTEC) strains producing VT1/VT2 verotoxins [4].
The routine use of antimicrobials in livestock, especially in the pig industry, for either 'prophylaxis' or 'metaphylaxis' represents a serious hazard for the selection of multidrug-resistant (MDR) Vet. Sci. 2020, 7, 48 2 of 12 Enterobacteriaceae strains [5,6]. In fact, the effectiveness of treatments against E. coli is threatened by the dramatic increase of extended-spectrum beta-lactamases (ESBL)-producing isolates worldwide [7][8][9]. A high number of ESBL have been described and designated according to the strains of bacteria or plasmid that produce them by using letters that usually evoke the name of the target antibiotic, but they are also named after substrates, their biochemical properties, location of their discovery and the name of the patient or the discovering investigators [10]. Genes encoding ESBL enzymes are usually located in plasmids and can be easily horizontally transferred among different bacterial populations. Most ESBL derive from the first enzymes discovered, TEM (named after the patient (Temoneira) providing the first sample) and SHV (sulfhydryl reagent variable), which confer resistance to beta-lactams such as ampicillin or amoxicillin. The most recent ESBLs are derived from the CTX-M (active on cefotaxime-Munich) enzymes, which confer resistance to third-and fourth-generation cephalosporins [10,11]. Cephalosporins constitute one of the largest families of antimicrobials widely used in both human and veterinary medicine [12]. On the other hand, in veterinary medicine, colistin sulfate is indicated for the treatment of Enterobacteriaceae infections and has been widely used as a preventive mass-medication of colibacillosis in piglets [13]. Although colistin is one of the antimicrobial agents with the lowest rates of resistance, the emergence of the mcr-1 (mobile colistin resistant) plasmid in the Enterobacteriaceae population has resulted in the appearance of strains with acquired resistance to this antibiotic [14,15]. In order to stop the increase of colistin-resistant strains, the European Medicines Agency established in 2016 that all EU members should restrict and reduce the use of colistin in animals for treating infections with a target level of 5 mg/PCU, where PCU refers to the 'population correction unit' and takes into account the animal population as well as the estimated weight of each particular animal at the time of treatment with antimicrobials. [16].
During recent decades, the importance of E. coli in porcine neonatal diarrhea seems to have decreased [17], possibly due to the successful vaccination plans implemented in the pig farms [18]. Therefore, antimicrobial treatment might not be required in most cases of diarrhea, especially in those where the role of the E. coli is not well defined. However, the systematic mass-medication of piglets suffering from diarrheic processes is still common in many cases.
The aim of the present work was to characterize the virulence factors and to determine the genotypic and phenotypic antimicrobial resistance patterns of E. coli isolates from cases of neonatal diarrhea in conventional pig farms of Catalonia (Spain). Moreover, in order to assess the clinical relevance of all isolates, samples from non-diarrheic pen-mates were also analyzed and compared to the clinical cases.

Microbiological Testing
A total of 122 E. coli isolates obtained from 24 conventional Spanish farms suffering neonatal diarrhea outbreaks during 2017 and 2018 were recovered from a previous study [19]. All the studied farms were located in Catalonia (NE of Spain), the Spanish region with the greatest number of production farms and one of the highest pig-density (242 animals/km 2 ) regions in Europe. The sampling procedure included 5-10 samples from diarrheic animals and 3-5 samples from apparently healthy pen-mates for each farm. One gram of fecal sample was obtained directly from the animals using rectal swabs and submitted for diagnostic testing to the Laboratori Veterinari de Diagnosi de Malalties Infeccioses of the Universitat Autònoma de Barcelona (Spain). Besides E. coli, all samples were tested for a panel of enteric infectious agents which comprised C. perfringens types A and C toxins (Cpα, Cpβ, Cpβ2); C. difficile toxins (TcdA, TcdB); rotavirus A (RVA), B (RVB) and C (RVC); porcine epidemic diarrhea virus (PEDV) and transmissible gastroenteritis virus (TGEV), and these general results were described elsewhere [19].
For E. coli isolation, stool samples were cultured in standard Blood Agar (BD GmbH, Germany) and MacConkey Agar (Oxoid, UK) and aerobically incubated at 37 • C for 24 h. A total of 122 pure isolates were collected: 94 from diarrheic samples of 1-week-old diseased piglets and 28 from feces of apparently healthy pen-mates. All strains were confirmed as E. coli using conventional biochemical tests. Only E. coli strains obtained as pure cultures from a direct seeding of feces were used for the characterization, since the lack of microbial diversity was considered abnormal.

Characterization of E. coli Virulence Factors and Toxins
The presence of E. coli toxins (Sta, Stb, LT, EAST1, VT1 and VT2), fimbrial adhesins (F4, F5, F6, F18 and F41) and non-fimbrial adhesin eae genes was analyzed by qualitative PCR, as described by Toledo et al. (2012) with slight modifications [20]. Briefly, PCR assays were carried out using the Biotaq DNA polymerase kit (Biotaq, Ecogen) using previously described primers [19]. The final 25 µL mixture consisted of: 1× PCR Buffer, 0.8 mM of the dNTP mix, 3 mM of MgCl 2 , 1 mM of each primer and 1 U of Taq polymerase. As template, 2.5 µL of the DNA sample was used. The PCR program consisted of 5 min at 94 • C, followed by 30 cycles of 1 min at 94 • C, 1 min of annealing at 57 • C and 1 min of extension at 72 • C, and a final extension step of 7 min at 72 • C. Positive and negative (ultrapure water) controls were included in each run and particular care was taken to prevent carry-over contamination.

Antimicrobial Resistance Phenotyping Analyses
Antimicrobial susceptibility of E. coli isolates was tested using the disk diffusion method [21]. Briefly, one to four colonies were suspended in 5 mL of distilled sterile water to achieve a turbidity of 0.5 in the McFarland scale. The dilution was then seeded onto Mueller-Hinton (Oxoid, UK) plates. Each isolate was tested for the following antimicrobial groups, using commercial disks: aminopenicillins (ampicillin, amoxicillin and amoxicillin/clavulanic acid); cephalosporins (ceftiofur, cephalexin, cefquinome, ceftriaxone); quinolones (ciprofloxacin, enrofloxacin, flumequine); aminoglycosides (gentamicin, neomycin, streptomycin, apramycin); tetracyclines (tetracycline, doxycycline); sulfonamides (sulfonamide, sulfamethoxazole/trimethoprim); florfenicol; colistin; lincospectin. Concentration of the commercial disks and breakpoints for each antimicrobial are shown in Table 1. Additionally, minimum inhibitory concentration (MIC) test was performed in order to evaluate antimicrobial susceptibility to colistin, using the broth microdilution method in 96-well plates [22]. E. coli ATCC 25,922 was used as the quality control strain. Briefly, the tested colistin concentrations ranged from 0.25 to 8 µg/mL. The tested strains were considered resistant when their MIC value was higher than the wild type cut-off value (MIC > 2µg/mL) [23].

Genetic Characterization of Antimicrobial Resistance Genes
Molecular diagnosis of AMR genes was performed for the ESBL genes bla SHV , bla CTX-M , bla CMY1 , bla CMY2 and bla TEM ; carbapenemase bla OXA and mcr-1 to mcr-5 colistin-resistance genes [25]. PCR conditions were homogenized for all the reactions as follows: 5 min at 94 • C, followed by 25 cycles of 1 min at 94 • C, 1 min of annealing at 55 • C and 1 min of extension at 72 • C, and a final extension step of 7 min at 72 • C. Amplified PCR products were Sanger sequenced for verification at the Genomic and Bioinformatics Service of the Universitat Autònoma de Barcelona (Barcelona, Spain). Sequences and chromatograms were manually explored to trim bad-quality bases using the BioEdit software. The obtained sequences were aligned using Clustal Omega and blasted against those described at the National Center for Biotechnology Information (NCBI) database.

Statistical Analysis
Chi-squared or Fisher exact tests were used for comparison between proportions when appropriate. Results were considered statistically significant when p < 0.05.

Microbiological Identification and Characterization of E. coli Isolates
Most of the E. coli isolates presented a very low prevalence of virulence factors and toxins (<5%). with the exception of EAST1, which was present in more than 70% of the strains (Table 2). Moreover, 74% of the strains were negative for all toxins and intimin (eae) genes. As such, it was not possible to classify them into any pathotype.  Table 3 shows the number of isolates presenting virulence factors and toxins grouped by pathotype and their total frequency in diarrheic and not diarrheic cases. The most common pathotypes within the classified strains were EPEC (n = 17) and ETEC (n = 15), with the frequency of EPEC being two times higher (21% vs. 10.6%) in non-diarrheic piglets. As expected, the frequency of ETEC was higher in animals presenting diarrhea than in those without clinical signs (15% vs. 3.5%). The presence of VTEC (4.3%) was only detected in cases of diarrhea. Similarly, E. coli strains positive for F4, LT, Sta, VT1 and VT2 genes were only isolated from clinical cases. Table 3. Description of the profiles of virulence factor production and classification of the toxigenic strains into pathotypes of E. coli strains.

Antimicrobial Resistance Phenotyping and Genotyping
According to the results obtained by the disk diffusion method, the overall prevalence of AMR in E. coli isolates was very high for most of the nine tested antimicrobial classes. The highest rates of resistance were found for sulfonamides (87.3%), tetracyclines (87.3%) and aminopenicillins (80.1%), followed by quinolones (49%), cephalosporins (32.9%) and aminoglycosides (27.9%). By contrast, low resistance levels were observed for colistin (6.5%) and apramycin (5.9%). Regarding the MIC results, 6 out of 122 (5%) were resistant to colistin (MIC > 2 µg/mL): 5 of them had an MIC value of 4 µg/mL, and one isolate had an MIC of 8 µg/mL. Comparing the colistin results obtained by disk diffusion and the MIC, only one strain was considered resistant to colistin by both methods.
The comparison study between E. coli isolates from diarrheic and non-diarrheic samples showed statistical differences in the proportion of AMR E. coli isolates ( Figure 1). Remarkably, piglets presenting diarrhea had higher prevalence of E. coli strains resistant to quinolones and aminoglycosides (gentamicin and neomycin) compared to the non-diarrheic isolates (Fisher exact test p < 0.05). 星Total n (%) 1星3 (3.2%) 0星0 1星3 (2.5%)

Antimicrobial Resistance Phenotyping and Genotyping
According to the results obtained by the disk diffusion method, the overall prevalence of AMR in E. coli isolates was very high for most of the nine tested antimicrobial classes. The highest rates of resistance were found for sulfonamides (87.3%), tetracyclines (87.3%) and aminopenicillins (80.1%), followed by quinolones (49%), cephalosporins (32.9%) and aminoglycosides (27.9%). By contrast, low resistance levels were observed for colistin (6.5%) and apramycin (5.9%). Regarding the MIC results, 6 out of 122 (5%) were resistant to colistin (MIC > 2 µg/mL): 5 of them had an MIC value of 4 µg/mL, and one isolate had an MIC of 8 µg/mL. Comparing the colistin results obtained by disk diffusion and the MIC, only one strain was considered resistant to colistin by both methods.
The comparison study between E. coli isolates from diarrheic and non-diarrheic samples showed statistical differences in the proportion of AMR E. coli isolates (Figure 1). Remarkably, piglets presenting diarrhea had higher prevalence of E. coli strains resistant to quinolones and aminoglycosides (gentamicin and neomycin) compared to the non-diarrheic isolates (Fisher exact test p < 0.05).  On the other hand, E. coli harboring CTX-M+ genes presented a wider resistance pattern than those CTX-M-negative strains (CTX-M-). Specifically, 81.5% of the CTX-M+ strains showed resistance to more than 10 antimicrobials, while in the case of CTX-M-the percentage was 37.9%. Statistical differences between CTX-M+ and CTX-M-strains were found for amoxicillin/clavulanic acid, cephalosporins, quinolones and gentamicin drugs (Figure 2). By contrast, colistin resistance genes were detected in a low proportion of cases: four strains harbored the mcr-1 gene (3%), and one strain was positive to mcr-4. Only one of the mcr-1+ strains, isolated from a diarrheic animal, exhibited phenotypic resistance to colistin in the MIC test (8 µg/mL).

% Resistant strains
On the other hand, E. coli harboring CTX-M+ genes presented a wider resistance pattern than those CTX-M-negative strains (CTX-M-). Specifically, 81.5% of the CTX-M+ strains showed resistance to more than 10 antimicrobials, while in the case of CTX-M-the percentage was 37.9%. Statistical differences between CTX-M+ and CTX-M-strains were found for amoxicillin/clavulanic acid, cephalosporins, quinolones and gentamicin drugs (Figure 2). By contrast, colistin resistance genes were detected in a low proportion of cases: four strains harbored the mcr-1 gene (3%), and one strain was positive to mcr-4. Only one of the mcr-1+ strains, isolated from a diarrheic animal, exhibited phenotypic resistance to colistin in the MIC test (8 µg/mL).  Table 1.

Discussion
The aim of this work was to characterize the E. coli strains involved in porcine neonatal diarrhea outbreaks by genotyping and phenotyping isolates grown as pure cultures from both clinical cases of diarrheic and non-diarrheic healthy pen-mates.
Surprisingly, the general prevalence of strains encoding genes for virulence factors, as well as the presence of the typical E. coli pathotypes, was low in our pig population studied. Merely 26% (32/122) of the strains could be classified into a specific pathotype. EPEC was the most frequent one, representing 14% of the total isolates, followed by ETEC (12%) and VTEC (3%). The presence of EPEC is usually reported as a causative agent of diarrhea, principally in post-weaned but also in suckling pigs, but in this study, the presence of this pathotype was found more prevalent in non-diarrheic piglets (21.4%) than in the diarrheic ones (11.8%). Similar results have been reported in previous works in agreement with the results presented in the present study. Mainly, EPEC prevalence ranging from 3% to 26% in diarrheic animals and 9% to 17% in non-diarrheic animals has been described [26][27][28]. Finally, EAST1 gene was highly detected in both diarrheic (70%) and non-diarrheic groups (86%) but again with no relationship with the clinical onset. The EAST1 toxin is produced by several E. coli pathotypes [29], and although initially it was believed to be involved in the pathogenesis of diarrhea, most authors have concluded that this toxin may not be a virulence determinant, because it has been commonly isolated from healthy animals as well [30,31].
On the contrary, ETEC and VTEC were more frequently isolated from diarrheic samples, suggesting a possible association between these pathotypes and the emergence of disease in 1-weekold piglets.  Table 1.

Discussion
The aim of this work was to characterize the E. coli strains involved in porcine neonatal diarrhea outbreaks by genotyping and phenotyping isolates grown as pure cultures from both clinical cases of diarrheic and non-diarrheic healthy pen-mates.
Surprisingly, the general prevalence of strains encoding genes for virulence factors, as well as the presence of the typical E. coli pathotypes, was low in our pig population studied. Merely 26% (32/122) of the strains could be classified into a specific pathotype. EPEC was the most frequent one, representing 14% of the total isolates, followed by ETEC (12%) and VTEC (3%). The presence of EPEC is usually reported as a causative agent of diarrhea, principally in post-weaned but also in suckling pigs, but in this study, the presence of this pathotype was found more prevalent in non-diarrheic piglets (21.4%) than in the diarrheic ones (11.8%). Similar results have been reported in previous works in agreement with the results presented in the present study. Mainly, EPEC prevalence ranging from 3% to 26% in diarrheic animals and 9% to 17% in non-diarrheic animals has been described [26][27][28]. Finally, EAST1 gene was highly detected in both diarrheic (70%) and non-diarrheic groups (86%) but again with no relationship with the clinical onset. The EAST1 toxin is produced by several E. coli pathotypes [29], and although initially it was believed to be involved in the pathogenesis of diarrhea, most authors have concluded that this toxin may not be a virulence determinant, because it has been commonly isolated from healthy animals as well [30,31].
On the contrary, ETEC and VTEC were more frequently isolated from diarrheic samples, suggesting a possible association between these pathotypes and the emergence of disease in 1-week-old piglets.
The AMR phenotypes of the E. coli isolates showed a wide resistance to sulfonamides, tetracyclines and aminopenicillins in both diarrheic and non-diarrheic piglets. High percentages of antimicrobial resistance for amoxicillin and ceftiofur have been previously reported in pig farms from Catalonia [32]. However, other studies conducted in other regions of Spain [33] or in Ireland [34] reported lower levels of resistance to tetracycline, sulfamethoxazole/trimethoprim, ampicillin and streptomycin. These differences could be explained by the fact that Catalonia is one of the areas with greater density of pigs in Europe and therefore is likely subjected to a high consumption of antibiotics. On the other hand, an interesting finding of this study was that E. coli isolates from diarrheic cases additionally presented extended AMR to quinolones and aminoglycosides (gentamicin and neomycin) compared to those from healthy pigs. These findings agree with other studies performed in Ontario [35] and in Denmark [36], which found consistently higher frequencies of resistance in diarrheic and ETEC isolates, respectively.
The genotyping study showed that the E. coli strains containing the bla CTX-M gene presented resistance to third-and fourth-generation cephalosporins, amoxicillin/clavulanic acid, quinolones and aminoglycosides. The relationship between ESBLs and quinolone resistance mechanisms has been documented by several authors [37][38][39], and it has been suggested that ESBL genes and plasmid-mediated quinolone resistance are located within the same plasmid [40]. Furthermore, strains harboring the bla CTX-M gene presented higher values of phenotypical resistance than those without this gene, mainly to all the tested beta-lactams, quinolones and gentamicin.
Regarding the genotyping of the CTX-M+ strains, all the allelic variants described in this study have previously been reported in pigs, with CTX-M-14 also being the most prevalent [41,42] In the case of colistin, resistance was assessed using two different methods, the disk diffusion method and the MIC test, finding low concordance between both techniques. According to several authors, the disk diffusion method does not seem to be a reliable method for testing the efficacy of colistin due to the poor diffusion of polymyxins in agar. Therefore, the use of other tests, such as MIC, is more adequate for testing the susceptibility of colistin [43].
In general, the prevalence of porcine colistin-resistant E. coli strains reported in Europe is low: 0.9% in Poland [44], 1.5% in France [45], 2.9% in Germany [46] and 3.7% in Switzerland [47]. These results are in concordance with the 5% of prevalence obtained in this study using the MIC test. Nevertheless, García et al. (2018) have reported an extremely high prevalence of colistin resistance (77%) in the northwest area of Spain [48]. Variations in animal health management, biosecurity and sanitation procedures in pig production between different regions of the same country could explain this different colistin-resistance prevalence.
Likewise, the prevalence of E. coli isolates harboring the mcr-1 gene was low (3%) in this study, in agreement with results reported in Germany [46,49] and Belgium [50] with prevalence below 10%. Nevertheless, prevalences of 26% and 88% have been published in studies conducted in the north-west area of Spain [48] and in France [51], respectively. Only one E. coli isolate of this study presented the mcr-4 gene; however, the frequency of this gene in the Spanish pig farms is variable, ranging from 21% [8] to 73% [48] depending of the region. Regarding other mcr genes, none of the isolates of this study were positive for mcr-2, mcr-3 or mcr-5. The mcr-2 gene has been detected in porcine and bovine E. coli isolates from Belgium [50], mcr-3 in E. coli isolates from pigs in UK [52] and mcr-5 in 3.4% of isolates in Spain [48].
In summary, E. coli CTX-M+ strains isolated from diseased piglets showed antimicrobial multi-resistance to β-lactams, tetracyclines, sulfonamides and lincospectin extended to quinolones and aminoglycosides. This high prevalence of E. coli MDR strains in the pig population can represent a serious threat for animal and human health [53].

Conclusions
More appropriate diagnostic approaches, including genetic and phenotypic analysis of AMR profiles, should be implemented in animal health to optimize the use of antibiotics for treating diseased animals when necessary and with an effective antimicrobial agent, as the presence of highly resistant strains harboring ESBL and mcr genes circulating in pig farms can represent a potential threat to animal and human health.