α - L -Fucosidases from Bursaphelenchus xylophilus Secretome—Molecular Characterization and Their Possible Role in Breaking Down Plant Cell Walls

: The pinewood nematode (PWN) Bursaphelenchus xylophilus , the causal agent of the pine wilt disease (PWD), enters above-ground parts of the tree, migrates through the resin canals and feeds on plant cells causing extensive damage. In order to penetrate the cell wall and establish a parasitic relationship with host trees, the PWN needs to secretea mixture of active cell wall degrading enzymes. In maritime pine, Pinus pinaster , which is high susceptible to PWN, xyloglucan is the major hemicellulosic polysaccharide in primary cells. The xyloglucan backbone is susceptible to hydrolysis by numerous endoglucanases, some of them specific to xyloglucan. However, to completely degrade xyloglucan, all substitutions on the glucan backbones must be released, and L fucose residues in xyloglucan branches are released by α - L -fucosidases. In the present study, the molecular characterization of two α - L -fucosidases found in PWN secretome was performed. Moreover, a novel α - L -fucosidase was identified and its cDNA and gene sequence were determined. The three-dimensional structures of these α - L -fucosidases were predicted and the transcript levels were analyzed, thus providing new insights into fundamental PWN biology and the possible role of these proteins as cell wall degrading enzymes.


Introduction
The pinewood nematode (PWN), Bursaphelenchus xylphilus, the causal agent of the pine wilt disease (PWD), is recognized worldwide as a major forest pest. Unlike other plant-parasitic nematodes such as cyst, root knot or root lesion nematodes, which infect the plant roots and tubers, the PWN enters above-ground parts of the tree, migrates through the resin canals and feeds on plant cells. Therefore, the pine tree cell wall, which is mainly composed of cellulose (40%-50%), hemicellulose (≈25%) and lignin (25%-35%), constitutes the first barrier that PWN must overcome. Cell wall proteins and pectins are among the other minor compounds of the cell wall [1]. The secretion of active cell wall degrading enzymes (CWDE) is essential for PWN to invade this cell wall and establish a parasitic relationship with the host tree.
Cellulose is a homopolymer of β-linked glucose units forming glucan chains linked to each other by hydrogen bonds, which form the cellulose microfibril. Its hydrolysis is achieved by the synergistic action of cellulases. In the B. xylophilus genome, the 11 genes identified as cellulases are endo-β-1,4glucanases belonging to the glycoside hydrolase family 45 (GH45) and are not found in any other nematode genus in which the complete genome is known [2,3]. Seven of these cellulases are also found in B. xylophilus secretome [4,5].
The importance of fungal and bacteria α-L-fucosidases in plant polysaccharide degradation, namely, on xyloglucan has been reported previously [10,12]. A secreted α-L-fucosidase belonging to the GH29 family from the plant pathogenic fungus Fusarium graminearum was shown to be active in releasing fucose from a xyloglucan fragment [13]. Nevertheless, the importance of α-L-fucosidases from plant-parasitic nematodes and their possible role in breaking down plant cell walls have been underestimated. However, the inability of a nematode to degrade any component of the plant cell wall may result in unsuccessful infection, thus, all possible carbohydrate active enzymes that could act on plant cell wall degradation must be further explored. Recently, a α-L-fucosidase gene from B. xylophilus (Bxy-fuca) was characterized and associated with nematode development, lifespan and reproduction [14]. In our previous studies, two α-L-fucosidases were detected in B. xylophilus secretome [5] and in the present study, the molecular characterization of these two α-L-fucosidases was performed, thus providing new insights into the possible role of these proteins in degrading the plant cell walls and in PWN pathogenicity.

Pinewood Nematodes
Nematodes from a Portuguese B. xylophilus isolate (BxPt17AS) were maintained in cultures of the fungus Botrytis cinerea grown on malt extract agar medium at 25 °C. Mixed developmental nematode stages grown for 15 days on fungal cultures were collected with distilled water using a 20 µm sieve and washed three times with sterile water or RNase free water. These nematodes were then used for DNA and RNA extraction or treatment with P. pinaster extracts. Aqueous P. pinaster extracts from the stems of non-inoculated seedlings were prepared as previously described [5] and used to incubate the nematodes, simulating pine stimulus. Nematodes were soaked in 5 mL of pine extracts for 16 h at 28 °C, washed three times in 1x M9 buffer and concentrated via centrifugation. Nematode pellets of approximately 100 µL were used for RNA extraction.

Fucosidase cDNA and Genomic Coding Sequence
Bursaphelenchus xylophilus RNA was extracted from various stages of B. xylophilus (eggs, juveniles, females and males) as previously described [5] and DNase treated with TURBO DNA-free kit (Invitrogen, Carlsbad, CA, USA). The cDNA synthesis was carried out using iScript™ cDNA Synthesis Kit (BioRad, Foster City, CA, USA). DNA extraction was performed with the DNeasy Blood and Tissue Mini kit (Qiagen, Venlo, Netherlands), following the manufacturer's instructions. Synthetized cDNA and extracted genomic DNA were used for PCR amplification of selected fucosidase using primers BxFuF (TGCTTATCAGCTGCTGTCTT) and BxFuR (TGGCTGTTTAATGGCAAAAGGG). Primers were designed, based on transcript sequence (protein ID All_gs454_002563), to flank the predicted protein cDNA coding sequence fragment and the putative full-length protein-encoding genomic sequence, including partial sequences of 5´and 3´regions. Reactions were performed with GoTaq DNA polymerase (Promega, Fitchburg, WI, USA) in a Thermal Cycler (BioRad) with an initial denaturation step of 95 °C for 2.5 min followed by 35 reaction cycles of 95 °C for 30 s, annealing at 52 °C for 30 s and extension at 72 °C for 1.5 (cDNA template) or 2.5 min (DNA template) and a final extension at 72 °C for 5 min. Amplified products were purified using the Qiaquick PCR Purification Kit (Qiagen) and sequenced in both strands in an Automatic Sequencer 3730xl under BigDyeTM terminator cycling conditions at Macrogen, Inc. using specific primers (Supplementary Materials, Table S1).

Sequence Analyses
Sequence analyses were conducted using BioEdit [18]. To identify introns in the genomic sequence, Softberry FGENESH [19] was used and the gene scheme was generated using Exon-Intron Graphic maker [20] available from WormWeb.org. Molecular mass values and isoelectric points of deduced amino acid sequence were predicted using the deduced protein sequence by analysis with ProtParam tool software available at the ExPASy website [21]. InterPro [16] was used for the protein family search, predicting the presence of domains and important sites, and SignalP 4.1 Server [22] was used to predict the presence of the signal peptides. Deduced amino acid sequences were BLASTp against the NCBI non-redundant protein database [23] and UniProt Knowledgebase (reviewed) database (UniProtKB/Swiss-Prot) [24], where records were manually annotated with information extracted from literature and curator evaluated computational analysis. Sequence alignments were obtained using BioEdit between deduced amino acid sequences and predicted fucosidases found after BLASp searches and identified in plant-parasitic nematodes' available genomes, as described above. This alignment was used to estimate phylogenetic relationships in MEGA 6 [25] by the maximum likelihood method based on the Jones-Taylor-Thornton (JTT) matrix-based model.

Protein Homology Modeling
The 3-D structures were modeled in silico with the protein structure homology-modeling server SWISS-MODEL [26]. The best template for protein 3-D structure prediction was selected with the SWISS-MODEL template library search. The optimality of the predicted structures was estimated with global model quality estimation (GMQE) and qualitative model energy analysis (QMEAN) [27], both implemented in SWISS-MODEL. The resulting GMQE score was expressed as a number between 0 and 1, reflecting the expected accuracy of a model built with that alignment and template and the coverage of the target. Higher numbers indicate higher reliability, whereas for QMEAN, scores of −4.0 or below indicate that the predicted structures are very low quality. Additionally, the obtained structures were validated with a Ramachandran plot analysis implemented on the RAMPAGE server [28]. The structures with a higher number of residues in favored regions (>90%) and lower number of residues in outlier region were selected. The binding of possible ligands to predicted protein structures was previewed using the protein-small molecule docking web service SwissDock [29] and visualized using UCSF Chimera 1.4 software [30].

Fucosidases' Transcript Levels
The relative transcript abundance of fucosidases in B. xylophilus grown on fungus and after P. pinaster extract stimulus was assessed by reverse transcription quantitative real time PCR (RT-qPCR) with SsoAdvance Universal SYBR Green supermix (BioRad), according to standard protocols, using the CFX96 Touch™ Real-Time PCR Detection System (BioRad). The cDNAs of extracted RNAs were synthesized as previously described and used for qPCR. The amplification kinetics of each transcript was normalized with the amplification kinetics of the actin and 18S genes, chosen as endogenous controls and amplified with primers as previously described [5]. Other primers used in qPCR were designed based on B. xylophilus transcript sequences (Supplementary Table S1). qPCRs were done at 98 °C for 30 s, followed by 40 cycles of 98 °C for 5 s and 60 °C for 30 s. Melting curve analyses were performed and validation experiments were first carried out to ensure equivalent amplification efficiency for all transcripts. The RT-qPCRs were conducted for three biological repetitions, with three technical replicates for each qPCR. Amplification efficiencies and Ct values were determined by the CFX ManagerTM Software 2.1 (BioRad) and the mean Ct values were used for relative transcript level analysis (ΔΔCt method) with statistically significant differences tested using the pair-wise fixed reallocation randomization test© in REST software [31].
The BxFUCA1 protein sequence deduced from transcriptome and genome data, is 462 amino acids long with a calculated molecular mass of 52.9 KDa and a theoretical isoelectric point (pI) of 6.23. The protein sequence includes a N-terminal signal peptide of 15 amino acids with no predicted transmembrane domain. After removing this signal peptide sequence, a molecular mass of 51.3 KDa and pI of 6.25 were estimated.
Nucleotide sequences of 2332 bp long and 1491 bp long were determined for BxFUCA2, and submitted to Genbank under the accession numbers MN935185 and MN935184 for genomic amplicon and cDNA amplicon, respectively. The genomic sequence revealed the presence of eight introns ( Figure 2).  BxFUCA1 and BxFUCA2 predicted proteins share 67.4% of amino acid sequence identity. The top 10 BLASTp hits of each predicted protein against NCBI non-redundant protein database were hypothetical proteins and predicted α-L-fucosidases from different nematode species (Caenorhabditis elegans, C. nigoni, C. briggsae, Steinernema carpocapsae, Ancylostoma caninum, A. duodenale, A. ceylanicum, Nippostrongylus brasiliensis and Haemonchus contornus) with sequence identities of 49%-54%. The BLASTp search against UniProt reviewed database revealed BxFUCA1 and BxFUCA2 sequence identities of 46%-53% with 12 reviewed α-L-fucosidases, included in the phylogenetic analysis as reference α-L-fucosidases. Predicted α-L-fucosidases from plant-parasitic nematodes genome sequencing data, available on WormBase Parasite database, were also included in this study. Phylogenetic analysis exhibited a clear separation of fusosidases from animal origin with those of bacteria or fungi origin and a nematode-specific divergence inside animal origin fucosidases is denoted. Additionally, all predicted plant-parasitic nematode fucosidases, excluding B. xylophiulus fucosidases, form a separate cluster from other animal-parasitic and free-living nematodes. Considering B. xylophiulus fucosidases, phylogenetic analysis revealed that BxFUCA1 is identical to BXY_0325100.1 gene product and that BxFUCA2 differs from the three predicted proteins deduced from B. xylophilus genome, thus forming a separate cluster from all the other nematode α-Lfucosidases ( Figure 4). Predicted α-L-fucosidases from the bacteria Thermotoga maritima and the fungus Fusarium graminearum are also included. The percentage of trees in which the associated taxa clustered together is shown next to the branches. Initial tree(s) for the heuristic search were obtained by applying the neighbor-joining method to a matrix of pair-wise distances estimated using a Jones-Taylor-Thornton model. The tree is drawn to scale with branch lengths measured in the number of substitutions per site. The analysis involved 44 amino acid sequences. All positions with less than 50% site coverage were eliminated.

Structural Analyses
To better understand the structure and function of these proteins, 3-D structures were generated by homology modeling. BxFUCA1 and BxFUCA2 structures were predicted using the crystal structure of the GH29 family α-L-fucosidase from Fusarium graminearum (PDB ID: 4ni3.1.A) as a template. Both proteins presented two domains. For BxFUCA1, a C-terminal domain composed of seven anti-parallel β-sheets and a N-terminal domain formed by five parallel β-sheets arranged around the central axis and surrounded by eight α-helices forming the catalytic pocket, were predicted. The BxFUCA1 3-D structure is similar to the Bxy-FUCA structure predicted by Tang et al. [14]. For BxFUCA2, the C-terminal domain was composed of eight anti-parallel β-sheets and the Nterminal catalytic domain was formed by five parallel β-sheets arranged around the central axis and surrounded by eight α-helices (Figure 4). In addition to the GMQE and QMEAN values of 0.64 and −3.84 for BxFUCA1 and 0.64 and −3.63 for BxFUCA2 predicted structures, the Ramachandran plots validated these structural models with more than 90% of the residues located in the most favored regions (Supplementary Materials, Figure S1). Docking of a α-L-fucose molecule (CHEBI ID: 42548) to the catalytic pocket of BxFUCA1 structure was previewed with a full fitness of −2164.72 Kcal/mol and an estimated ΔG of −6.36 kcal/mol. For BxFUCA2, the docking of this fucose molecule to the catalytic pocket was also previewed with a full fitness of −2135.23 Kcal/mol and an estimated ΔG of −6.41 kcal/mol ( Figure 5).

Transcript Level Analysis
The relative transcript abundance of BxFUCA1 and BxFUCA2 in B. xylophilus grown on fungus and after P. pinaster extract stimulus was assessed by RT-qPCR. Both genes were found to be expressed in fungus and after pine extract stimulus, with no significant differences found between the two conditions after statistical analysis (p > 0.1). The transcript level of BxFUCA2 in B. xylophilus under pine tree extract was higher than BxFUCA1 ( Figure 6).

Discussion
Several studies have shown that plant-parasitic nematodes possess a set of synergistically active CWDE to allow their invasion and migration through the plant host tissues. The study of nematode CWDE not only enhance our understanding on the nematode-host interaction but also provides information on novel sources of enzymes with potential use in the biomass industry. There is a sustained focus on the discovery and application of carbohydrate-active enzymes for the selective deconstruction of plant biomass for the production of food, chemicals, biomaterials and biofuels [12].
A genome wide analysis of the CWDE from the genome sequenced phyto-parasitic nematode species identified 119 putative CWDE in B. xylophilus genome, most of them glycoside hydrolases. The number of GH genes and their family varies between analyzed species, indicating a possible evolution of plant-parasitic nematodes according to their feeding behavior [3]. Although not mentioned by the authors as predicted CWDE in published B. xylophilus genome [2,3] or even first published secretome [4], three genes and proteins with possible fucosidase activity belonging to GH29 were identified in B. xylophilus genome and secretome and were found listed in the supplementary tables of the respective publications. These proteins were identified as GH29 and from these, only the one coded by the BXY_0325100.1 gene (Bxy-fuca) was further characterized and described as playing a crucial role in the development, lifespan and reproduction of B. xylophilus [14]. The two proteins identified as α-L-fucosidases in our previous study on B. xylophilus secretome [5] were further characterized in this study, and a new α-L-fucosidase of B. xylophilus (BxFUCA2) was described.
The characterization of the complete cDNA coding sequences of BxFUCA1 and BxFUCA2 and analysis of their deduced amino acid sequences confirmed that these proteins are α-L-fucosidases belonging to the GH29 CAZy family, and they contain the typical domains and conserved sites of these proteins. A signal peptide was detected in predicted amino acid sequences, suggesting that these proteins are secreted, which is in accordance with previous data on B. xylophilus secretome [5]. Phylogenetic analysis of these proteins showed that fucosidases from B. xylophilus form a separated cluster, which is distant from the other plant-parasitic nematodes fucosidases, indicating a possible evolution of these proteins according to nematodes' host and feeding behavior. Unlike other plantparasitic nematodes, such as Globodera sp. or Meloidogyne sp. that infect the plant roots and/or tubers, B. xylophilus enters the above-ground parts of the pine tree, migrates through the resin canals and feeds on cells.
To date, GH29 crystal structures have been determined only for bacterial and fungi species and not yet for any eukaryotic species. The 3-D structure of BxFUCA1 and BxFUCA2 were predicted by homology-modeling with the available structure of the secreted α-L-fucosidase from the plant pathogenic fungus Fusarium graminearum. This fungal fucosidase was found to be active in releasing fucose from xyloglucan [13], the major hemicellulosic polysaccharide in pine cell walls [7]. A similar role for B. xylophilus-secreted BxFUCA1 and BxFUCA2 could be inferred, acting on the fucosylated side groups of xyloglucan, breaking down the self-aggregation of xyloglucan chains, the binding of xyloglucan to cellulose and exposing the xyloglucan backbone for further degradation. These proteins may also act on other fucosylated substrates like pectin or glycoproteins, contributing to deconstruction of the plant cell wall. The possible involvement of BxFUCA1 in hemicellulose digestion along with its crucial role in the development, lifespan and reproduction of B. xylophilus, has recently been suggested in a work describing the Bxy-fuca gene [14].
Our findings suggest that α-L-fucosidases belonging to the GH29 family, from B. xylophilus and other plant-parasitic nematodes, should be included in CWDE studies as this family of proteins is involved in many crucial biological processes that include fucosylated glycoconjugates, and it has been recognized as being involved in xyloglucan degradation. Xyloglucan constitutes a significant reserve of metabolically accessible monosaccharides for diverse phytopathogenic, saprophytic and gut symbiotic micro-organisms. To overcome the intrinsic stability of the diverse glycosidic bonds in this plant polysaccharide, bacteria and fungi have evolved extensive repertoires of xyloglucan-active enzymes [10] and plant-parasitic nematodes, like B. xylophilus, also had to do it.
The molecular characterization of the two α-L-fucosidases identified in B. xylophilus secretome increases our knowledge on this important group of enzymes, and provides new insights about the importance of these proteins in B. xylophilus pathogenicity. Moreover, the limited homology of these proteins with other known fucosidases makes them promising targets for the further development of target-specific strategies for PWN management.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Figure S1: Ramachandran analysis of predicted three-dimensional structure of two α-L-fucosidase from Bursaphelenchus xylophilus. Ramachandran plot of BxFUCA1 (a) and BxFUCA2 (b) predicted structures, Table S1: Primers used for BxFUCA2 cDNA and DNA coding sequence amplification, sequencing and qPCR and for BxFUCA1 qPCR.