Validation of expression profiles of differentially expressed transcripts identified in cattle-derived and buffalo-derived Theileria parva isolates by RNA sequencing

Abstract

The intracellular apicomplexan protozoan parasite Theileria parva is a causative agent of cattle theileriosis which manifests in two disease syndromes, namely East Coast fever (ECF) and Corridor disease. Although ECF was eradicated from South Africa, cattle theileriosis still persists in the form of Corridor disease. Moreover, it is not known if the T. parva parasites present in buffalo in South Africa could cause ECF if they were to become established in cattle. This has made it essential to identify genetic differences that would allow successful discrimination of cattle-derived (causative agents of ECF) and buffalo-derived (causative agents of Corridor disease) T. parva parasites. Consequently, Next Generation Sequencing (NGS) was utilized to analyse T. parva transcriptomes from two isolates representing cattle-derived and buffalo-derived parasites, in order to identify gene expression profiles that may characterize cattle-derived and buffalo-derived T. parva isolates. However, RNA-sequencing (RNA-seq) experiments can be influenced by variability caused by technical effects including multiple template preparation stages, diverse sequencing chemistries and complex data processing of NGS experiments; it is thus crucial that data from these experiments is validated using other technologies. Thus, the aim of this study was to use quantitative real-time polymerase chain reaction (qPCR) for validation of differentially expressed genes (DEGs) identified from the RNA-seq study using NGS. Three groups of genes representing different expression profiles, including: 1. constitutively expressed genes; 2. up- and down-regulated genes and 3. genes exclusively expressed in one isolate or the other, were selected for validation. Prior to validation of expression profiles for the selected set of genes using qPCR, endogenous control genes had to be selected in order to normalize qPCR gene expression data. Since there is no information available on the evaluation of the expression stability of these genes in T. parva isolates, the expression stability of five candidate reference genes, β-actin, glyceraldehyde-3-phosphate dehydrogenase, 28S rRNA, cytochrome b and fructose bisphosphatase aldolase (F6P), was evaluated for identification of reliable reference genes. The outcome of the stability rankings for each gene varied according to the program showing that the criteria for stability ranking differ from program to program. It is for this reason that the RefFinder tool, used in this study, integrates the different programs and gives a recommended comprehensive ranking. Therefore, based on this comprehensive analysis between the two T. parva isolates investigated, 28S rRNA and β-actin genes were selected as most suitable reference genes for this study. Intra- and inter-assay variation analysis of the selected reference genes showed that there was no significant variation in the expression of these genes between the two T. parva isolates with the p values being less than 0.05 and the coefficient of variation percentage being low (<2) for all the genes tested. Thus, we propose that genes coding for 28S rRNA and β-actin proteins be employed as endogenous control genes in studies that involve gene expression analysis of T. parva. Validation of expression profiles from RNA-seq data obtained using NGS was performed using qPCR. In this study, the comparative CT method for qPCR data analysis was employed to analyse the expression profiles of selected genes. The use of this method requires initial validation by ensuring that the target genes have approximately the same amplification efficiency as the endogenous control genes. Therefore, the amplification efficiencies of target genes and endogenous control genes were evaluated by constructing validation plots from standard curves generated from selected constitutively expressed and differentially expressed genes, in comparison to the standard curves of the two endogenous control genes. Initially, cDNA was prepared from total RNA isolated from bovine and buffalo lymphoblastoid cell cultures infected with Theileria parva (Muguga) and Theileria parva (7014), respectively, previously used for RNA-seq by NGS. The quantity of the parasite cDNA from the two isolates was interpolated from the standard curve and standardized to a concentration of 36.03 ng/μl, to eliminate concentration bias in downstream gene expression analysis. This study passed the comparative CT method validation experiment since the absolute slopes of ΔCT vs. Log input cDNA for the selected target genes were all less than 0.1 as required. Twenty DEGs, constituting up- and down-regulated genes and genes exclusively expressed in one isolate or the other, and 10 stably expressed (constitutive) genes were selected for validation of expression profiles from RNA-seq data obtained using NGS. Discrepancies between RNA-seq and qPCR analyses were observed from all three groups of target genes but mostly in the constitutively expressed group of genes; in this group only 40% of the qPCR results corroborated with RNA-seq findings while 60% demonstrated variations in expression with four genes down-regulated and one up-regulated in T. parva 7014 relative to T. parva Muguga. Since most of the disagreements in the two datasets were down-regulated expression, this finding suggests that RNA-seq was more sensitive in detecting low abundant RNA transcripts.