The combination of GWAS and RNA-seq data analysis provided an opportunity for making a multi-faceted analysis of biological markers for RFI. Using the two methods in tandem, we obtained mutual support for results from multiple analyses.
The priority in our study was the identification of QTLs for RFI and the biological pathway on which they have an impact, thus RNA-seq data was mostly used for verification and interpretation of GWAS results. We have achieved this in four different ways: identification of transcription factors in QTLs with targets in DE genes, finding cis-eQTLs in QTL regions for RFI, finding protein-protein interactions between genes located in QTLs and DE genes, and SNP calling based on RNA-seq data.
GWAS analysis revealed three QTL regions on SSC9, which overlap with a number of QTL found in previous studies according to the Animal QTL Database (34). In total QTL from 26 articles were overlapping with the QTL identified in the present study. In two papers overlapping QTL were associated with such production traits as body weight at day 21 (QTL-1,2,3) (35), lipid (QTL-2) and protein (QTL-1,2,3) accretion rate, and average daily gain (QTL-1,2,3) (36).
Our lead SNP for QTL-2 was located in an intron region of the DNAH11 gene. SNP calling based on RNA-seq data furthermore revealed two deleterious variants in conserved domains of DNAH11. Moreover, SNPs located in QTL were associated with the expression of three genes, however, only one gene, AGMO, shows a p-value lower than 0.05 in the differential expression analysis.
Analysis of DE in the RNA-seq data revealed 187 and 20 genes for RFIcont and RFIdisc, respectively. Expression of five DE genes, ENSSSCG00000038005, MCOLN2, ENSSSCG00000036427, ENSSSCG00000051238, and PARD3B, were associated with SNPs having p-values in GWAS for RFI below 0.05. These polymorphisms may potentially have an impact on RFI through the expression of the mentioned genes, despite of the fact that they were not revealed by GWAS as genome-wide significant SNPs. One of the eSNPs, 2:394039, was located in a QTL for Feed Conversion Ratio according to a previous studies (36).
Interactions between GWAS and DE results were additionally studied using transcription factor–target networks. In this way, one of the transcription factor coding genes located in QTL-1, AHR, can explain the differential expression of 65 DE genes. The AHR gene has been considered a candidate gene for reproductive traits in previous studies (37–39). SNP calling based on RNA-seq data revealed four missense mutations in the protein-coding part of the AHR gene. Thus, we conclude that AHR is a promising candidate gene for an RFI genetic determinant.
Protein-protein interaction revealed a number of interactions between genes located in QTL and DE genes. The most interesting finding was interactions between the QTL gene, AGR2, and four DE genes, HK2, PDIA3, PDIA4, and PDIA6. The last three DE genes were found in one common cluster of genes involved in glycoprotein processing, whereas AGR2 by itself belongs to the disulfide isomerase family and is involved in cysteine-rich glycoprotein mucin processing (40, 41).
In conclusion, we have identified a number of pathways and genes of potential importance for RFI. Of particular interest is for instance the gene encoding the AHR transcription factor, which is located in QTL-1 and targets a large number of DE genes. Another gene of interest is DNAH11, in an intron of which the QTL-2 lead SNP was located. In both genes, missense mutations were found. Further studies in independent cohorts should be performed to confirm the relevance of the candidate genes for feed efficiency in pigs.