Characterizing genetic variation on the Z chromosome in Schistosoma japonicum reveals host-parasite co-evolution

Abstract

Background: Schistosomiasis is a neglected tropical disease that afflicts millions of people worldwide; it is caused by Schistosoma, the only dioecious flukes with ZW systems. Schistosoma japonicum is endemic to Asia; the Z chromosome of S. japonicum comprises one-quarter of the entire genome. Detection of positive selection using resequencing data to understand adaptive evolution has been applied to a variety of pathogens, including S. japonicum. However, the contribution of the Z chromosome to evolution and adaptation is often neglected.

Methods: We obtained 1,077,526 high-quality SNPs on the Z chromosome in 72 S. japonicum using re-sequencing data publicly. To examine the faster Z effect, we compared the sequence divergence of S. japonicum with two closely related species, Schistosoma haematobium and S. mansoni. Genetic diversity was compared between the Z chromosome and autosomes in S. japonicum by calculating the nucleotide diversity (π) and Dxy values. Population structure was also assessed based on PCA and structure analysis. Besides, we employed multiple methods including Tajima's D, F_ST, iHS, XP-EHH, and CMS to detect positive selection signals on the Z chromosome. Further RNAi knockdown experiments were performed to investigate the potential biological functions of the candidate genes.

Results: Our study found that the Z chromosome of S. japonicum showed faster evolution and more pronounced genetic divergence than autosomes, although the effect may be smaller than the variation among genes. Compared with autosomes, the Z chromosome in S. japonicum had a more pronounced genetic divergence of sub-populations. Notably, we identified a set of candidate genes associated with host-parasite co-evolution. In particular, LCAT exhibited significant selection signals within the Taiwan population. Further RNA interference experiments suggested that LCAT is necessary for S. japonicum survival and propagation in the definitive host. In addition, we identified several genes related to the specificity of the intermediate host in the C-M population, including Rab6 and VCP, which are involved in adaptive immune evasion to the host.

Conclusions: Our study provides valuable insights into the adaptive evolution of the Z chromosome in S. japonicum and further advances our understanding of the co-evolution of this medically important parasite and its hosts.

Keywords: Schistosoma japonicum; Adaptive evolution; Genetic diversity; Z chromosome.

Fig. 1

Population genetic parameters across the genomes of Schistosoma japonicum to detect the supporting evidence for faster and more adaptive Z chromosome evolution. Mann-Whitney test was used for significant differences between the Z and autosomes denoted by *P < 0.05, **P < 0.01, ***P < 0.001. A Sequence divergence was compared in a pairwise manner between S evidence japonicum and S. haematobium and between S. japonicum and S. mansoni. B Adaptive substitutions rate (α) across the genomes of S. japonicum. Compare the adaptive substitution rate of Z chromosome and autosomes using all S. japonicum individuals (left of dash). Compare the adaptive substitution rate of the Z chromosome and autosomes for each population of S. japonicum (right of dash). C Nucleotide diversity (π) in five S. japonicum sub-populations calculated based on variation of the Z chromosome and autosomes. D Boxplot of absolute divergence (Dxy) between TW and other subpopulations. The dotted line represents a Dxy value of 0.4. E Heatmap of absolute divergence (Dxy) between each pair of subgroups for the sex chromosome (Z chromosome). The darker the red color, the higher the value

Fig. 2

Genetic diversity and population structure of sampled Schistosoma japonicum. A Map of sample collection for S. japonicum. B Principal component analysis (PCA) of S. japonicum male adults. The graph on the left is based on the genetic variation from all individuals’ Z chromosomes; the graph on the right is based on other individuals except TW (10) and PH (14). C Population structure of 72 S. japonicum individuals based on their Z chromosome variants. D Phylogenetic tree constructed based on the Z chromosome variations of 72 S. japonicum samples

Fig. 3

Positively selected gene DYS in TW population identified by CMS and iHS. A and B Positively selected signatures of DYS (Sj3555) from the CMS and iHS analysis. The dashed lines represent the empirical thresholds for the selected region (top 0.1% empirical distribution of P-value of iHS test = 5.0169; top 1% empirical distribution of CMS score = 12.831). The genome region harboring DYS (dystrophin) is marked in the figure. C Haplotype heatmap for SNP variants within the DYS gene region. Two nonsynonymous variants are marked with arrows. D F_ST distribution in the DYS gene region (ChrZ: 50351736–50406101). Five nonsynonymous variants within the TW version DYS are highlighted in red. E Alternative allele frequency for three non-synonymous mutations in the TW and C-L populations. F NJ phylogenetic tree based on the dystrophin full-length or C-terminal protein sequence. Caenorhabditis elegans was used as the outgroup

Fig. 4

Candidate gene LCAT may be related to the slower development of Schistosoma japonicum in Taiwan compared with other regions. A Positively selected gene LCAT (Sj2846) identified by CMS analysis. The dashed lines represent the empirical thresholds for the selected region (top 1% empirical distribution of CMS score = 12.831). The genome region harboring the gene LCAT (lecithin-cholesterol acyltransferase) is marked in this figure. B Haplotype network based on the SNPs in the LCAT region (ChrZ:24,935,941–24,982,741). Each circle represents a haplotype, and its size suggests the number of individuals harboring the haplotype. C Kernel density distribution of the PBS statistic (top 1% value = 2.756928) for the entire Z chromosome in the TW sub-population. The dotted line marks an SNP site in LCAT (ChrZ-24960821). D Extended haplotype decay around the LCAT-ChrZ-24,960,821 allele in the TW and C-L populations. E Relative mRNA expression levels of LCAT in females and males in developmental stages after infection of the definitive host. F Relative mRNA expression levels of LCAT in RNAi-treated parasites were analyzed by qPCR (mean ± standard error). GFP was used as the control group. Three biological replicates were performed. *P < 0.05, **P < 0.01, ***P < 0.001. G RNAi of LCAT causes parasite hypercontraction, Scale bar, 2000 μm. H Reproductive organs from LCAT and GFP RNAi parasites under confocal laser scanning microscopy. O, ovary; V, vitelline gland; T, testis. Three biological replicates were performed

Fig. 5

The candidate gene Rab6 may be related to the compatibility of intermediate hosts of Schistosoma japonicum in the C-M strain. There are 13 samples from the C-M populations and 12 samples from the C-L populations. A Positively selected signatures identified by XP-EHH (between C-M and C-L populations; top 0.5% value = 2.64) and iHS (within the C-M population; top 0.1% P = 3.76). The dashed lines represent the empirical thresholds for the selected region. The candidate window harboring gene Rab6 is highlighted in red. There are 13 samples from the C-M populations and 12 samples from the C-L populations. B Zoomed view of gene Rab6 and a heatmap of linkage disequilibrium (LD) measured by the squared Pearson’s correlation coefficient (r²) for SNP variants. C Haplotype network based on 97 SNPs in the Rab6 region. Each circle represents a haplotype, and its size suggests the number of individuals harboring the haplotype. D The relative mRNA expression levels of Rab6 in the four larval stages, including egg, miracidium, sporocyst, and cercaria. E The relative mRNA expression levels of Rab6 in RNAi-treated parasites were analyzed by qPCR (mean ± standard error). GFP was used as the control. Three biological replicates were performed. F RNAi of RAb6 causes parasite hypercontraction, with the GFP-treated group as a control. Scale bar, 2000 μm. G The morphology of the control group and experimental group intestines under confocal laser scanning microscopy. Scale bar, 25 μm. The intestinal region has been highlighted with red arrows; G, gut. Three biological replicates were performed

Fig. 6

The candidate gene VCP identified by F_ST analysis may be related to the host immune evasion of Schistosoma japonicum. There are 13 samples from the C-M populations and 12 samples from the C-L populations. A Positively selected signatures identified by F_ST between the C-M and C-L populations. The dashed lines represent the empirical thresholds for the selected region. The candidate window harboring gene VCP is highlighted in red. B Haplotype network based on 27 SNPs in the VCP gene region (there are 27 SNPs in the whole VCP gene region of six populations). Each circle represents a haplotype, and its size suggests the number of individuals harboring the haplotype. C The relative mRNA expression levels of VCP in RNAi-treated parasites, with GFP as the control group, were analyzed by qPCR (mean ± standard error). Three biological replicates were performed. D RNAi of VCP causes parasite hypercontraction. The GFP-treated group was used as a control. Scale bar, 2000 μm. E Reproductive organs from VCP and GFP RNAi parasites under confocal laser scanning microscopy. O, ovary; T, testis. Three biological replicates were performed. F Edu labeling showing the expression of EdU⁺ proliferative cells (pink) in GFP (RNAi) and VCP (RNAi) parasites. Three biological replicates. G The relative mRNA expression levels of VCP in the four larval stages, including egg, miracidium, sporocyst, and cercaria