Keywords

Biological process; Bos indicus; Intramuscular fat; Lipid droplet; Variants; SNP

Introduction

Genetic variations in genes involved in the lipids storage in muscle may represent useful markers in marker-assisted selection programs, which may interfere with important phenotypic traits such as intramuscular fat deposition (IMF) in cattle. The main biological processes involved with IMF are adipogenesis (differentiation of adipocytes in preadipocytes) and lipogenesis (fat synthesis), controlled by different genes. Some of these genes are acetyl coenzyme-A carboxylase (ACC), adiponectin (adipoQ) e perilipin 2 (PLIN2) (Hoashi et al., 2008; McManaman et al., 2011), which were associated with the amount of body fat deposited in mammals. PLIN2 is a member of the PAT Family, which is involved in lipid droplet formation in the liver and some peripheral tissues. PLIN2 is also known as adipose differentiation-related protein (ADRP), or adipophilin, whose expression was detected in several cell types including fibroblasts, endothelial and epithelial cells and tissues, such as mammary gland and adrenal cortex, and can be used as a marker of cells lipid accumulation (Kimmel et al., 2010). Specifically, the PLIN2 gene encodes a 50 kD protein that is expressed in most tissues. It is involved in the uptake of long chain fatty acids and the formation and stabilization of lipid droplets (Brasaemle et al., 1997).

Previous studies have shown that the low expression of PLIN2 in mice fed a high-fat diet was associated with lower lipid accumulation in the liver (Imai et al., 2007). In studies with pigs Gandolfi et al. (2011) also reported that PLIN2 expression was associated with the availability of extracellular lipids to skeletal muscle cells. However, there is little information about genetic variation in PLIN2 of Nellore cattle (Bos indicus) and its association w 86 ith important traits such as intramuscular fat deposition.

Nellore breed is the most important breed raised in Brazil and is a Bos indicus species, which is genetically less predisposed to deposit intramuscular fat than Bos taurus species (Lehnert et al., 2007; Wang et al., 2009). The detection of markers that can be used in the future in the breeding program, as well, a better understanding of the genetic architecture of important traits such as IMF, is an important issue for beef quality and production. Thus, the main goal of this study was to identify and characterize SNPs in PLIN2 gene of Nellore cattle through a next-generation sequencing approach.

Materials and Methods

Animals

Six animals were used, three animals with the highest and three with the lowest estimated 99 genomic values (GEBVs) for the intramuscular fat traits measured in the Longissimus dorsi 100 muscle. These animals were selected based on data previously estimated and published by Cesar 101 et al. (2014).

Extraction of genomic DNA and primer design

Genomic DNA was previously extracted from blood samples, quantified and evaluated for quality, as described by Cesar et al. (2014). To design the primers, approximately 3 kb fragments of the PLIN2 gene were selected considering the size of the gene and the distance between the exons to cover the entire length of the gene including up and downstream regions.

To this end, five different primers were designed, with an overlap of at least 100 base pairs (bp) between the drawn primers. The design of the primers was performed based on the data deposited in NCBI from the program Primer 3 (http://109bioinfo.ut.ee/primer3-0.4.0/). The quality of the candidate sequences and the drawn primers were tested using the Net primer program (http://www.premierbiosoft.com). Subsequently, the best combination of forward and reverse primers (Table 1) was chosen. These were then tested using the BLAST tool (Basic Local Alignment Search Tool, http://blast.ncbi.nlm.nih.gov/Blast.cgi) for confirmation of similarity with the Bos indicus species and the desired region.

Table 1. Sequence and size of amplicons of each primer used.

Amplicons design

For the amplification of the fragment, a polymerase chain reaction (PCR) with each pair of primer designed was performed using the Veriti^® thermal cycler (Applied Biosystems, USA). In general PCR amplification was performed under the following conditions: initial denaturation at 98°C for 30 seconds, followed by 30 cycles the in denaturation at 98°C for 10 seconds, annealing at 64°C for 30 seconds, and extension at 72°C for 3 minutes and ending with a final extension of 72°C for 3 minutes. All primer pairs followed the same amplification conditions. For PCR, 20 μl of the reaction containing 0.3 μM of each primer, the enzyme Taq polymerase Emeraldamp Max Hs (Takara Bio, USA) and 40-60 ng of the template DNA were used. Visualization and verification of the amplicon sizes were performed in 1% agarose gel electrophoresis.

Amplicons purification

The amplicons were purified by means of Agencourt AMPure XP magnetic beads (Beckman Coulter, USA), with the recommended volume of Agencourt AMPure XP (beads) for the sample volume, following the recommended by the manufacturer.

New generation sequencing and data analysis

For DNA libraries preparation, Nextera^® XT DNA and Nextera^® XT Index (Illumina, San Diego, USA) were utilized. In the first stage of libraries preparation the amplicons were fragmented by enzymatic reaction (transposomes) and simultaneously t 132 he adapters were added at their ends. Sequencing was performed on MiSeq platform (Illumina, San Diego, USA) using the MiSeq Reagent Kit v2 (500cycle), paired-end format (both directions of the DNA strand) and the size of 250 bases.

The quality of the sequencing data was analyzed and visualized through the FASTQ program version 0.11.5. The Seqyclean version 1.9.11 program was used to remove possible contaminant sequences (reference Univec database), reads with less than 200 bases and quality Phred score below 24. Finally, the mapping was performed using BWA software version 0.7.17 against the bovine reference genome (Bos taurus UMD3.1), formatted using the Picard program version 2.6.0 and the SNP calling was performed by the GATK 3.6 program using the haplotype caller option. The SNPs identified were filtered according to the quality parameters (QUAL> = 30, DP> = 10, QD> 2.0, FS> 60.0, MQ> 40.0, MQRankSum> -12.5 and ReadPosRankSum> -8.0).

Annotation of SNPs and identification of non-tolerable effects

After filtering, SNPs’ functional annotation was performed against the Bos taurus UMD3.1 reference genome from the Ensembl database release 94, using the VEP tool (Variant Effect Predictor, http://www.ensembl.org/info/docs/tools/vep/index.html) (McLaren et al., 2016). The aim of the annotation was to describe SNPs location in the genome (regions of exons, introns, intergenic, 3'UTR, etc.) consequence of mutation, classification based on the amino acid change (synonyms and non-synonyms) and possible effect on protein function. To predict whether the non-synonymous variants affect the protein function the Sorting Intolerant From Tolerant (SIFT) algorithm (Ng & Henikoff 2003) was used, which calculates a scaled probability (SIFT score) of observing a specific amino acid at a position from multiple sequence alignment of homologous proteins. When the score lies at or below the 0.05 threshold the variant is considered not tolerated and when the variant score lies 155 above the 0.05 threshold the variant is considered tolerated (Ng & Henikoff 2003).

The dataset supporting the conclusions of this article is available in the European Variants Archive (EVA) repository (EMBL-EBI), under accession PRJEB32049 and analyses ERZ857797. [https://www.ebi.ac.uk/eva/].

Results and Discussion

PLIN2 amplification and sequencing data

The genomic region that harbor PLIN2 gene (chr8: 25,129,104 – 25,143,157) was amplified by conventional PCR using four different primers from six different animals with extreme values of genomic estimated breeding value (GEBV). Three of the animals presented negative GBV (1:- 0.37, 2:-0.51, and 3:-0.29) and three others positive GEBV (4:0.65, 5:0.58, and 6:0.36) according previous report by Cesar et al., 2014. The purified PCR products from each sample were sequenced by next generation sequencing (NGS) technology. An average 587,401 reads were obtained per sample sequenced on the MiSeq Illumina equipment (Table 1). The number of reads remained after quality control for next analysis is shown in Table 2.

Table 2. Number of sequenced reads, number of reads after filtering by quality and percentage of reads mapped against the Bos taurus UMD3.1 reference genome

Mapping and variant calling

After sequencing quality control the reads were mapped against the Bos taurus UMD3.1 reference genome by the BWA program version 0.7.17. Approximately 99,97% of the reads could be mapped against the reference genome (Table 2). After mapping, the variants call file (single nucleotide polymorphism mutation, SNP) was filtered based on the criteria described in the methods section and 134 SNPs identified in the genomic region of PLIN2 gene (Table 3) remained for further analyzes.

Table 3. Annotation of 134 SNPs identified on PLIN2 genomic region from Nellore steers, 328 according to region and consequence.

Variant effect prediction and annotation

Among the identified SNPs five are novel comparing to Ensembl database release 94. Three of them are in intronic regions and two in the downstream region (5 Kb), which are different of the novel SNPs recently reported by Yue and collaborators (2019). These authors also identified five novel variants using four different Chinese native breeds, which were annotated in exonic region and were associated with growth traits in.

The novel SNPs identified herein were classified as with moderate impact (non-disruptive variant that could change protein effectiveness) on variant consequence by the VEP tool (McLaren et al., 2010), which could promote severe consequences in transcripts and explain phenotypic differences among individuals.

Most of the variants are located in intronic and downstream gene (5 kb) regions (> 88%), eight percent in 3' UTR region of the gene, and seven percent in exonic region of PLIN2 gene (Table 4). More than 71% of the SNPs in exonic regions were classified as non-synonymous, which could play an important functional role in PLIN2 expression and consequently in differences on fat deposition. The non-synonymous variants correspond to the mutations of a nucleotide that can cause alteration of the translated amino acid and synonymous corresponds to the mutations of a nucleotide with no corresponding alteration of the translated amino acid (Hassan et al., 2018).

Table 4. Number, region and estimated consequence of the variants identified on PLIN2 genomic region from Nellore steers by VEP tool.

Herein, the VEP tool was also used to predict the SIFT score of the five non-synonymous variants, which predicted that 80% of the non-synonymous variants are non-tolerated (Table 2). As is well known, PLIN2 is primarily expressed during the early process of adipocyte differentiation, which promotes fat droplet formation in muscle cells. Because of the importance of IMF quantity for meat flavor, tenderness, juiciness, ove 201 rall consumer acceptance, and human health (Killinger et al., 2004), the identification of the genetic architecture and the molecular mechanisms that control IMF deposition has become an important point for meat and human health research. As reported by our group, PLIN2 expression level was associated with trans-eQTL at chr3: 87,253,086 in Nellore cattle, which means that animals with different genotypes showed different expression level among the 193 animals (Cesar et al., 2018). However, in previous genome-wide association study we did not identify a significant association between this trans eQTL and intramuscular fat deposition using 286 animals from the same population (Cesar et al., 2014).

On the other hand, in pig population the authors revealed that PLIN2 gene expression analyses showed a positive correlation with higher intramuscular fat deposition (Davoli et al., 2011). Whereas in human, the PLIN2 overexpression has been shown to decrease the expression level of PPARa target genes as well as the transcriptional activity of mitochondrial genes (Bosma et al., 2012). This effect resulted in higher intramyocellular lipid storage, showing that genetic markers involved in IMF deposition could help to better understand of the biological processes involved in fat deposition, an important trait for beef quality and production (cost of production, feed convertion). As well as biological processes involved in metabolic diseases such as diabetes type 2 and obesity in humans.

The present study allowed the identification of possible causal mutations associated with the intramuscular fat deposition characteristic in the PLIN2 gene. However, the number of animals is a limiting factor to verify the importance of these mutations in the deposition of intramuscular fat among the animals tested. Because we could not achieve association with good accuracy with small number of animals, so the present work had the main goal identify putative mutations in PLIN2 gene, which allow future association studies in the same or different population.

Conclusion

This is the first study performed to identify novel variants of SNP type in PLIN2 gene in Nellore cattle, which is an important gene for IMF, using next generation sequencing. Herein, we could identify five novel variants of SNP type in PLIN2 genomic region, which can potentially change the protein effectiveness and could be used as genomic markers in animal breeding program. However, further association studies are necessary to verify if these novel SNPs are associated with fat deposition traits in Nellore cattle.

Acknowledgments

This work was supported by the EMBRAPA (Macroprograma 1, 01/2005); São Paulo Research Foundation - FAPESP [grant numbers 2014/11871–5, 2014/22884–0, 2012/23638–8, and 2016/26030-1]; National Council for Scientific and Technological Development – CNPq [grants from Luiz L Coutinho and Luciana CA Regitano].

Conflict of Interest

The authors confirmed that there is no conflict of interests.