Updated annotation and meta-analysis of Brugia malayi transcriptomics data reveals consistent transcriptional profiles across time and space with some study-specific differences in adult female worm transcriptional profiles

Abstract

Genomics, transcriptomics, and proteomics have significantly advanced our understanding of obligately host-associated microbes, where interrogation of the biology is often limited by the complexity of the biological system and limited tools. This includes the causative agents of many neglected tropical diseases, including filarial nematodes. Therefore, numerous transcriptomics studies have been undertaken on filarial nematodes. Most of these transcriptomics studies focus on Brugia malayi, which causes lymphatic filariasis and is a laboratory model for human filarial disease. Here, we undertook a meta-analysis of the publicly available B. malayi transcriptomics data enabling the direct cross comparison of samples from almost a dozen studies. This reanalysis highlights the consistency of transcriptomics results across many different studies and experimental designs from across the globe for over a decade of research, across many different generations of a sequencing technology, library preparation protocols, and differential expression tools. Males and microfilariae across samples had similar expression profiles. However, female samples were clustered into two differential expression patterns that were significantly different from one another. Largely, we confirm previous results for all studies reanalyzed including tissue-specific gene expression and anti-Wolbachia doxycycline treatment of microfilaria. However, we did not detect previously reported differential expression upon in vitro or in vivo treatment with ivermectin, albendazole, and DEC, instead identifying a consistent lack of transcriptomic change upon exposure to these anthelminthic drugs. Updated annotation has been provided that denotes poorly supported genes including those overlapping rRNAs.

Fig 1. Heatmap of the 9,727 Differentially Expressed Genes in the Meta-Analysis.

The dendrogram at the top of the heatmap was generated using pvclust. The red values are approximately unbiased (au) and the green are the bootstrap support (bp) values, both of which are generated by pvclust. The values shown are for illustrative purposes with the full dendrogram available in S3 Fig. The histogram at the bottom shows the distribution of all the z-score values in the heatmap. The heatmap uses a z-score normalization of log₂(TPM) values for 9,727 differentially expressed genes between 237 samples reanalyzed from 12 projects. The legend at the top is broken into three sections: project color, if the sample was drug treated, and sample life stage. The left hand legend is broken into two sections: the outer section denotes the WGCNA cluster and the inner section denotes if cluster matches the main or inverse WGCNA cluster expression pattern. Samples are labeled with first author, title, and bioproject from [–29].

Fig 2. Differential Expression of Soma and Germline Samples.

The dendrogram at the top of the heatmap was generated using pvclust. The red values are approximately unbiased (au) and the green are the bootstrap support (bp) values, both of which are generated by pvclust. The histogram at the bottom shows the distribution of all the z-score values in the heatmap. There are 7,191 differentially expressed genes between the proliferative zone (germline), meiotic zone (germline), and body wall (soma) samples. The heatmap shows a z-score normalization of log₂(TPM) values for the 7,191 differentially expressed genes divided between three WGCNA clusters. The top legend denotes the site of tissue origin. The left hand legend is broken into two sections: the outer section denotes the WGCNA cluster and the inner section denotes if cluster matches the main (grey) or inverse (black) WGCNA cluster expression pattern.

Fig 3. Differential Expression of Tetracycline Treatment.

The dendrogram at the top of the heatmap was generated using pvclust. The red values are approximately unbiased (au) and the green are the bootstrap support (bp) values, both of which are generated by pvclust. The histogram at the bottom shows the distribution of all the z-score values in the heatmap. Between the two treatment groups, there are 1,423 DE genes sorted into two WGCNA clusters. A z-score normalization of log₂(TPM) values was used for the heatmap. The top legend denotes the treatment groups. The left hand legend is broken into two sections: the outer section denotes the WGCNA cluster and the inner section denotes if cluster matches the main (grey) or inverse (black) WGCNA cluster expression pattern.

Fig 4. Lack of Differential Expression Upon Ivermectin Treatment.

The linear mixed model and PCA plots show a lack of effect from Ivermectin treatment [19]. The linear mixed model was generated using all counts. The PCA plots were generated using z-score normalized log₂(TPM) values for genes where the sum of the TPM values for each gene is greater than 0. (A) linear mixed model for study I (PRJNA303987) (B) linear mixed model (PRJNA303986). (C) Principal Components Analysis for Study I (PRJNA303987). The samples are colored based on IVM dosage and time period. (D) Principal Components Analysis for Study II (PRJNA303986). The samples are colored based on IVM dosage and time period.

Fig 5. Expression of Candidate Drug Target Genes.

Differential expression of potential drug target genes. The genes and clusters were pulled directly from Fig 1. The dendrogram at the top of the heatmap was generated using pvclust. The red values are approximately unbiased (au) and the green are the bootstrap support (bp) values, both of which are generated by pvclust. The values shown are for illustrative purposes with the full dendrogram available in S2 Fig. The histogram at the bottom shows the distribution of all the z-score values in the heatmap. The heatmap uses a z-score normalization of log₂(TPM) values for 39 genes of interest, pulled directly from Fig 1. The legend at the top is broken into three sections: project color, if the sample was drug treated, and sample life stage. The left hand legend is broken into two sections: the outer section denotes the WGCNA cluster and the inner section denotes if cluster matches the main or inverse WGCNA cluster expression pattern. On the right there is the gene name followed by the common name found in the source data file. Samples were labeled with first author, title, and bioproject from [–29].