Abstract

In the current study, allele specific expression analysis was performed in two subspecies cows (Bos taurus and Bos indicus) at SNP and gene levels. RNA-Seq data of 21,078,477 and 20940063 paired end reads from pooling of whole blood samples (Leukocyte) from 40 US Holstein (Bos Taurus) and 45 Cholistani cows (Bos indicus) obtained from SRA database in NCBI. Quality control and trimming of row RNA-Seq data were processed by FASTQC and Trimmomatic softwares. The transcriptome was assembled by TopHat2 software in two cow’s population by aligning and mapping the RNA-Seq reads on bovine reference genome. The SNPs were discovered by Samtools software and ASE analysis was performed by Chi-square test. Results showed that 50183 and 137954 SNPs were discovered on the assembled transcriptome of Holstein and Cholistani cow samples, respectively, and 15308 SNPs were common in both breeds. 10158 SNPs from 50183 (20%) in Holstein and 31523 SNPs from 137954 (23%) in Cholistani cows were identified as ASE-SNPs. Reference allele and alternative allele count in Holstein and Cholistani cows were 3041 and 7155, respectively. Among 131 discovered SNPs in 41 genes with different expression in Holstein and Cholistani cows, 31 ASE-SNPs (5 in Holstein; 26 in Cholistani cows) were discovered.

Key words
SNP discovery; transcriptome; Cholistani cows; Holstein cows

INTRODUCTION

Allele Specific Expression (ASE) is the phenomena that two alleles of the same loci are expressed differently (Gu & Wang 2015GU F & WANG X. 2015. Analysis of allele specific expression-A survey. Tsinghua Sci Thechnol 9: 513-529.), and its a powerful method that measures the expression of each allele through SNP in RNA samples. ASE is an important aspect of gene regulation and one of the important genetic factors that lead to phenotypic variation can be used to identify the variance of gene regulation factors (Gaur et al. 2013GAUR U, LI K, MEI S & LIU G. 2013. Research progress in allele-specific expression and its regulatory mechanisms. J Applied Genet 54: 271-283., Mayba et al. 2014MAYBA O, GILBERT H, LIU J, HAVERTY P, JHUNJHNWALA S, JIANG Z, WATANABE C & ZHANG Z. 2014. MBASED: Allele-specific expression detection in cancer tissues and cell lines. Anim Genet 15: 405-426.). Although the majority of genes are expressed equally from both alleles, some genes are differentially expressed. Besides the gene expression differences between species, the inter individual differences in gene expression are often highly heritable and can be highly context-specific (Wayne et al. 2004WAYNE ML, PAN YJ, NUZHDIN SV & MCINTYRE LM. 2004. Additivity and trans-acting effects on gene expression in male Drosophila simulans. Genet 168: 1413-1420., Gibson & Weir 2005GIBSON G & WEIR B. 2005. The quantitative genetics of transcription. Trends Genet 21: 616-623., Hughes et al. 2006HUGHES KA, AYROLES JF & REEDY MM. 2006. Segregating variation in the transcriptome: cis regulation and additivity of effects. Genet 173: 1347-1355., Lemos et al. 2008LEMOS B, ARARIPE LO, FONTANILLAS P & HARTL DL. 2008. Dominance and the evolutionary accumulation of cis- and trans-effects on gene expression. Proc Natl Acad Sci USA, p. 14471-14476., Ayroles et al. 2009AYROLES JF, HUGHES KA & ROWE KC. 2009. A genome wide assessment of inbreeding depression: gene number, function, and mode of action. Conserv Biol 23: 920-930., McDaniell et al. 2010MCDANIELL R, LEE BK & SONG L. 2010. Heritable individual specific and allele-specific chromatin signatures in humans. Sci 328: 235-239.). ASE may accumulate with genetic divergence and possibly with adaptation to different environments and are responsive to dynamic developmental processes (Von Korff et al. 2009VON KORFF M, RADOVIC S & CHOUMANE W. 2009. Asymmetric allele specific expression in relation to developmental variation and drought stress in barley hybrids Plant J 59: 14-26.). ASE assays can be used to identify cis, trans and cis-by-trans regulatory variation (Main et al. 2009MAIN BJ, BICKEL RD, MCINTYRE LM, GRAZE RM, CALABRESE PP & NUZHDIN SV. 2009. Allele-specific expression assays using Solexa. BMC Genom 10: 422-430.).

RNA sequencing (RNA-Seq) is a powerful new method for mapping and quantifying transcriptomes developed to analyze global gene expression. In other words, RNA-Seq is a next generating sequencing based technology for studying of whole transcriptome and gene expression. This technique provides insights at multiple levels into the transcription of the genome as it yields sequence, splicing and expression-level information, so provides a far more precise measurement of levels of transcripts and their isoforms than other methods (Wang et al. 2009WANG Z, GERSTEING M & SNYDER M. 2009. RNA-seq: a revolutionary tool for transcriptomics. Nature Rev Genet 10: 57-63.). It simultaneously enables study of transcriptomics sequences and very accurate quantitative gene expression (digital expression). Hence, these data are very suitable for high-throughput study of expression level of all transcribed genes and their SNPs. Recently, RNA-Seq has also been used as an efficient and cost-effective method to systematically identify SNPs in transcribed regions in different species (Cloonan et al. 2008CLOONAN N, FORREST A, KOLLE G, GARDINER B & FAULKNER G. 2008. Stem cell transcriptome profiling via massive-scale mRNA sequencing. Nat Methods 5: 613-619., Morin et al. 2008MORIN RO, CONNOR M, GRIFFITH M, KUCHENBAUER F & DELANEY A. 2008. Application of massively parallel sequencing to micro RNA profiling and discovery in human embryonic stem cells. Genome Res 18: 610-621., Chepelev et al. 2009CHEPELEV I, WEI G, TANG Q & ZHAO K. 2009. Detection of single nucleotide variations in expressed exons of the human genome using RNA-Seq. Nucl Acids Res 37: 106-113., Cirulli et al. 2010CIRULLI ET, SINGH A, SHIANNA KV, GE D & SMITH JP. 2010. Screening the human exome: a comparison of whole genome and whole transcriptome sequencing. Genome Biol 11: 57-64.). Transcription is the first step in translation of genome to function underlying genetic codes. Therefore, transcriptase might fill the gap between genotype and phenotype and help understanding the mechanisms from sequence to function (Wang et al. 2009WANG Z, GERSTEING M & SNYDER M. 2009. RNA-seq: a revolutionary tool for transcriptomics. Nature Rev Genet 10: 57-63.).

Previous studies discovered SNPs in bovine milk transcriptome using RNA-Seq (Canovas et al. 2010CANOVAS A, RINCON G, ISLAS-TREJO A, WICKRAMASINGHE S & MEDRANO JF. 2010. SNP discovery in the bovine milk transcriptome using RNA-Seq technology. Mammalian Genom 21: 592-598., Wickramasinghe et al. 2012WICKRAMASINGHE S, RINCON G, ISLAS-TREJO A & MEDRANO JF. 2012. Transcriptional profiling of bovine milk using RNA sequencing. BMC Gen 13(1): 45., Banabazi et al. 2016BANABAZI M, NEJATI-JAVAREMI A, IMUMORIN M, GHADERI-ZEFREHI S & MIRAEI-ASHTIANI SR. 2016. Single nucleotide polymorphisms (SNP) on transcriptome of Holstein cows shared with illumina bovine SNP arrays. Online J Vet Res 20(3): 177-182., Pareek et al. 2016PAREEK CS ET AL. 2016. Single Nucleotide Polymorphism discovery in bovine pituitary gland using RNA-Seq technology. PLoS ONE 11: e0161370.). It has been detected 19,175 genes expressed in milk samples corresponding to approximately 70% of the total number of analyzed genes. The SNP detection analysis revealed 100,734 SNPs in Holstein samples, and a large number of those corresponded to differences between the Holstein breed and the Hereford bovine genome (Canovas et al. 2010CANOVAS A, RINCON G, ISLAS-TREJO A, WICKRAMASINGHE S & MEDRANO JF. 2010. SNP discovery in the bovine milk transcriptome using RNA-Seq technology. Mammalian Genom 21: 592-598.).

Chitwood et al. (2013)CHITWOOD JL, RINCON G, KASIER GG, MEDRANO JF & ROSS PJ. 2013. RNA-seq analysis of single bovine blastocysts. BMC Genom 14: 350-365. were analyzed transcriptomics data to identify SNP in individual blastocyst expressed genes, and individual SNP were examined to characterize allele specific expression. Expressed biallelic SNP variants with allelic imbalances were observed in 473 SNP, where one allele represented between 65-95% of a variant’s transcripts.

In recent years, single nucleotide polymorphisms (SNP) have been the most important and efficient tool in animal breeding. About 40% of the SNPs in the genes cause a change in an amino acid. SNPs are either transition or transversion. Transitions are interchanges of two-ring purines (A↔G) or one-ring pyrimidines (T↔C), while transversions are interchanges of purine to pyrimidine and viceversa (G↔C، G↔T، A↔C ،A↔T). Arefnezhad et al. (2015)AREFNEZHAD B, KOHRAM H, MORADI SHAHRE BABAK M, SHAKERI M, DONG Y, ZHANG X, WANG W & HOSSEINI SALEKDEH GH. 2015. Genetic Variant Detection of Caspian Horse Using High-throughput Sequencing Technology (in Persian). J Agric Biotech 4(6): 101-116. reported that transition and transversion nucleotide replacement were 1155417 and 512986 in Caspian horse, respectively, and replacement ratio of transition to transversion (Ts/Tv) for SNPs was 2.25.

The importance of understanding transcriptomic variation is obvious as the role of gene expression in shaping phenotypes is well documented. In particular, the transcriptomic variation among cattle breeds may provide mechanistic knowledge on their differentiation on phenotypes including appearance, physiological, behavioral, and production traits. There is accumulating evidence that variation in gene expression, presumably controlled by genomic variations within regulatory elements, contributes to phenotypic variation (Passador-Gurgel et al. 2007PASSADOR-GURGEL G, HSIEH WP, HUNT P & DEIGHTON N. 2007. Quantitative trait transcripts for nicotine resistance in Drosophila melanogaster. Nat Genet 39: 264-268.). There are substantial phenotypic difference between Holstein and Cholistani cattle. In particular, they differ remarkably in their resistance to thermal stress, parasites, and diseases (Huang et al. 2012HUANG W, NADEEM A, ZHANG B, BABAR M, SOLLER M & KHATIB H. 2012. Characterization and comparison of the leukocyte transcriptomes of three cattle breeds. PLoS ONE 7(1): e30244.).

In the current study, SNP discovery and Allele Specific Expression analysis were performed in two subspecies cows (Bos taurus and Bos indicus) at SNP level and gene level. We used mRNA-Seq to characterize and compare the Leukocyte transcriptomes of US Holstein and Cholistani cows. These variations may provide partial explanations for differential phenotypes between cattle breeds, particularly between Bos taurus and Bos indicus cattle.

MATERIALS AND METHODS

RNA-Seq data of 21,078,477 and 20940063 paired end reads with 75 bp length resulted from pooling of whole blood samples (Leukocyte) of 40 Holstein cows at the University of Wisconsin, Dairy Cattle Center, USA, and 45 Cholistani cows at Gujait Peer Farm, Bahawalpur, Punjab, Pakistan, respectively, (Huang et al. 2012HUANG W, NADEEM A, ZHANG B, BABAR M, SOLLER M & KHATIB H. 2012. Characterization and comparison of the leukocyte transcriptomes of three cattle breeds. PLoS ONE 7(1): e30244.) obtained from SRA database in NCBI for Holstein cows (http://www.ncbi.nlm.nih.gov/sra/SRX317197) and Cholistani cows http://www.ncbi.nlm.nih.gov/sra/SRS454433). Animal care procedures and all analysis were approved by Ethic Committee (Razi University, Kermanshah, Iran).

mRNA sequencing was run on Illumina Genome Analyzer IIx (Illumina Inc., San Diego, CA). Data were converted from Sra format to Fastq format by fastq-dump command from Ubuntu linux version of Sratoolkit 2.5.4-1. Data quality control was checked by FastQC (v0.11.3) likewise trimmed for linked adaptors and bad quality reads by Trimmomatic 0.33 (Bolger et al. 2014BOLGER AM, LOHSE M & USADEL B. 2014. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30(15): 2114-2120.).Adaptors were considered according to sequencing instrument as default (TruSeq2-PE.fa) and the minimum read length was set at 50 bp. Trimmed reads were aligned on UMD3.1 reference genome (release 81) based on annotation data by Tophat2 (Kim et al. 2013KIM D, PERTEA G, TRAPNELL C, PIMENTEL H, KELLEY R & SALZBERG S. 2013. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol 14(4): R36.), which applies Bowtie2 (Langmead & Salzberg 2012LANGMEAD B & SALZBERG SL. 2012. Fast gapped-read alignment with Bowtie 2. Nat Meth 9(4): 357-359.) as the aligner. The transcriptome was assembled by TopHat2 software in two cow’s population by aligning and mapping the RNA-Seq reads on bovine reference genome. The SNPs were discovered by Samtools software (v. 0.1.19) and ASE analysis was performed by Chi-square test (P.value < 0.01).

RESULTS

Quality control and preparation of RNA-Seq data

After data editing, the removed and low quality reads in both breeds were almost equal and relatively low. For example, amongst the 20940063 initial reads in Cholistani cows, 19379487 reads had high quality and 1560576 reads had low quality, therefore, 5.7% reads were removed from the analysis.

The length of whole transcriptome assembled, for example 52798651 bases in Holstein, indicates around 2% of the whole genome (around 2.6 Mbp) expressed as mRNA. In Cholistani cows, read mapping rate for forward and reverse reads were 81.3 and 79.9%, respectively, and multiple alignments rate was about 9.4%. Overall read mapping was 80.6% and concordant pair alignment was 70.1%. In Holstein cows, read mapping rate for forward and reverse reads were 66.3 and 55.4%, respectively, and multiple alignments rate was about 7.2%. Overall read mapping was 60.8% and concordant pair alignment was 51.3%.

SNP and ASE-SNP discovery

After quality control and filtering, 50183 and 137954 SNPs were discovered on the assembled transcriptome of US Holstein and Cholistani cow samples, respectively, and 15308 SNPs were common in both breeds. The number of discovered SNPs in Cholistani cows (Bos Indicus) was approximately three times higher than Holstein (Bos Taraus) cows (Table I).

Based on the results of Chi-square (χ2) test on 3041 and 7155 loci in the Holstein and Cholistani cows, respectively, it was found that number of reference and alternate alleles were equal.

Totally, in Holstein cows 10158 from 50183 SNPs (20%) were identified as ASE-SNPs. From 10158 loci, number of imbalance alternate and reference alleles were 5006 (49%) and 5152 (51%), respectively. There is generally some bias toward reference allele. This indicates the reference genome has been applied well for mapping RNA reads on both subspecies.

In Cholistani cows, 31523 from 137954 SNPs (23%) were identified as ASE-SNPs. Among 31523 loci, number of imbalance alternate and reference alleles were 21153 (67%) and 10370 (33%), respectively.

SNP and ASE-SNP types on SNP level and gene level

In the present study, 12 SNP types were identified (4 transition and 8 transversion) and the most commonly SNPs were transition SNPs, including 69.6% in Holstein cows and 70.6% in Cholistani cows (Table 2). Replacement ratio of transition to transversion (Ts/Tv) for SNPs was 2.3 and 2.4 in Holstein and Cholistani cows, respectively. The results obtained by Arefnezhad et al. (2015)AREFNEZHAD B, KOHRAM H, MORADI SHAHRE BABAK M, SHAKERI M, DONG Y, ZHANG X, WANG W & HOSSEINI SALEKDEH GH. 2015. Genetic Variant Detection of Caspian Horse Using High-throughput Sequencing Technology (in Persian). J Agric Biotech 4(6): 101-116. confirmed this concept.

In ASE-SNPs, the percentage of transition increased from 69.6% to 71% and 70.6% to 73% in Holstein and Cholistani cows, respectively. Replacement ratio of transition to transversion (Ts/Tv) for ASE-SNPs increase from 2.3 to 2.4 and 2.4 to 2.7 in Holstein and Cholistani cows, respectively (Table II).

In transcriptome of US Holstein and Pakistanian Cholistani cows’ population, 24616 genes have been discovered which 41 genes identified with different expression (Salimpour et al. 2016SALIMPOUR M. 2016. Differential gene expression analysis between the Holstein and cholistani (a Pakistani breed) population using RNA sequencing (RNA-seq) (in Persian). Ms.C. Thesis, College of Agriculture and Naturan Resources, University of Tehran, Iran.). In the current study, from 24616 discovered SNPs in whole genome of Holstein and Cholistani cows population, 131 SNP were found on mentioned 41 genes at Salimpour et al. (2016) report (21 SNPs in Holstein cows and 110 SNPs in Cholistani cows).

From 131 discovered SNP in 41 genes with different expression in Holstein and Cholistani cows population, 31 SNPs were identified as ASE-SNP (5 ASE-SNPs in Holstein cows and 26 ASE-SNPs in Cholistani cows) as shown in Table III.

DISCUSSION

Based on the results of current study the number of discovered SNPs in Cholistani cows (Bos Indicus) was approximately three times higher than Holstein cows (Bos Taraus). Because, for the alignment of both species; which Holstein is a bos taurus and Cholistani (zebo) is a Bos indicus; used a same reference genome with Herford origin, which is also a Bos taurus cow. In addition, stringent settings of tophat2 program were not used in alignment, as with large number of mismatch between the nucleotides on the transcriptome of Cholistani cows and reference genome, alignment may still be successful. Therefore, in SNP discovery analysis, all these mismatches were considered as SNP. Also, above mentioned settings increase relative alignment and mapping rate. Some additional discovered SNPs on the tanscriptome of Cholistani cow are due to 20% higher alignment and mapping rate in Cholistani compared to Holstein cows (70.1% versus 51.3%). The number of discovered SNPs did not correlate with chromosome length (Table I). So, transcription across the genome does not occur with a homogeneous distribution with the same coverage. In other words, some regions contain more candidate genes or important genes that transcription is more intense and deeper in those regions. So, these regions have a larger share of the assembled transcriptomes. Also, the SNPs in these regions have high frequency and remain after filtration.

By SNP screening process, Allelic specific expression (ASE) was identified in both American Holstein and Pakistani Cholistani cows. Gene’s expression levels in Cholistani and Holstein cows have been shown in Table III. Results showed that there are significant different between these two subspecies (P.value < 0.01). Gene ontology (GO) enrichment and pathway analysis revealed that these genes are involved in 20 pathways. A large number of genes are involved on immune response pathways, the electron transport chain and the pathway of translate. These pathways maybe effect on different levels of heat stress and disease resistance. Results showed that most of the genes in metabolic pathways had high expression in Zebo while these genes had low or no expression in Holstein cows, likewise many of these genes are involved on immune pathways in Cholistani cows. Some factors effect on gene expression difference in mentioned two sub-species including: mutation in genes (as Single Nucleotide Polymorphism), epigenetic effects including allele specific expression in this article, environmental effects and gene expression regulatory effects (gene interactions as gene- network). Banabazi et al. (2016)BANABAZI M, NEJATI-JAVAREMI A, IMUMORIN M, GHADERI-ZEFREHI S & MIRAEI-ASHTIANI SR. 2016. Single nucleotide polymorphisms (SNP) on transcriptome of Holstein cows shared with illumina bovine SNP arrays. Online J Vet Res 20(3): 177-182. were found 53478 and 145443 SNPs across the genome on the transcriptome of Holstein and Cholistani cows respectively; that 178 SNPs (24 SNPs in Holstein cows and 154 SNPs in Cholistani cows) were found in 41 detected gene with different expression in current research.

Based on the results there was no SNP in some genes. Generally, a portion of difference in gene expression is due to SNPs and also it could be caused due to regulation of gene expression under different condition or due to epigenetic effects, such as allelic specific expression.

The expression difference between two alleles in a single-nucleotide position causes phenotype diversity and probably explains the large part of variances between these two bovine subspecies, especially in diversity, susceptibility to disease and parasites, tolerating environmental stress such as biological and non-biological stresses in different environmental conditions. While, differential gene expression analysis or even allelic specific expression in gene level may not be able to explain phenotype diversity.

ACKNOWLEDGMENTS

This research was supported by Razi University and Animal Science Research Institute of IRAN (ASRI).

REFERENCES

Publication Dates

History