SNP Discovery from Transcriptome of Cashmere Goat Skin

Tissue samples

The goats were housed under the same environment with free access to feed and water at the Aerbasi White Cashmere Goat Breeding Farm (Inner Mongolia, China). The skin samples were collected from the lateral body of four two-years-old goats at four distinct hair cycle stages. March, June, September, and December represent telogen, early anagen, cashmere fast-growing period, and catagen, respectively. A total of sixteen samples were collected and stored in liquid nitrogen.

RNA isolation and detection

Total RNA was extracted using the Trizol method and treated with RQ1DNase (Promega, Madison, WI, USA) to remove DNA. The quality and quantity of the purified RNA were determined by measuring the absorbance at 260 nm/280 nm (A260/A280) using smartspec plus (BioRad, Hercules, CA, USA). RNA integrity was further verified by 1.5% Agarose gel electrophoresis, and then further sequencing was conducted with pooled RNA samples from each stage.

cDNA library construction and high-throughput sequencing

Construction of cDNA library was independently prepared using 10 μg of total RNA extracted from samples at each of the four stages. Polyadenylated mRNAs were purified and concentrated with oligo (dT)-conjugated magnetic beads (Invitrogen, Carlsbad, CA, USA) before being used for directional RNA-seq library preparation. Purified mRNAs were iron fragmented at 95°C followed by end repair and 5′ adaptor ligation. Reverse transcription was performed with Reverse transcriptase primer harboring 3′ adaptor sequence and randomized hexamer. The cDNAs were purified and amplified and PCR products corresponding to 200 to 500 bps were purified, quantified and stored at −80°C until used for sequencing. For Illumina sequencing, the cDNA libraries were prepared following the manufacturer’s instructions and applied to Hiseq 2000 system (Gene Denovo Co, Guangzhou, China) for 100 nt paired-end sequencing by Majorbio Inc (Shanghai, China).

Reads mapping and SNP detection

The genome sequence of Capra hircus breed Yunnan black goat and the gene annotation were downloaded from Goat Genome Database (GGD) (http://caprinae.kiz.ac.cn/). Sequence reads without Ns passed the quality filter and the adapter sequences from read ends were removed. The remaining reads were aligned to the genome using TopHat2 (http://ccb.jhu.edu/software/tophat/index.shtml) allowing at most one mismatch. SNPs were detected using mpileup in SAMtools (http://samtools.sourceforge.net/, free software) (tools for alignments in the SAM format), SAMtools can identify single base substitutions as well as small insertions and deletions; however, only SNPs were considered in the current analysis. A detected SNP was kept only if the read depth was more than 10 and the quality score was over 20. Gene expression was calculated measuring the read per kilobases per million mapped reads (RPKM) values, which were normalized by total reads.

Collecting SNP statistics

SNP frequency was calculated by dividing 1 kb of reference by number of SNPs (Num_per_killobase). SNP information such as SNP distribution among genes, SNP distribution across genomic regions, SNP depth and Minor allele frequency (MAF) were extracted from the result file of SAMtools. Functional effects of SNPs (e.g., intergenic, 5′UTR, 3′UTR, CODING) were carried out using SnpEff (Cingolani et al., 2012).

Functional annotation of detected SNPs

Genes containing SNPs were annotated by Blast2GO (https://www.blast2go.com/) with Gene Ontology (GO) terms under the categories of biological processes, molecular functions, and cellular components at level 2 (Götz et al., 2008). Annotated information was imported into BGI WEGO program (http://wego.genomics.org.cn/cgi-bin/wego/index.pl) in WEGO native format. KEGG pathway analyses were conducted using the online KEGG Automatic Annotation Server (http://www.genome.jp/kegg/).

SNP validation

To validate the accuracy of detected SNPs, we randomly selected 30 SNPs for SNP validation by Sanger sequencing using DNA as templates. Primers were designed to amplify the flanking sequence of selected SNPs using Primer 5 (Premier biosoft, Palo Alto, CA, USA) with a fragment length of 200 to 500 bp. Primers were synthesized by Sango Biotech (Shanghai, China) and are listed in Table S2. Thirty individuals were used for construction of the DNA pool, which were sequenced by Taihegene Biotechnology Co., Ltd. (Beijing, China).

Transcriptome sequencing

RNA sequencing obtained a total of 19 G data, each stage average 4.7 G data. The sequencing results demonstrated that 77% to 82% reads with high quality of 97 nucleotides (nt) were mapped to the reference goat genome and 91% to 93% unique-mapped-reads were obtained from total reads at the four stages. QC results in each sample are listed in Table 1 and tell us that the probability that the corresponding base call is incorrect is less than 0.001. Quality score across all bases of March skin sample is shown in Figure 1.

Table 1

Quality control results and the high-quality clean reads of each sample

Figure 1

Quality score across all bases of March skin sample.

SNP statistics

Using the reference goat genome, we identified a total of 56,231 putative SNPs for the four stages. The detailed SNP information is provided in Table S1. Among the 56,231 putative SNPs, 63.6% were found in coding sequence (CDS) regions, 28.2% in intergenic regions, and 8.2% in introns (Figure 2). As shown in Table 2, of the 56,231 putative SNPs, 35,763 SNPs (63.6%) were identified in exons, which were highly represented, while 4,611 SNPs (8.2%) were identified in the introns, which were lowly represented.

Figure 2

Single nucleotide polymorphism occurrence over coding sequence, intron, and intergenic genomic regions.

Table 2

Classification of putative SNPs

A total of 20201 expressed genes (RPKM>0) were detected and 18,433 were annotated among them. We also analyzed the occurrence statistics of SNP in annotated genes. A total of 10,057 genes containing SNPs were identified. On the average, there are 5.59 SNPs per gene. Among these genes, genes with 1 SNP were more common and those with less than 10 SNPs accounted for 93.6% of all SNP containing and expressed genes. The number of genes with 11 to 15 SNPs per gene is 445 and there were 20 genes harboring more than 30 SNPs per gene. The histogram of the number of SNPs per gene is shown in Figure 3. In order to investigate the mutation rate among genes, the SNP frequency within genes was calculated. The top twenty annotated genes with the highest SNP frequency are listed in Table 3.

Figure 3

Single nucleotide polymorphism (SNP) distribution among genes. The horizontal axis represents number of SNPs per gene.

Table 3

The top twenty annotated genes with the highest SNP frequency

The read depth for each SNP was calculated based on the sequence data. The majority of SNPs (57.2%) have a read depth between 10 and 50, while 20% SNPs have read depths ranging from 51 to 100 (Figure 4).

Figure 4

Statistics of single nucleotide polymorphism (SNP) read depth in Cashmere goat skin transcriptome. The horizontal axis represents the read depth of SNPs. The vertical axis represents the number of SNPs with the corresponding read depth.

MAF was the frequency of the less common allele of a polymorphic locus in a given population. The average MAF was 18% in putative SNPs identified from the skin transcriptomes of Cashmere goats (Figure 5). As only SNPs with MAF more than 0.05 are regarded as true SNPs. SNPs with more than 0.05 MAF accounts for the 79.9% and the average MAF of the considered true SNPs was 22%.

Figure 5

Statistics of minor allele frequency (MAF) of discovered single nucleotide polymorphisms (SNPs) in Cashmere goat skin transcriptome. The horizontal axis represents the SNP MAF in percentage, while the vertical axis represents the number of SNPs with given MAF. The average MAF is 0.18.

Functional annotation of SNP

After the GO annotation, 9,392 expressed genes containing SNPs were assigned with one or more GO terms. The plotted GO annotations of these annotated genes are shown in Figure 6. For the cellular component, the major categories were cell and cell part. For molecular function, genes involved in the binding and catalytic activity are highly represented. Under the GO category of biological process, cellular process was the most represented GO term, followed by metabolic process and biological regulation. The GO categories of expressed genes containing SNPs were in the same proportion to the GO categories of all expressed gene.

Figure 6

Gene ontology of all annotated genes in Cashmere goat skin transcriptome and the expressed single nucleotide polymorphism-containing genes.

KEGG pathway analysis categorizes gene functions by biochemical pathways. A total of 2381 expressed genes were assigned with one or more KEGG annotation and were assigned to 304 KEGG pathways. Of these KEGG pathways, metabolic pathways including carbon fixation, nitrogen metabolism, lipid metabolism contains the most annotated genes. Forty genes were classified into the Wnt signaling pathway. TGF-beta signaling pathway contains 27 annotated genes and 36 genes were classified into the NF-kB signaling pathway. The top 10 KEGG pathway classification of the top 200 highly expressed gene contained SNP (Figure 7).

Figure 7

The top 10 KEGG pathway classification of the top 200 highly expressed gene contained single nucleotide polymorphism.

SNP validation

Of the 30 SNPs we selected to evaluate for the validation rate, 22 could amplify target sequences. Within these amplified sequences, 14 SNPs were validated (Table 4). The estimated predicting accuracy reached 64%.

Table 4

The result of SNP validation