The transcriptome from asexual to sexual in vitro development of Cystoisospora suis (Apicomplexa: Coccidia).. Transcriptome C.suis

The aim of this project is to perform RNAseq of C. suis harvested at different developmental stages - especially asexual stages (merozoites) and early and late sexual stages (microgametes, macrogametes and early oocyst) - to better understand the transcriptional changes underlying the developmental process. Experimental design, sampling and RNA-seq library preparation For the sampling of sexual stages released from host cells we collected cell culture supernatant every day, from day of cultivation (doc) 6 to day 14. The material was washed twice with phosphate-buffered saline (PBS; Gibco) and pelleted by centrifugation at 600 × g for 10 min. The numbers of merozoites, sexual stages and oocysts were counted in a Neubauer-counting chamber for each given time point. For each day, seven biological replicates were harvested and the mean numbers of each stage per biological replicate were calculated. Pellets from the same wells were pooled to increase the number of parasites per sample and the analysis was performed for three time points: 1. pool of days 6, 7 and 8 (merozoites, type I and II) = Time point 1. 2. pool of days 9, 10 and 11 (merozoites type II and early sexual stages, i.e. gamonts) = Time point 2. 3. pool of days 12, 13 and 14 (mainly sexual stages, gametes, and unsporulated oocysts) = Time point 3. Total RNA was isolated from infected cell cultures using an RNeasy® Mini kit (Qiagen, Hilden, Germany) and treated with RNase-free DNase (Qiagen) according to the manufacturer’s instructions to remove any DNA contamination. Total RNA was quantified using a NanoDrop® 2000 (Thermo Fischer Scientific, Waltham, MA, USA), and samples were sent for library preparation using a reverse stranded protocol with poly-A enrichment. Sequencing libraries were prepared at the Core Facility Genomics, Medical University of Vienna, using the NEBNext Poly(A) mRNA Magnetic Isolation Module® and the NEBNext Ultra® II Directional RNA Library Prep Kit for Illumina according to manufacturer's protocols (New England Biolabs, Ipswich, Massachusetts, USA). Libraries were QC-checked on a Bioanalyzer 2100® (Agilent Technologies, Santa Clara, CA, USA) using a High Sensitivity® DNA kit for correct insert size and quantified using Qubit dsDNA HS® assay (Invitrogen, Waltham, Massachusetts, USA). Pooled libraries were sequenced on a NextSeq500® instrument (Illumina, San Diego, California, USA) in 1 × 75 bp single-end sequencing mode. Approximately 21.5 million reads were generated per sample. RNA-Seq data analysis Sequencing reads were mapped against the concatenated fasta sequences of C. suis (version 48 from ToxoDB)114 and S. scrofa (version 1.11 from Ensembl, GCA_000003025.6)115,116 using STAR (version 2.7.3a with option –outSAMmultNmax 1)117 and the combined annotations of each genome (version 48 for C. suis and 11.1.98 for S. scrofa). Only the reads mapping to the C. suis genome were subsequently used for quantification and further analysis. Quality control was performed with FastQC118 and QualiMap119. RNA degradation was taken into account via the TIN (Transcript Integrity Number) values, which were measured for each gene and library with the RSeQC120. It was used to assess gene body coverage (module geneBody_coverage.py with with option -l 500) and to calculate transcript integrity numbers (TIN scores, module tin.py with option -c 20), TIN is considered an accurate and reliable measurement of RNA integrity at the sample level120,121. Gene expression was quantified with featureCounts (version 1.5.0a)122 with options -s 2 -Q 20 –primary. Identification and analysis of differentially expresed genes All statistical analysis were performed in R (version 4.1)123. Given the repeated measures design of our experiment (briefly, gene expression was measured for seven samples at each of three timepoints) we employed a linear mixed model framework124 to account for the covariance structure in the data. Differential gene expression analysis between the three time points was performed via linear mixed models with the function dream (R package variancePartition, version 1.18.3)124,125, which is a wrapper for the function lmer in package lme4. Replicate ID was fitted as random intercept, and hypothesis testing was carried out for a fixed categorical effect of time with the three time points as factor levels. We further included the median TIN scores, calculated across all genes in each library, as a continuously distributed (nuisance) covariate in our model. We filtered for genes with a minimum count of 30 and four counts per million reads in at least four out of seven replicates of each time point. The remaining counts were quantile normalized before differential gene expression analysis with the function voomWithDreamWeights (R package variancePartition). The p-values were adjusted for multiple testing according to Benjamini and Hochberg's false discovery rate (FDR) correction126. Genes with FDR > 0.05 and absolute log2FC > 1 were considered significantly differentially expressed. Gene ontology enrichment analysis To explore the broader biological context of the identified genes, gene ontology (GO) enrichment analysis was performed via topGO127 with the Fisher’s Exact Test and the GO annotations from ToxoDB (version 50). The “Weighted01” algorithm which accounts for the GO hierarchy was applied.