After initial variant calling, variants not belonging to exonic regions as per Annovar’s algorithm (5), were filtered out in order to focus on exonic variants. This is since exonic variants are more likely to affect the final protein product, possibly through alterations of gene function, expression, translation or even post-translational modifications. Following this filtration, the study was left with around 76,000 variants, a number that was further brought down to 49,965 after a second filtration based on Wilcoxon test scores as calculated by the BaseQSum calculator functionality of GATK’s variant annotator (9).
These 49,965 exonic variant sites were found to be distributed in a total of 11,568 genes. These variants were categorized according to their effect on translation, i.e. synonymous, non-synonymous, stop-gain, frame-shift or unknown. To statistically analyze the association study, the u-based z approximation derived from the rank-sum test was used instead of the more conventional p-values.
Once more, GATK’s variant annotation functionalities were used to calculate these values which are outputted as the ReadPosRankSum values. Negative values were omitted owing to those representing presence of the alternate alleles near the ends of reads which are often erroneous in these analyses. The remaining positive values were used to plot the frequency of variants per chromosome and used to speculate on possible chromosome-wise association between the variants and possible desired traits (Figure 1).
Genes that prior studies had implicated in the meat tenderization process, encode proteins including calpain-1 (Capn1), calpain-2 (Capn2), calpastatin (Cast), caspase 3 (Casp3), caspase 9 (Casp9), αB-crystallin (Cryab), heat shock protein 27 (Hsp27), heat shock protein 40 (Hsp40) and heat shock protein 70 (Hsp70) (Saccà et al., 2019). Among these, this study found variants in HSP70, CASP3, CAPN1, CAPN2 and CAST encoding genes. Table 1 shows the positions and types of variants found in these genes.
Table 1. Variants found in genes previously implicated in meat tenderization process.
Gene
|
Chromosome
|
Position
|
Base change
|
Amino acid change
|
CAST
|
NC_030814
|
14434086:exon3:c
|
G83A:p
|
R28>K
|
CAST
|
NC_030814
|
14434086:exon18:c
|
C1245G:p
|
V415>V
|
CAST
|
NC_030814
|
14434086:exon28:c
|
G2052A:p
|
E684>E
|
CAST
|
NC_030814
|
14434086:exon25:c
|
C1863T:p
|
P621>P
|
CAST
|
NC_030814
|
14434086:exon14:c
|
G972A:p
|
P324>P
|
HSP70.1
|
NC_030830
|
22440674:exon1:c
|
G309A:p
|
V103>V
|
HSP70.1
|
NC_030830
|
22440674:exon1:c
|
C993T:p
|
I331>I
|
HSP70.1
|
NC_030830
|
22440674:exon1:c
|
C1653T:p
|
S551>S
|
HSP70.1
|
NC_030830
|
22440674:exon1:c
|
C1528T:p
|
L510>L
|
HSP70.1
|
NC_030830
|
22440674:exon1:c
|
C24A:p
|
G8>G
|
CAPN1
|
NC_030836
|
43822608:exon17:c
|
C1803A:p
|
R601>R
|
CAPN2
|
NC_030823
|
25579185:exon12:c
|
G1522A:p
|
D508>N
|
CAPN2
|
NC_030823
|
25579185:exon13:c
|
T1560A:p
|
L520>L
|
CAPN2
|
NC_030823
|
25579185:exon17:c
|
G1764A:p
|
G588>G
|
CAPN2
|
NC_030823
|
25579185:exon18:c
|
C1833T:p
|
Y611>Y
|
CAPN2
|
NC_030823
|
25579185:exon3:c
|
C327T:p
|
A109>A
|
CASP3
|
NC_030834
|
30502447:exon3:c
|
C267T:p
|
N89>N
|
Further analysis looking into the number and nature of non-synonymous variants that occurred in this study’s specimen and the changes they resulted in for the encoded protein revealed the final list to contain a total of 20,879 non-synonymous variants that were distributed in 6850 genes. The results of the gene ontology search also provided some promising findings. In particular, the functional clustering analysis showed that the processes most likely impacted by the identified genes were fatty acid metabolism and fatty acid degradation. Both of these processes may hold significant contributions to the leanness and eventual quality of meat produced from the concerned specimen (10). The third pathway with the most significant scores from the list of genes was the PPAR pathway (peroxisome proliferative-activated receptor) which too is concerned with fatty acid metabolism. Table 2 shows these results alongside the associated statistical scores as calculated by DAVID’s in-built algorithm.
Table 2. Pathways affected by the identified SNPs and their associated statistical scores.
Annotation Cluster 1
|
Enrichment Score: 3.16
|
Count
|
p-Value
|
Benjamini
|
KEGG_PATHWAY
|
Fatty acid degradation
|
17
|
3.40E-06
|
1.80E-04
|
KEGG_PATHWAY
|
Fatty acid metabolism
|
15
|
1.00E-03
|
9.70E-03
|
KEGG_PATHWAY
|
PPAR signalling pathway
|
13
|
9.20E-02
|
3.40E-01
|
Annotation Cluster 2
|
Enrichment Score: 2.8
|
Count
|
p-Value
|
Benjamini
|
KEGG_PATHWAY
|
Protein digestion and absorption
|
36
|
2.90E-08
|
3.80E-06
|
KEGG_PATHWAY
|
ECM-receptor interaction
|
22
|
5.20E-04
|
5.50E-03
|
KEGG_PATHWAY
|
Focal adhesion
|
34
|
2.30E-02
|
1.20E-01
|
KEGG_PATHWAY
|
Amoebiasis
|
20
|
3.10E-02
|
1.50E-01
|
KEGG_PATHWAY
|
PI3K-Akt signaling pathway
|
32
|
9.40E-01
|
9.90E-01
|