Expression profiles of five tissues
Brain, liver, lung, muscle and placenta samples were analysed from of the 4 genetic combinations with 3 males and 3 female samples in each group for a total of 120 samples. Between 60-100M 100bp PE reads, or 90-130M 75bp PE reads per sample passed quality control. Reads were aligned to the extended Brahman reference genome (UOA_brahman_1 plus non-PAR Y chromosome from UOA_angus_1) using hisat2 with default settings, giving an average mapping rate of 89%. The total number of expressed genes among samples ranged from 16,368 to 17,013 and showed no substantial variation between tissues. There was a high correlation coefficient between pure breed Bi and Bt expression of the same genes in each tissue (Supplementary figure 1a-e). There were 14,143 genes expressed in all tissues (Supplementary figure 1f) with 5 genes consistently highly expressed in all five tissues: Insulin-Like Growth Factor 2 (IGF2), Eukaryotic Translation Elongation Factor 1 Alpha 1 (EEF1A1), Collagen Type III Alpha 1 Chain (COL3A1), Actin Beta (ACTB) and the paternally expressed gene 3 (PEG3).
Each tissue formed a tight cluster by multi-dimensional scaling analysis of the 5 tissues which are distinct and well separated from the others, irrespective of the genetic group (Figure 1a). A multi-factor model was used to account for and remove tissue effects, after which a PCA separated the samples by genetic groups in the first principle component (x-axis) and by sex in the second principle component (y-axis) (Figure 1b). The pattern of expression for each tissue from each genetic type showed the same pattern within sex, with the 2 purebred groups separated for all tissues, and the reciprocal crosses less well separated (Supplementary Figure 2a-e). The 20 most highly expressed genes in each tissue are reported in Supplementary Table 1.
Differential gene expression between purebred groups
There were 1,085, 1,495, 1,935, 2,515 and 2,645 genes for which the normalized average number of mapped reads (CPM) differed significantly between purebred Bt and Bi brain, placenta, lung, liver and muscle, respectively. We designated these as differentially expressed genes (DEGs). A similar number of genes showed higher or lower average transcript abundance in Bt as Bi for all five tissues. Muscle had the largest number of DEGs among the tissues studied, but about 84% of these showed a fold change (FC) < 2, while in other tissues ~62-72% showed a FC < 2. The most significantly enriched gene ontology (GO) biological process and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways in muscle included collagen metabolic process (GO:0032963); collagen fibril organization (GO:0030199); amino sugar and nucleotide sugar metabolism (bta00520) and glycine, serine and threonine metabolism (bta00260). Genes in all four of these pathways had higher expression in Bt than in Bi.
Among the DEGs, ~10% in each tissue were lncRNAs. About 92% of DE lncRNAs had the opposite transcriptional direction to differentially expressed genes located within 100 kb.
DEGs common to all five tissues in pure-bred groups
There were 110 DEGs between Bi and Bt in common for all five tissues, comprising 50 annotated protein-coding genes, 42 genes lacking annotation in the reference genome and 18 lncRNAs (Figure 2a). Alignment of the unannotated protein-coding genes to known genes in other cattle and ruminant reference genomes facilitated the annotation of 37 of the unnamed DEGs, based on >90% sequence identity. Of the 87 genes for which annotation was obtained and that were DE in all five tissues between the purebred animals, 84 had consistent relative abundance between subspecies Bt and Bi with respect to genotype in all tissues. The 3 exceptions were Aldehyde Oxidase 1 (AOX1), Choline Dehydrogenase (CHDH), Syntaxin 11 (STX11), whose expression was in a different direction (Bt vs Bi) in the liver compared with the other 4 tissues.
GO pathway analysis of the set of 87 annotated genes showed that they were significantly enriched in 10 GO terms, including oxidation-reduction process (GO:0055114), intracellular protein transport (GO:0006886), glycogen catabolic process (GO:0005980), positive regulation of protein autophosphorylation (GO:0031954) (Figure 2b).
Tissue-specific genes between purebred groups
Genes that were DE between purebred Bt and Bi in only one of the five tissues examined were considered tissue specific DEGs. Using a FDR cut-off of <0.05, liver and muscle had a substantially higher number of unique DEGs compared with brain, placenta, or lung (see Figure 2a). After filtering for an absolute FC ≥ 2, two thirds of the tissue-specific DEGs were removed, leaving 187, 328, 289, 388 and 191 tissue-specific DEGs in each of the tissues respectively, and liver and muscle no longer had substantially higher numbers of unique DEGs. GO biological process pathway enrichment analysis for these filtered tissue-specific DEGs identified 54 GO terms (Supplementary Table 2). The liver-specific DEGs were enriched for 6 GO terms including ion binding and primary metabolic processes. Muscle was enriched for 9 GO terms including the collagen fibril organization pathway. Brain was also enriched for 9 GO terms that included pathways involved with detection of stimulus and nervous system processes. Lung was enriched for 10 GO terms, most of which were related to fundamental biological processes, including regulation of molecular function and cellular response to endogenous stimulus. Placenta was enriched for 20 GO terms which were linked to proton-transporting V-type ATPase and V1 domain small molecule metabolic process.
Differential gene expression between crossbred groups
Comparison of transcript abundance between the reciprocal cross-bred groups (Bt x Bi and Bi x Bt) did not reveal a substantial number of DEGs (<20/tissue at FDR <0.05), except for liver which had 2,473 DEGs. However, only 143 (5.8%) of the liver DEGs had a fold change greater than 2. We performed GO biological process pathway enrichment analysis and KEGG pathway enrichment analysis for the protein coding DE genes with >2-fold change. The GO analysis showed that DEGs were significantly enriched in 6 GO terms, including: macromolecule metabolic process (GO:0043170), primary metabolic process (GO:0044238), cellular metabolic process (GO:0044237), metabolic process (GO:0008152), nitrogen compound metabolic process (GO:0006807) and organic substance metabolic process (GO:0071704) which are all involved in metabolic processes. The only significantly enriched KEGG pathway was metabolic pathways (path: bta01100).
Pairwise comparisons of the DEGs in liver for the 4 genetic groups were performed to explore relationships in expression patterns between pure bred and crossbred concepti. The sire dominated the liver expression pattern in Bt-sired crossbred (Bt x Bi) liver which had 1,276 DEGs when compared to purebred Bi liver, versus 219 DEG when compared with purebred Bt liver. However, the dam breed appears to dominate expression pattern in Bi-sired crossbreds, with 317 DEGs in the Bi x Bt crossbred compared with purebred Bt, but 150 DEGs when compared with purebred Bi liver transcripts.
Expression pattern of DEGs from the purebred groups in comparison with crossbred groups
The expression pattern of the 6,456 DEGs between tissues of purebred animals was examined in the reciprocal crossbred groups. Of these DEGs 5,784 (~90%) showed an additive expression pattern where both paternal and maternal genomes contributed to the gene expression levels in the crossbred groups (Figure 3a), as suggested by the transcript abundance falling approximately midway between that of the two purebred classes. However, transcript abundance of some DEGs (672) was more consistent with parent-of-origin driven expression (Figure 3b-i). Different types of such effects were observed, predominantly maternal/paternal dominance and Bt or Bi allele derived dominance. The abundance of DEGs between crossbred groups fell into three general categories: co-dominant, dominant and recessive expression patterns, with dominance in some cases driven by either the male or the female (Figure 3) the number of genes falling into each category are given in Table 1.
GO analysis of the DEGs that overlapped between tissues showed that they were significantly enriched in 19 GO terms including positive regulation of cellular metabolic process, positive regulation of nitrogen compound metabolic process and membrane-enclosed lumen. The transcript levels of the DEGs involved in these significantly enriched pathways had exclusively higher expression in the purebred Bt compared with the purebred Bi.
Boxplots illustrating the different expression patterns observed among the 4 genetics groups: Bt x Bt, Bi x Bt, Bt x Bi and Bi x Bi (sire breed given first). Y-axis is expression level (counts per million) in a log2 scale. a) Taurus driven additive expression, irrespective of parent. b) Maternal genome driven indicine dominance. c) Maternal genome driven taurine dominance. d) Paternal genome driven indicine dominance. e) Paternal genome driven taurine dominance. f) Taurine dominant – activation. g) Taurine dominant - inhibition. h) Indicine dominant - activation. i) Indicine dominant – inhibition. j) complex inheritance.
Table 1: Number of genes showing a parent of origin effect on expression patterns in five tissues.
|
Brain
|
Liver
|
Lung
|
Muscle
|
Placenta
|
Maternal genome driven - Taurine
|
0
|
12
|
2
|
16
|
0
|
Maternal genome driven - Indicine
|
0
|
11
|
4
|
11
|
5
|
Paternal genome driven – Taurine
|
6
|
24
|
9
|
16
|
7
|
Paternal genome driven - Indicine
|
5
|
36
|
13
|
13
|
12
|
Taurine dominant – activation
|
27
|
65
|
40
|
33
|
25
|
Taurine dominant – inhibition
|
7
|
30
|
24
|
10
|
12
|
Indicine dominant – activation
|
13
|
38
|
17
|
15
|
10
|
Indicine dominant - inhibition
|
2
|
22
|
3
|
15
|
5
|