Firstly, we obtained the whole genome of B. xylophilus (SAMEA7282713) (https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_904067135.1/) data from Wormbase (https://parasite.wormbase.org). To obtain a representative set of SULP orthologues, we first screened all publicly available nematode genomes deposited in NCBI database (https://www.ncbi.nlm.nih.gov/gene) (Howe et al., 2016) for the presence of complete, well-annotated SULP families. After applying strict criteria for sequence completeness and annotation quality, only three genomes (Brugia malayi, Caenorhabditis elegans, and Caenorhabditis briggsae) were retained. We used BLASTX to compare the queried protein sequence with the protein database of B. xylophilus to obtain the candidate BxSULPs sequence. BLAST searches were performed with an E-value cut-off of 1 × 10^-5, ≥ 30% query coverage and ≥ 25% identity. After that, we queried the Pfam database (https://pfam.xfam.org/) to search for the conserved structural domains possessed by the SULP protein and downloaded its Hidden Markov Model (El-Gebali et al., 2019). The protein sequences with the Hidden Markov Model were then searched from the protein database of B. xylophilus by HMMER3.0 software. HMMER 3.0 was run using the trusted cut-off (TC) defined in Pfam for the Sulfate_transp and STAS domains, with an inclusion threshold of E ≤ 1 × 10^–5 and a minimum alignment coverage of 70% of the Pfam model. Finally, the candidate sequences were characterized using NCBI-conserved domains (https://www.ncbi.nlm.nih.gov/Structure/cdd) (Lu et al., 2020b) and InterPro (https://www.ebi.ac.uk/Interpro) (Cantelli et al., 2021) to eliminate duplicates and acquire the protein sequences of BxSULPs. The physicochemical properties of BxSULPs were analyzed using Expasy 3.0 (https://web.expasy.org/protparam) (Duvaud et al., 2021).

In order to evaluate the extent of phylogenetic similarity among the 10 protein sequences, we conducted phylogenetic and conserved motif analyses of BxSULPs. Evolutionary relationships of BxSULP were analyzed using MEGA 11.0 software through multiple sequence alignments of 10 BxSULPs and the optimal model for the maximum likelihood tree was determined. Based on the prediction results of the maximum likelihood tree in the best model, the LG model was selected to construct the phylogenetic evolutionary tree of BxSULP. This model had a BIC (Bayesian Information Criterion) value of 11645.256 and an AICc (Corrected Akaike Information Criterion) value of 11362.515. The model was chosen because it provided the lowest BIC value (11645.256), indicating that it is the best fit for the data among the tested models (Luo et al., 2010). A maximum likelihood evolutionary tree was constructed with the Bootstrap test performed 1000 times (Koichiro et al., 2021). Subsequently, the constructed phylogenetic evolutionary tree was landscaped using iTOL (https://itol.embl.de/tree/) (Letunic and Bork, 2024). Gene information for Bx-sulps was obtained from the gene structure annotation file, and gene structure analysis was conducted using the GSDS 2.0 (http://gsds.gao-lab.org/) (Hu et al., 2015). Additionally, the distribution of conserved motifs in BxSULPs was examined using MEME Suite 5.5.6 (https://meme-suite.org/meme/) (Bailey et al., 2015) with default parameters and a motif value of 8. The results from MEME were combined with the maximum likelihood tree to create gene structure maps with annotated conserved motifs.

The gene density information of B. xylophilus from the annotation file of gene structure was extracted by Gene Density Profile function of TBtools. Then, we integrated the annotation file and the density file to map the distribution of BxSULP on chromosomes by Gene Location Visualize from GTF/GFF function (Chen et al., 2020, 2023). Subsequent analysis included predicting the secondary structure of BxSULPs using online tools like Prabi (https://doua.prabi.fr/software/cap3) (Thibaut et al., 2015), NetPhos 3.1 (https://services.healthtech.dtu.dk/services/NetPhos-3.1/) (Blom et al., 1999), and Protter 1.0 (http://wlab.ethz.ch/protter/start/) (Ulrich et al., 2014). Swiss-Model (https://swissmodel.expasy.org/interactive) (Marco et al., 2014) was employed to forecast the protein tertiary structure, and SAVES v6.1 (https://saves.mbi.ucla.edu) was used to create Ramachandran diagrams.

B. xylophilus was isolated from symptomatic Pinus thunbergii logs collected in Zhejiang Province, China. Xylem samples were cut into 2–3 cm thick disks and placed on a modified Baermann funnel (25 °C, 24 h) (Shinya et al., 2013). Nematodes that migrated into the collection tube were recovered and examined under a light microscope. Individual adult females were carefully picked out and identified to species according to the morphological diagnostic criteria (Abelleira et al., 2011). Only those specimens that exhibit the typical morphological characteristics of B. xylophilus are retained for subsequent cultures and experiments.

Synchronized eggs were obtained from embryos by incubating mixed-stage B. xylophilus in petri dishes at 25 °C in the absence of light for 10 minutes. This process facilitated the adhesion of eggs to the substrate through surface glycoproteins. The top layer containing water and nematodes was then removed to harvest the synchronized eggs. Subsequently, these eggs were allowed to hatch at 25 °C in the dark to yield second instar larvae (L2) in a nutrient-free setting. These synchronized L2 larvae were then transferred onto a Botrytis cinerea lawn on a PDA plate. Nematodes were harvested at 24, 48, and 72-hour intervals to obtain third instar larvae (L3), fourth instar larvae (L4), and adults, respectively. Dauer third-instar larvae (D3) and dauer fourth-instar larvae (D4) were isolated from infected pinewood by Baermann funnel. Male and female nematodes were manually sorted under microscope using worm pickers (Tang et al., 2020).

Gene expression in distinct developmental stages was analyzed by reverse-transcription quantitative polymerase chain reaction (RT-qPCR) by synchronizing B. xylophilus in different developmental stages. To explore the molecular response mechanisms of Bx-sulps to various treatments, expression levels of 10 Bx-sulps were determined across various developmental stages of B. xylophilus, including eggs. The statistical significance of the data was determined using a one-way analysis of variance (ANOVA) using GraphPad Prism 10.6.0.890 with a student’s t-test with p < 0.05 as the significance threshold.

Differentially expressed genes were chosen for validation via RT-qPCR (Tanaka et al., 2019). Primer sequences were generated with Primer 6.0 software and subsequently synthesized by Sangon Biotech (Shanghai, China) ( Table 1 ). Amplification was performed on a QuantStudio 1 Plus instrument (Thermo Fisher Scientific, USA) following the manufacturer’s instructions for the One Step PrimeScript™ RT-PCR Kit (RR064A, Beijing), and relative gene expression was calculated using the -ΔΔCt method. Gene-expression data were normalized to the Bx28S rRNA reference gene.

Following a BLASTX search in the Wormbase database and subsequent validation of conserved domains using Conserved Domains and InterPro, we identified 10 distinct SULP proteins in B. xylophilus, designated as BxSULP1 to BxSULP10 based according to their ascending order on chromosomes ( Table 2 ). To further elucidate the evolutionary relationships, we conducted a comparative analysis of the protein sequences of SULP found in BxSULP1 to BxSULP10 with those of analogous proteins from other nematodes.

Comparative analysis of BxSULPs with homologs from nematodes like C. briggsae and C. elegans revealed both conserved and divergent evolutionary patterns. High BLASTP scores (above 100) and low E-values (below 1×10^-10) indicated significant sequence similarity, implying shared conserved structural and functional domains, notably in the Sulfate Transporter and Anti-sigma factor antagonist (STAS) domain crucial for sulfate binding and transport. Notably, BxSULP8 exhibited a remarkable BLASTP score of 462 and an exceptionally low E-value of 8.4×10^–148 compared to a STAS domain-containing protein in C. elegans. This high degree of sequence similarity suggested that BxSULP8 was highly conserved and likely played a crucial role in sulfate transport, similar to its orthologs in other nematodes. Given the critical function of the STAS domain in sulfate binding and transport, BxSULP8 might be particularly important for maintaining sulfate homeostasis in B. xylophilus. This warrants further functional investigations to elucidate its specific role in nematode biology. These differences underscored the importance of understanding the evolutionary and functional contexts of SULP across diverse nematode species, providing insights for targeted potential control strategies and targeted interventions.

The average number of amino acids in BxSULPs was 202.2, with an average molecular mass of 22.71 kDa. The size of BxSULP varied, ranging from 17.63 kDa (BxSULP1) to 29.47 kDa (BxSULP2). Notably, BxSULP10, the largest protein, may exhibit more intricate functionalities or interactions. The isoelectric points (pI) of BxSULPs spanned from 6.51 (BxSULP6) to 10.19 (BxSULP9). Proteins with pI values near neutrality are anticipated to demonstrate enhanced stability under physiological conditions, whereas those with extreme pI values might serve specialized roles. BxSULP9 could participate in processes occurring in high pH environments. BxSULP6, characterized by the highest aliphatic index (119.51), may exhibit greater resistance to thermal denaturation compared to other BxSULPs, a critical feature for its potential function in extreme nematode environments. BxSULP2, BxSULP5, BxSULP9, and BxSULP10 were hydrophilic, while the remaining proteins were hydrophobic. BxSULP4, BxSULP7, and BxSULP8 demonstrated stability, contrasting with the instability observed in the other proteins. This suggested that protein may possess more robust structures, potentially ensuring sustained functionality over time ( Table 3 ).

We performed phylogenetic analysis on the full-length protein sequences. The conserved STAS and Sulfate_transp domains were subsequently mapped onto the alignment to verify motif conservation. The phylogenetic analysis revealed distinct clustering patterns among SULPs, indicating conserved and divergent evolutionary trajectories. The grouping of BxSULPs with homologs from closely related nematodes like C. briggsae and C. elegans suggested significant sequence conservation and shared ancestry. For example, BxSULP7 clustered closely with CbrSULP6 and CeSULP6, implying a strong evolutionary relationship, while functional equivalence remains to be experimentally determined. However, notable divergence was observed among SULPs from different nematode species. BxSULP2 exhibited a unique branching pattern, suggesting potential functional or adaptive divergence specific to B. xylophilus, possibly linked to its distinct ecological niche. The phylogenetic tree categorized BxSULP into distinct clades, denoted by different colors. Purple represents the SULP I family, red the SULP II family, and green the SULP III family, highlighting separate evolutionary lineages within the SULP family. SULP4, SULP5, SULP6, SULP8 of B. xylophilus, SULP3 of C. briggsae and SULP8 of B. malayi form a well-supported clade. Similarly, SULP1, SULP9 of B. xylophilus, SULP7 of C. elegans and SULP7 of C. briggsae clustered together in another distinct clade. SULP3, SULP7, SULP10 of B. xylophilus, SULP1, SULP2, SULP4, SULP5, SULP6, SULP8 of C. elegans, SULP2, SULP4, SULP5, SULP6, SULP8 of C. briggsae and SULP6 of B. malayi clustered in a separate branch, supporting the concept of gene duplication and subsequent functional diversification. These clades were also strongly supported by high bootstrap values, affirming their robustness ( Figure 1 ).

In the context of sulfate permeases, the motifs are likely involved in critical functions such as sulfate binding, transmembrane transport, and protein stability ( Supplementary Figure S1 ). Analysis determined the compositions of introns and exons in Bxsulps, revealing that the number of exons ranged from 7 to 16 across the 10 genes. Specifically, BxSULP3, BxSULP8 and BxSULP10 had only 7 exons, while BxSULP5 exhibited the highest number with 14 exons. Examination of structural motifs in the BxSULPs unveiled shared motifs among these proteins, as corroborated by the analysis of conserved motif structures ( Supplementary Figure S2 ). The figure illustrated the conserved motifs identified within BxSULPs, emphasizing regions of sequence similarity across the proteins. These motifs are short, conserved sequences that are likely involved in critical functions such as sulfate binding, transmembrane transport, and protein stability. The clustering of these motifs indicates regions of high conservation, suggesting potential functional roles.

The chromosomal localization of the 10 Bx-sulps in B. xylophilus was determined using TBtools v2.056, precisely mapping the locations of these channel-encoding genes across various chromosomes ( Figure 2 ). Analysis of chromosomal localization revealed that the genes encoding the 10 SULP proteins were distributed among five chromosomes of B. xylophilus. Specifically, on chromosome 1, Bx-sulp1 was situated in the 3′ region. Bx-sulp2 and Bx-sulp3 were positioned in the middle of chromosome 2. Furthermore, Bx-sulp4 and Bx-sulp5 were identified in the 3′ region of chromosome 3. Within chromosome 5, Bx-sulp6 and Bx-sulp7 were located in the central region, and Bx-sulp8 was found in the 3′ region. Additionally, Bx-sulp9 and Bx-sulp10 were situated in the 5′ region of chromosome 6, with no genes present on chromosome 4. The dispersion of Bx-sulps across multiple chromosomes indicates their non-clustered distribution throughout the genome. This dispersion had implications for gene regulation, as genes located on distinct chromosomes could be governed by diverse regulatory mechanisms and environmental factors.

Analysis of the three-dimensional structures of 10 BxSULPs revealed similar proportions of structural elements, including α-helices (32.56%–60.22%), β-folds (6.11%–23.12%), β-turns (0.60%–11.05%), and random coils (28.73%–53.82%). The transmembrane segment of BxSULPs predominantly comprised α-helices, facilitating their integration into the cell membrane to form channel structures for sulfate transport. Despite variations in the number and positioning of transmembrane segments across different proteins, a general conservation of these regions was observed. While BxSULPs exhibited sequence diversity, their three-dimensional structures displayed substantial similarity, indicating the preservation of core functional features essential for sulfate transport ( Table 4 ; Figure 3 ). The Ramachandran plot was used to assess the three-dimensional structures of 10 proteins. The dihedral angles of BxSULPs residues were found to be located in the yellow favored regions, and the spatial structure was found to be over 90% stable, indicating a high degree of confidence in the tertiary structure. The Ramachandran plot analysis confirms the high stability of the three-dimensional structures of BxSULPs, with the majority of residues falling within the core regions, indicating a high degree of confidence in the predicted protein conformations ( Figure 4 ).

The topological structure analysis indicated that BxSULP2, BxSULP3, BxSULP4, BxSULP6, BxSULP7, BxSULP8, and BxSULP10 were transmembrane proteins characterized by the presence of transmembrane domains enabling their integration into the cell membrane, primarily facilitating sulfate transmembrane transport. Conversely, BxSULP1, BxSULP5, and BxSULP9 were intra-membrane proteins. Although these proteins lacked transmembrane domains, they can be secreted outside the cell via signaling peptides. These signal peptides were short, hydrophobic regions at the N-terminus of the protein that guided the protein through the endoplasmic reticulum and into the extracellular space. These intra-membrane proteins operate within the extracellular matrix or engage in intercellular signaling upon secretion. Once secreted, these proteins can interact with other molecules in the extracellular matrix or participate in intercellular signaling pathways. Their functions might include binding to extracellular ligands, modulating cell-cell interactions, or influencing the local environment around the cell. Notably, BxSULP1 and BxSULP3 exhibited the highest number of N-glycosylation sites, with 3 sites each, enhancing their stability and functionality in membrane-associated or secreted contexts. This suggested a pivotal role in sulfate transport and interaction with the extracellular environment. The elongated signal peptide of BxSULP7 suggested more precise localization within the cell membrane, thereby enhancing its efficacy in sulfate transport. Despite variations in transmembrane structure, glycosylation patterns, and signal peptides among the BxSULPs, they exhibited a degree of structural similarity. These structural variations suggested functional partitioning. The multi-pass trans-membrane proteins were predicted to form the core sulfate-transport channels. In contrast, the secreted is formed that contain signal peptides might function as extracellular sulfate-binding or sulfate-donating molecules, modulate local sulfur availability or serving as signaling ligands during host invasion and diapause ( Figure 5 ).

To analyze gene expression in different developmental stages of B. xylophilus, we employed reverse-transcription quantitative polymerase chain reaction (RT-qPCR). By analyzing gene expression data from B. xylophilus under various nematode states and stress conditions, the expression data of 10 Bx-sulps were obtained. Significant variations in Bx-sulps expression across different nematode states ( Figure 6 ). Specifically, Bx-sulps exhibited low expression levels during egg incubation and fourth instar larva (L4). Only two genes, Bx-sulp3 and Bx-sulp5, were highly expressed in second instar larva (L2). Bx-sulp1, Bx-sulp2 and Bx-sulp3 were highly expressed in third instar larva (L3). Notably, Bx-sulp4, Bx-sulp8, Bx-sulp9, and Bx-sulp10 showed heightened expression in dauer third-instar larva (D3), crucial for maintaining essential physiological functions during diapause and enhancing resistance. The overexpression of Bx-sulps potentially bolstered resistance during diapause, with sulfate likely playing a role in cell membrane stabilization and enhancing cell tolerance to extreme conditions. Moreover, Bx-sulp1, Bx-sulp6 and Bx-sulp7 were prominently expressed in female adult nematodes and male adult nematodes, suggesting involvement in sex differentiation. While L2, L3, and L4 stages exhibited similar sulp genes expression patterns, the transition to adulthood, either male or female, was marked by significant changes in gene expression profiles.