Life-stage specific transcriptomes of a migratory endoparasitic plant nematode, Radopholus similis elucidate a different parasitic and life strategy of plant parasitic nematodes

Radopholus similis is an important migratory endoparasitic nematode, severely harms banana, citrus and many other commercial crops. Little is known about the molecular mechanism of infection and pathogenesis of R. similis. In this study, 64761 unigenes were generated from eggs, juveniles, females and males of R. similis. 11443 unigenes showed significant expression difference among these four life stages. Genes involved in host parasitism, anti-host defense and other biological processes were predicted. There were 86 and 102 putative genes coding for cell wall degrading enzymes and antioxidase respectively. The amount and type of putative parasitic-related genes reported in sedentary endoparasitic plant nematodes are variable from those of migratory parasitic nematodes on plant aerial portion. There were no sequences annotated to effectors in R. similis, involved in feeding site formation of sedentary endoparasites nematodes. This transcriptome data provides a new insight into the parasitic and pathogenic molecular mechanisms of the migratory endoparasitic nematodes. It also provides a broad idea for further research on R. similis.


RNA sequencing and transcripts annotation.
In this study, RNAs of eggs, juveniles, females and males of R. similis were sequenced using Illumina Hiseq 2000 which yields 62,440,260, 62,440,260, 62,440,260 and 62,440,260 reads respectively. 64761 unigenes were assembled in total (accumulated length of 129158666nt) with the average length of 1994nt, 3136nt for N50, which had 22,372 distinct singletons. Juveniles had the largest number of unigenes, followed by eggs, females and males (Table S2, Fig. S1). In total, 48432 unigenes were annotated by using NR, NT, Swiss-Prot, KEGG, COG, GO database respectively (Table S3). GO annotation of transcripts between different life stages have also been compared (Figs S2-S5).
Sequence similarity analysis. R. similis transcriptome data was blasted against the protein data of 9 nematodes including C. elegans, A. suum, B. pahangij, L. loa, M. incognita, M. floridensis, M. hapla, B. xylophilus and G. pallida. 8677 proteins from M. floridensis have similar sequences in R. similis transcriptome data, which were more than in other 8 nematodes (Fig. 1). 36 transcripts showing similarity with Wolbachia genes have been screened out. None of these sequences got the best hits with sequence in Wolbachia using NCBI online blastx (https://blast.ncbi.nlm.nih.gov/Blast.cgi) against NR database (Table S4). Other 744 transcripts were found to have similarity with genes from citrus (Table S5) and need further studies.
Transcriptional changes during R. similis life cycle. In this study, 11443 differentially expressed genes have been identified among eggs, Juveniles, females and males of R. similis (Fold change ≥ 5, FDR < 10 −5 ) (Fig. 2a). The gene expression model profiles of R. similis from eggs to females and males can be clustered into eight groups (Fig. 2b,c). There were three out of eight in egg to female and two out of eight from egg to male significantly enriched expression profiles (p < 0.05) respectively.
There were 2613 genes up-regulated in juveniles compared to eggs. These genes were mainly related to immunity, digestion and infection (Fig. S6a). In addition, significantly gene up-regulations can also be found in the gluconeogenesis, insulin-like pathway and steroid hormone pathway. It has been proved that these pathways were involved in the dauer formation of C. elegans 23,24 . There were 3546 genes down-regulated, which were mainly related to metabolism, growth, proliferation, transcription and protein synthesis (Fig. S6b). These genes expression patterns are consistent with the process from egg to an infective juvenile stage. 1350 genes were up-regulated in females compared to juveniles, and these genes were mainly involved in protein synthesis, immunology, sex differentiation and egg laying processes (Fig. S6c). In oocyte-forming and embryonic development-related metabolic pathway can also found up-regulated genes, which can be explained as reproductive organs development and early embryonic development in females. There were 5649 differentially expressed genes during juveniles developed into males, while, 3214 genes were up-regulated in males. Most of these up-regulated genes were www.nature.com/scientificreports www.nature.com/scientificreports/ associated with signal transduction and energy metabolism (Fig. S6e). In contrast, 2435 down-regulated genes were significantly enriched in xenobiotics metabolism and infection pathway (Fig. S6f).
Identifying plant parasitism genes. CAZYmes and cell wall degradation. Based on Carbohydrate-Active enZYmes Database, 4221 unigenes were assigned to five families of CAZYmes (File S2) using CAZymes Analysis Toolkit. 86 unigenes were found to have sequence similarity with cell wall degrading enzymes (Table S6), including fifteen sequences that could have been acquired by HGT (Table S6) and contamination from bacteria was eliminated. Among these 86 putative genes coding for cell wall modified proteins, 52 genes belong to GH5 family. 23 genes were grouped into other 6 cell wall degrading enzymes families: five genes to GH16 family, eight genes to GH30 family, two genes to GH43 and GH53 respectively, one gene to GH28 family and five genes to polysaccharide lyase. In addition, eleven genes having similarity with expansin-like protein coding genes were also found in R. similis transcriptome (Table 1). Among the 52 genes belonging to GH5 family, 44 genes could be divided into six clusters according to their sequence similarity.
The types of cell wall modified enzymes predicted in R. similis were much more than in other six PPNs. There was also a difference in the type of cell wall modified enzymes between R. similis and Aphelenchidae nematodes, B. xylophilus, A. besseyi and A. ritzemabosi ( Table 1). The cellulases predicted in R. similis belong to GH5 family, whereas most of predicted cellulases in Aphelenchidae nematodes belong to GH45 family. Only four GH5 family cellulases were predicted in A. ritzemabosi. However, pectin lyase was predicted only in B. xylophilus among these three Aphelenchidae nematodes. Expression analysis revealed that all putative genes coding for cell wall modifying protein were low expressed in eggs (Fig. 3). Most of them got the highest expression level in juveniles or  www.nature.com/scientificreports www.nature.com/scientificreports/ females than in other stages. It's worth noting that one predicted beta-1,4-endoglucanase (CL4657.Contig2) got the highest expression level in males. The expression pattern for this gene has been confirmed by qPCR.
Nematode resistance to oxidative stress. 102 genes were identified as antioxidant of nematodes, which can be divided into glutathione S-transferase (GST), superoxide dismutase (SOD), glutathione peroxidase and peroxiredoxin (Tables 2, S7). There were 47 sequences hit to glutathione S-transferase (GST). These predicted GSTs could be divided into eight classes: Delta, Zeta, Theta, Kappa, Sigma, Omega, Phi and Tau, most of them (n = 29) were grouped into Sigma GSTs. The expression patterns of these genes in R. similis were inconsistent. Each stage has several genes with the highest expression level (Fig. S7). A total of 27 genes coding for superoxide dismutase (SOD) were predicted and could be grouped into two classes: Cu-Zn SOD (n = 21) and Fe-Mn SOD (n = 6). There were eleven putative genes up-regulated in juveniles and sixteen putative genes up-regulated in adults (Fig. S7b). There were two and eight sequences hit to SODs of fungi and plant, respectively (Table S7). In the meantime, ten and eighteen sequences annotated to peroxiredoxin and glutathione peroxidase respectively. There were eight and two putative genes coding for peroxiredoxin were up-regulated in juveniles and females respectively (Fig. S7c). For genes coding for glutathione peroxidase, ten and eight were up-regulated in juveniles and adults respectively (Fig. S7d).
The species and quantity of predicted antioxidant in R. similis were more similar to those in P. zeae, than C. elegans and other four PPNs (M. incognita, G. pallida, A. bessey and B. xylophilus). Among all antioxidant, GST was the most abundant in these seven nematodes except for M. incognita. In addition, the migratory plant parasitic nematodes seemed to have more antioxidant than the sedentary plant parasitic nematodes (Table 2).  www.nature.com/scientificreports www.nature.com/scientificreports/ Prediction of other parasitic related genes in R. similis. 109 candidate genes could be classified in nine kinds of previously characterized parasitism genes of PPNs 25 (Table 3, File S2). These parasitism genes include: calreticulin, chorismate mutase, fatty acid and retinol binding protein (FAR), 14-3-3 and SPRY domain-containing protein that affect the host signaling pathway and inhibit host defense response, ubiquitin extension protein that regulate the host ubiquitination process and proteases such as Cathepsin 25 . Among these putative parasitism genes in R. similis, cathepsin is the most abundant (n = 50). Besides, parasitic related genes predicted in R. similis were much more than in P. coffeae, P. thornei and P. zeae, respectively.

Prediction of R. similis chemosensation genes.
Genes coding for G proteins, signal transduction proteins and regulators in chemosensation process were predicted in R. similis (Table 4, File S2). Only four G proteins, GCY-35, GPA-2, GPA-3, GPA-5 and GPA-11, could found similar sequences in R. similis transcriptome (Table 4). Six putative signal transduction proteins have been predicted in R. similis, including odr-1 and daf-11. odr-1 expressed in AWC, AWB neurons which could regulate olfactory responses in C. elegans. daf-11 can be expressed on both AWC, AWB and ASK neurons; it can also be expressed in ASI and ASJ and functions as a regulator of dauer formation and recovery 26,27 . Meanwhile, R. similis transcriptome were aligned with twelve chemosensation regulators in C. elegans, and only adp-1 and gpc-1 did not have any similar sequences in R. similis.

Genes involved in RNAi pathway.
To predict genes involved in the RNAi process, transcriptome data of R. similis was blasted against 72 proteins involved in C.elegans RNAi process (File S2). The result showed that putative RNAi genes were annotated to endogenous pri-miRNA formation and exportation genes drsh-1, Dicer complex and RISC (RNA-induced silencing complex) component genes dcr-1, drh-1, siRNA amplification genes ego-1, mut-7, RNAi intercellular spreading gene rsd-3 and RNAi negative regulation genes eri-1/6/7. No sequence was found to have similarity with sid-1/2, which facilitate exogenously double-stranded RNA enter cells.
Maybe, there were only few genes associated with small RNA synthesis, dsRNA uptake and spreading, RNAi inhibitors in R. similis compared to animal parasitic nematodes including A. suum and B. malayi. However, the number of putative nuclear RNAi effectors in R. similis was similar to A. suum and B. malayi. Moreover, R. similis

Validation of gene expression profiles by qRT-PCR.
To elucidate the expression profiles of this sequencing, nineteen differentially-expressed genes were chosen for qRT-PCR. The stability of the reference gene was test (Table S9). These genes were from various functional categories including cell wall degradation, antioxidant and biological progress etc. (Table 5). Among these genes, thirteen genes shown significantly higher expression level in juveniles than other life stages (p < 0.05), which were annotated as cell wall degradation genes, FAR and two function unknown protein. Two genes shown significantly higher expression levels (p < 0.05) in males than other life stages were predicted to be involved in cell wall degradation and antioxidant. Three genes were highly expressed in eggs which annotated as DNA-binding protein, calreticulin and SPRY domain containing protein. Furthermore, one putative gene coding for serine carboxypeptidase has also been confirmed to have a significantly higher expression level in females (p < 0.05) than in other life stages (Fig. 4). This independent experiment and statistical analyzed data revealed the reliability of the sequencing data.

Discussion
R. similis is an important migratory endoparasitic nematode, the morphology and biological characteristics of this nematode in different developmental stages were different. This is the first study on transcriptome of R. similis in eggs, juveniles, females and males separately. The number of unigenes obtained from juveniles sequencing were highest, followed by eggs, females and males. Total 64761 unigenes were obtained from R. similis transcriptome including 11443 differentially expressed genes. The expression patterns of differentially expressed genes were consistent with the biological functions of each life stage, and reflect the changes in the morphological and infectious ability of this nematode. The infection-related genes were up-regulated at juvenile stage and down-regulated at male stage. The genes associated with genital development such as gonads were up-regulated at female stage.
Transcriptome sequencing has been used in effectors research of PPNs 9,28,29 . These effectors include cell wall degrading enzymes, antioxidant enzymes and other proteins involved in resistance to the host defense responses. Many genes coding for cell wall modified enzymes were predicted in R. similis, as compared to other nematodes especially in the GH5 family. There were six clusters containing 44 sequences in the predicted GH5 family. It is speculated that there are lots of alternative splicing sites in the GH5 family genes of R. similis. As one kind of migratory endoparasitic nematode, R. similis can infect more than 250 plant species. To loosen cell walls of different plants and migrate in plant tissue, R. similis have evolved a variety of cell wall degrading enzymes. Only one predicted gene, coding for cell wall modified enzyme was highly expressed in males indicted that males also needed to degrade the cell wall in order to facilitate its movement in host tissue. R. similis was under the pressure of the reactive oxygen species produced by the host all the time. So, they must synthesize peroxidase to neutralize ROS in host plant 28 . All putative genes coding for SODs, glutathione peroxidases, peroxiredoxins and part of genes coding GSTs were highly expressed in juveniles and adults. These expression patterns were consistent with the biological characteristics of R. similis. It was reported that a glutathione peroxidase with signal peptide was secreted to the body surface of nematodes against the host's ROS. This gene was expressed during the whole infective stages of G. rostochiensis 29 . It can be speculated that these genes predicted in R. similis may also be involved  Table 5. Differentially-expressed genes used in qPCR validation.
www.nature.com/scientificreports www.nature.com/scientificreports/ in the resistance to the host ROS process. Besides, some genes coding for GSTs, highly expressed in eggs of R. similis, may play a role in xenobiotics metabolism 30 . Six kinds of putative effector genes have also been identified in R. similis. These effectors have been reported previously to play a role in host signaling regulation, host defense response inhibition and degradation of the host protein process [31][32][33][34] . Many sequences are identified as genes coding for cathepsin, which play a role in degrading host cell wall proteins and host guard proteins, predigesting proteins before uptake or protein digestion in the body 25 .
By comparing the types and numbers of pathogenic genes in the PPNs of different parasitism types, it can be confirmed that the parasitic modes of the nematodes are closely related to the type and expression pattern of their pathogenic genes. For R. similis, the number of putative genes coding for cell wall modified enzymes, peroxidases and proteases was higher than that of the sedentary endoparasite nematodes. On the other hand, no gene identified as effector required for giant cell or syncytium formation induced by root knot and cyst nematode such as 7E12 35 CEPs CLE peptide 36 and 16D10 CLE 37 were found in R. similis. Similarly, no sequences similar to MAP-1 involved in host recognition of root-knot nematode 38 were found. In addition, no sequences similar of disrupt host defenses effectors as 10A06, Hs19CO7 and 30C02 39-41 from cyst nematodes were found in R. similis. This nematode probably don't have these effectors which have been reported in sedentary endoparasitic nematodes. The difference between R. similis and sedentary endoparasite nematodes also exist in the expression patterns of pathogenic genes. Predicted pathogenic genes in R. similis were highly expressed in both juveniles and females, while these genes in sedentary endoparasite nematodes mainly expressed highly in juveniles 28,29,40 . All these differences between R. similis and sedentary endoparasitic nematodes may be related to the different parasitic strategies. As a migratory parasitic nematode, R. similis move and feed in the hosts. It is always necessary to degrade host cell walls which directly lead to the necrosis of the host cell rather than altering the structure of the host cell and maintaining the cell activity as sedentary endoparasite nematodes do. R. similis also needs to face more oxidative stress from the host during migrating and feeding, and therefore require more peroxidases. The similar phenomenon has been reported previously in P. zeae 42 .
No putative genes coding for cell wall modified enzymes in GH45 family was found in R. similis and other nematodes parasitizing plant underground tissue. At present, GH45 family genes are only reported in Aphelenchoidea nematodes for PPNs 43 . It is supposed that the GH45 genes in the nematode may be horizontally transferred from fungus 44 . Nematodes parasitizing aboveground parts of plants can feed on fungi to complete their life cycle. Therefore, some species of these nematodes may acquire GH45 genes during long-term feeding on fungi. While, R. similis and other nematodes parasitizing plant underground tissue can only feed on live plant tissue, therefore these nematodes have no GH45 gene. In addition, antioxidant enzymes predicted in R. similis were more than in A. besseyi and B. xylophilus, which may be due to the complex soil environment for R. similis to survive.
C. elegans larvae will stop developing and form dauer larva to increase the adaptability to the environment when it is perceived the large number of populations, the insufficient food supply, or the unsuitable ambient temperature 24 . The formation and recovery of dauer larva is regulated by signal transduction pathways such as guanylyl cyclase pathway, TGFβ-like pathway, Insulin-like pathway and steroid hormone pathway. At the same time, change has also been found in the energy acquisition in the dauer larva stage of C. elegans. The metabolism of dauer larvae is adapted to utilize internal energy reserves, predominantly lipid in the form of triglycerides, but also glycogen 45 . Dauer larva in R. similis has not been reported. However, R. similis could survive for six to eleven months without host 46,47 . When infect host, R. similis complete its life cycle in 18 to 20 days at 24 to 27 °C 48 . In this study, we found that genes related to dauer larva energy metabolism in the gluconeogenesis pathway, and the genes involved in signal transduction of insulin-like pathway and steroid hormone pathway were up-regulated www.nature.com/scientificreports www.nature.com/scientificreports/ in the juvenile. This results indicated that R. similis has the basis of signal regulation and energy metabolism of develops into a dauer-like larva. It is speculated that the development of juveniles to adults stop when the external environment is not suitable for their survival and develop into dauer-like larva to overcome unsuitable environment.
By comparing the predicted RNAi related genes between R. similis and other nematodes, we found that some genes involved in the C. elegans RNAi pathway are conserved in nematodes such as smg-2/6 for dsRNA uptake and spreading, xpo-1/2 for endogenous pre-miRNAs formation and transportation, and xrn-2 for RNAi inhibition 10 . In C. elegans, SID-1/2 is a key transmembrane protein for exogenous dsRNA and Exo-siRNA entry into cells 49,50 and RNAi delivery between cells is required with the assistance of RSD-2/3/6 51 . However, none of these genes got similar sequence in R. similis. In addition, the number of predicted RISC genes in R. similis was less than in C. elegans. We did not find all genes that were important for entering and spreading of dsRNA. But RNAi can also been observed after treated with genes specific dsRNAs 7,34 , suggesting that alternative proteins or similar dsRNA propagation mechanisms or a class of effectors that play a role in both endogenous and exogenous RNAi 52 . R. similis may also have the similar mechanism with root-knot nematodes.
HGT from bacteria and fungi to PPNs are quite common [53][54][55][56][57][58][59] . In this study, 15 predicted cell wall degrading genes horizontal transfer from fungi and bacteria have been found. We also found that 744 transcripts have sequence similarity with genes from citrus. So far, only three such candidate genes from host plant been reported in G. rostochiensis 60 . Only some candidate genes giving that expression have been verified in this study and other genes deserve further investigation.
This study, first time revealed the gene expression of R. similis, excavated a variety of pathogenic genes. The study is crucially important for examining the parasitic and pathogenic molecular mechanisms of migratory endoparasitic plant nematodes (such as R. similis), and their interaction with hosts.