Comprehensive analysis of 124 transcriptomes from 31 tissues in developing, juvenile, and adult Japanese Black cattle.

Omic analyses of economically important animals, including Japanese Black cattle, are currently underway worldwide. In particular, tissue and developmental stage-specific transcriptome characterization is essential for understanding the molecular mechanisms underlying the phenotypic expression of genetic disorders and economic traits. Here, we conducted a comprehensive analysis of 124 transcriptomes across 31 major tissues from fetuses, juvenile calves, and adult Japanese Black cattle using short-read sequencing. We found that genes exhibiting high tissue-specific expression tended to increase after 60 days from fertilization and significantly reflected tissue-relevant biology. Based on gene expression variation and inflection points during development, we categorized gene expression patterns as stable, increased, decreased, temporary, or complex in each tissue. We also analysed the expression profiles of causative genes (e.g. SLC12A1, ANXA10, and MYH6) for genetic disorders in cattle, revealing disease-relevant expression patterns. In addition, to directly analyse the structure of full-length transcripts without transcript reconstruction, we performed RNA sequencing analysis of 22 tissues using long-read sequencing and identified 232 novel non-RefSeq isoforms. Collectively, our comprehensive transcriptomic analysis can serve as an important resource for the biological and functional interpretation of gene expression and enable the mechanistic interpretation of genetic disorders and economic traits in Japanese Black cattle.

1. Introduction

Omic analyses of economically important animals are currently underway worldwide. A comprehensive molecular characterization of their genetic components can provide better understanding of their beneficial and adverse traits, including genetic disorders, leading to improved animal husbandry. Japanese Black cattle symbolize these attempts; it is a beef breed, popular for the meat’s high degree of marbling due to intramuscular fat deposition (Fig. 1A). Strict selection for marbled meat under a closed breeding system in Japan has made Japanese Black cattle genetically distinct from other breeds. However, this closed system has resulted in intense inbreeding, fewer founder animals, and a decline in the effective population sizes of Japanese Black cattle., Inbred populations are susceptible to recessively inherited disorders. Notably, artificial insemination (AI) breeding conception rates in Japanese Black cattle have been declining, with the first to third AI conception rates decreasing from 66.4% to 52.6% between 1989 and 2018 in Japan. Previous studies have reported conception failure with embryonic mortality occurring before Day 60 post-fertilization., In addition to embryonic mortality, 4.61% of calves die before within 8 months of birth., Thus, analysing the transcriptome of the early embryonic stage before Day 60 post-fertilization and juvenile calves in Japanese Black cattle may provide valuable insights into the molecular mechanisms underlying these genetic disorders; however, to date, RNA sequencing (RNA-seq) data from major tissues of the developmental stages in cattle are not available.

Figure 1

Japanese Black cattle and Wagyu Genome Database. (A) A Japanese Black sire and cow of the Kagoshima Prefecture. Section between sixth and seventh rib in the Japanese Black cattle. Marbling in meat due to intramuscular fat deposition. (B) RNA-seq data from the Wagyu Genome Database browser (https://wagyu.hgc.jp). bigwig (bw) and TPM derived from short-read and long-read sequencing are shown. Samples and the accompanying data from short-read and long-read sequencing are listed.

Over the last decade, substantial genotyping data using single-nucleotide polymorphism arrays have led to the identification of genomic regions associated with genetic disorders, and economic traits in cattle populations. In addition, next-generation sequencing can rapidly discover candidate causative variants in the genomic regions. However, it is challenging to identify the variants underlying these traits, as genomic knowledge in cattle is limited when compared with that in humans and mice. For this reason, researchers have conducted RNA-seq across a wide range of adult cattle tissues and cells to provide insight into the molecular mechanisms underlying the phenotypic expression of genetic disorders and economic traits.

In this study, we analysed 124 RNA-seq datasets from 31 major tissues at three developmental stages: fetuses at 30, 60, and 90 days post-fertilization, calves, and adult Japanese Black cattle. In addition, to directly analyse the structure of full-length transcripts without transcript reconstruction, we performed RNA-seq analysis of 22 tissues using long-read sequencing. The RNA-seq data are open access and available at the DNA Data Bank of Japan (DDBJ) and Wagyu Genome Database (https://wagyu.hgc.jp). It should also be noted that this is a symbolic case from the viewpoint of data sharing in Japanese Black cattle. Especially for economically important animals, there has been a data-sharing barrier with the concern that it may impose a problem on the animal husbandry business. However, recent international legal coordination related to benefit-sharing has, at least in part, lowered this barrier, which has improved the scientific benefit of collaborative research. Here, we describe the first comprehensive transcriptome analysis of Japanese Black cattle to better understand their unique phenotypic features.

2. Materials and methods

2.1. Ethics

All animal experiments were performed in accordance with the Guidelines for the Care and Use of Animals of the Cattle Breeding Development Institute of Kagoshima Prefecture, University of the Ryukyus, and the Japan Livestock Technology Association. All the cattle used in this study were owned by the Cattle Breeding Development Institute of Kagoshima Prefecture.

2.2. Sample collection and RNA-seq using HiSeq3000 and MinION

We collected 124 samples from seven Japanese Black cattle for RNA-seq analysis, the details of which are summarized in Supplementary Table S1. Samples were collected in 50 mL conical tubes (Iwaki, Shizuoka, JPN, Cat. #2345-050), flash-frozen in liquid nitrogen, and stored at -80°C until RNA extraction. We extracted the total RNA using TRIzol (Thermo Fisher Scientific, Waltham, MA, USA, Cat. #10296010). The extracted RNA was converted to cDNA using the SMART-Seq v4 Ultra Low Input RNA Kit (Takara, Shiga, JPN, Cat. #Z4888N). cDNA was then sequenced using HiSeq3000 (Illumina, San Diego, CA, USA) and MinION (Oxford Nanopore Technologies, Oxford, UK) according to the manufacturer’s instructions.

2.3. Analysis of sequencing data

The sequencing data were mapped to the bovine reference genome (ARS-UCD1.2) using mimimap2 for short reads and STAR for long reads, followed by conversion of data file to a sorted BAM file using SAMtools. We used a gtf file (Bos_taurus.ARS-UCD1.2.105.gtf) as a reference for gene annotation (21,224 genes with coding sequences). We determined gene expression levels [transcripts per million (TPM)] using StringeTie, or RSEM to account for differences in sequence depth and gene length across samples. We classified the 124 samples into seven biological systems: nervous/endocrine, locomotor, cardiopulmonary, digestive, immune, urogenital, and fetus/fetus-support systems (Supplementary Table S1). To evaluate tissue gene expression specificity, we calculated the Z-score by scaling the TPM of each gene to have a mean of zero and variance of one within an individual animal. To evaluate the changes in gene expression during development, log2 [fold change (FC)] was calculated by dividing the maximum TPM by the minimum TPM, and the Z-score of log2 (FC) was calculated as follows:

To evaluate gene expression patterns during development, gene expression was classified into the following five categories based on any variation and inflection point: (i) ‘Stable’ expression pattern was represented by genes with a Z-score ≤2 of log2(FC), which was calculated by dividing the maximum TPM by the minimum TPM during development. (ii) ‘Increase’ expression pattern was not included in ‘Stable’ and was represented by an increase in TPM of more than 10% before the developmental stage. (iii) ‘Decrease’ expression pattern was not included in ‘Stable’ and was represented by a decrease in TPM of more than 10% before the developmental stage. (iv) ‘Temporary’ (Peak and Dimple) expression pattern was not included in ‘Stable’ and was represented by an inflection point (increase to decrease or decrease to increase) in the expression. (v) ‘Complex’ expression pattern represented a pattern other than 1–4.

2.4. Gene ontology analysis

The PANTHER classification system (PANTHER17.0) was used to assess the probability of enrichment of tissue-specific genes (Z-score ≥ 4, TPM > 10).

2.5. Analysis of detection of transcript isoforms from long-read sequencing

To detect transcript isoforms from long-read sequencing data, we compared all junctions in the long-read sequencing reads with the RefSeq transcripts, allowing for a margin of 20 base pair (bp) gaps, as previously described. Since MinION reads are inaccurate, isoforms, and splicing junctions were validated using short-read sequencing. The following filters were applied to remove low-confidence isoforms: TPM ≥ 10 and isoform read counts/total read counts assigned to the same gene >10%. Finally, we classified isoforms according to the following nine types: alternative first exon, alternative last exon, alternative 5′ splice site, alternative 3′ splice site, exon skipping, exon shuffling, intron retention, unannotated exon, and combination patterns.

2.6. Data availability

The datasets supporting the results of this study are included in this article and in its supplementary files. The sequencing data generated in this study have been submitted to DDBJ under accession number DRA014315 and have been deposited in the Wagyu Genome Database (WGDB; https://wagyu.hgc.jp/open/download/rnaseq/processed_data/, Fig. 1B) of the Japan Livestock Technology Association (Yushima, Bunyouku, Tokyo 113-0034, Japan) and managed by the WGDB Consortium. The expression data in the following files can be downloaded from the Genomic Expression Archive (GEA) at DDBJ under accession number E-GEAD-530 and the WGDB: https://wagyu.hgc.jp/open/download/rnaseq/summary/; 1_short_read_gene.txt, 2_short_read_TPM.txt, 3_short_read_tissue_specifity_Z_score.txt, 4_time_course.txt, 5_long_read_gene.txt, and 6_long_read_tissue_specifity_Z_score.txt.

3. Results and discussion

3.1. Japanese Black cattle gene expression catalogue development using short-read sequencing

We collected 124 samples from 31 major tissues at three developmental stages of the Japanese Black cattle: fetuses at 30, 60, and 90 days post-fertilization, calves at 10 and 23 days post-delivery, and adults at 788 and 1,845 days post-delivery (Supplementary Table S1). Using these tissues, we obtained 6,069,868,313 reads from 124 RNA-seq data using HiSeq3000 with an average number of reads of 48,950,551 and an average mapping rate of 97% on the ARS-UCD1.2 reference genome (Supplementary Table S2). Compared with that in other tissues, the number of reads (329,727,863) was higher and the mapping rate (77.9%) was lower in the allantois and amnions derived from fetuses 30 days of post-fertilization (Supplementary Table S2), suggesting that the quality of the RNA was low due to the long time gap following sample acquisition. We summarized the statistics describing the assembled transcripts and compared them with the annotations using Bos_taurus.ARS-UCD1.2.105. gtf as the reference file (Supplementary Table S3). Of the 181,180 transcripts, we detected 32,705 genes that corresponded to the annotation reference file (Supplementary Table S3, 1_short_read_gene.txt), of which 19,815 had coding sequences (Supplementary Table S3). To compare gene expression across developmental stages and tissues, we calculated the normalized transcript expression levels (TPM) for 19,815 coding and 12,890 non-coding genes (2_short_read_TPM.txt). Of the 19,815 coding genes, we found an average of 17,366 genes expressed at TPM > 1 and 14,667 genes expressed at TMP > 10 across seven developmental stages (Table 1). Similarly, an average of 14,821 genes were expressed at TPM > 1 and 11,035 genes were expressed at TMP > 10 across 31 tissues (Table 1). Figure 2A shows the number of genes represented in each tissue at each developmental stage (Fig. 2A, Supplementary Table S4). Despite differences in experimental conditions, the number of represented genes was similar among developmental stages and tissues (Fig. 2A), indicating that our RNA-seq dataset was suitable for comparing gene expression across developmental stages and tissues.

Table 1

Summary of represented genes with/without coding sequence across developmental stages and tissues in Japanese Black cattle using short-read sequencing

Developmental stage	Stages	Sex	No. of sample	Ave. no. of reads	No. of represented protein-coding genes (>1 TPM)	No. of represented non-protein-coding genes (>1 TPM)	No. of represented protein-coding genes (>10 TPM)	No. of represented non-protein-coding genes (>10 TPM)
Fetus	01_30d_fetal_stage	Female	2	258,077,334	14,904	2,149	10,453	543
	02_60d_fetal_stage	Female	17	48,178,503	17,143	4,293	14,420	1,052
	03_90d_fetal_stage	Male	18	52,813,165	17,736	4,749	15,025	1,134
Calf	04_10d_calf_stage	Female	21	48,034,135	17,996	5,303	15,887	1,449
	06_23d_calf_stage	Female	21	47,278,037	17,870	5,074	15,592	1,337
Adult	08_788d_adult_stage	Female	23	49,037,673	17,924	5,313	15,682	1,446
	09_1845d_adult_stage	Female	22	35,109,334	17,986	5,400	15,608	1,475
Average					17,366	4,612	14,667	1,205

Figure 2

Number of represented coding genes and tissue-specific gene expression in 124 Japanese Black cattle transcriptomes using short-read sequencing. (A) Number of represented coding genes across developmental stages and tissues. The order of the samples on the x-axis is described in Supplementary Table S4. (B) The tissue-specific expression of a coding gene. To evaluate the tissue specificity of gene expression and calculate the Z-score, we scaled the TPM of each gene to have a mean of zero and a variance of one within an individual animal. Number of represented coding genes with Z-score ≥2 and ≥4 are shown. The order of the samples on the x-axis is described in Supplementary Table S5. (C) Highly tissue-specific genes (Z-score ≥ 4) in the cerebral cortex (CYP46A1) and skeletal muscles (APOBEC2) in adults at 788 days of post-delivery, rumen (KRT4) and lymph node (CXCL9) genes in calves at 10 days of post-delivery. TPM is indicated on the y-axis; scale is the same in each gene. Blue (CYP46A1, APOBEC2 and KRT4) and orange (CXCL9) RefSeq transcripts are transcribed from the plus and the minus strand reference genome, respectively (A color version of this figure appears in the online version of this article).

3.2. Gene expression tissue specificity

We then calculated a Z-score to measure the tissue specificity of gene expression in each individual animal (Fig. 2B, Supplementary Table S5), except in the fetus 30 days of post-fertilization because of their small size and difficulty in tissue separation. Across developmental stages, a large proportion of genes were specifically expressed in the nervous, endocrine, and immune systems compared with other systems (Fig. 2B, Supplementary Table S5). Genes with Z-scores ≥4 (i.e. genes with highly tissue-specific expression) were rarely found in fetuses 60 days of post-fertilization and tended to gradually increase with development (Fig. 2B). This finding implies that the onset of tissue differentiation may occur after 60 days of fertilization in cattle (Fig. 2B, Supplementary Table S5), suggesting that failure of tissue differentiation after this point may cause embryonic mortality in cattle.

We detected tissue-specific genes for each analysed tissue based on the Z-score (3_short_read_tissue_specifity_Z_score.txt). Examples of highly tissue-specific genes (Z-score ≥ 4) are shown in Fig. 2C, including CYP46A1 (cerebral cortex), APOBEC2 (skeletal muscle and heart), KRT4 (rumen), and CXCL9 (lymph nodes). CYP46A1 is responsible for maintaining cholesterol homeostasis in the brain and its expression level is associated with brain function (Fig. 2C). APOBEC2 is exclusively expressed in the skeletal and cardiac muscles, and APOBEC2 deficiency leads to decreased body mass and myopathy (Fig. 2C). KRT4 is a keratin family gene, involved in keratinized epidermal structure development, including that of the rumen, (Fig. 2C). CXCL9 is a cytokine belonging to the CXC chemokine family that is involved in immune cell differentiation and activation (Fig. 2C). In addition to these examples, gene ontology enrichment analysis of tissue-specific genes (Z-score ≥ 4) confirmed their known tissue-associated functions (Supplementary Table S6). For instance, the cerebral cortex-specific genes were significantly enriched for synaptic signalling (P = 5.8E−24, enrichment fold = 5.5), pituitary gland-specific genes for regulating hormone levels (P = 6.23E−03, enrichment fold = 3.53), skeletal muscle-specific genes for muscle system process (P = 7.91E−17, enrichment fold = 7.57), lung-specific genes for arachidonate transport (P = 2.84E−04, enrichment fold = 14.35), rumen-specific genes for epidermis development (P = 5.31E−07, enrichment fold = 4.73), lymph node-specific genes for immune system process (P = 7.50E−18, enrichment fold = 4.07), kidney-specific genes for sodium ion transport (P = 3.43E−10, enrichment fold = 6.05), and mammary gland-specific genes regulating fat cell differentiation (P = 3.00E−05, enrichment fold = 10.03) (Supplementary Table S6). Given the consistency of our results with biological functions, our gene expression catalogue can serve as a reference for future studies seeking to understand gene functions and their variants in Japanese Black cattle.

3.3. Gene expression pattern during development

We then evaluated the changes in gene expression in each tissue during development. Gene expression patterns were classified into the following five categories based on the variation and inflection point of gene expression during development: ‘Stable’, ‘Increase’, ‘Decrease’, ‘Temporary’ (Peak and Dimple), or ‘Complex’ (Fig. 3A, Supplementary Table S7, 4_time_course.txt). For example, we classified patterns of ‘Increase’ (N = 157 genes), ‘Decrease’ (N = 119 genes), and ‘Temporary Peak’ (N = 201 genes) in the rumen during the fetal, calf, and adult stages (Fig. 3A, Supplementary Table S7). In this study, we did not detect an unambiguous ‘Temporary Dimple’ expression pattern. Examples of genes showing ‘Increase’, ‘Decrease’, and ‘Temporary Peak’ expression patterns in the rumen are shown in Fig. 3B. For instance, GSR exhibits increased expression during rumen development. GSR catalyses the reduction of glutathione disulphide to glutathione, an antioxidant essential for maintaining a reducing environment in the cell. This finding suggests that GSR may be critical for reducing oxidative stress from fermentation in the rumen of cattle as they develop. Conversely, PTN was expressed in fetuses 30, 60, and 90 days of post-fertilization, and its expression decreased as the rumen developed. PTN is a neurite outgrowth factor and mitogen, suggesting that PTN may be involved in enteric neuronal plexus development in the rumen mucosa and mucosal epithelial proliferation during early embryonic development. Some genes exhibited a temporary increase in expression at different developmental stages, and these patterns were represented as ‘Peak’ (Fig. 3A and B; Supplementary Table S7). For example, KRT13 exhibited a ‘Temporary Peak’ pattern 10 days of post-delivery. KRT13 is a keratin family protein that pairs with KRT4 and co-localizes in the non-cornified stratified epithelia. Consistent with the expression pattern of KRT13, KRT4 also exhibited a ‘Temporary Peak’ pattern 10 days of post-delivery (Fig. 3B). These results suggest that KRT4 and KRT13 may function cooperatively during rumen development in the early calf stage.

Figure 3

Classification of developmental gene expression patterns using short-read sequencing. (A) Classification of gene expression patterns of ‘Increase’ (N = 157 genes), ‘Decrease’ (N = 119 genes), and ‘Temporary Peak’ (N = 201 genes) in rumen during fetal, calf, and adult stages. Gene expression was classified into five categories based on variation and inflection point to determine the developmental expression pattern: (i) ‘Stable’ expression pattern was represented by genes with Z-score ≤2 of log2(FC), which was calculated by dividing the maximum TPM by minimum TPM during development. (ii) ‘Increase’ expression pattern was not included in ‘Stable’ and was represented by an increase in TPM of more than 10% before developmental stage. (iii) ‘Decrease’ expression pattern was not included in ‘Stable’ and was represented by a decrease in TPM of more than 10% before developmental stage. (iv) ‘Temporary’ (Peak and Dimple) expression pattern was not included in ‘Stable’ and was represented by having an inflection point (increase to decrease, or decrease to increase) in the expression. (B) GSR exhibits ‘Increase’ pattern in rumen during development. PTN exhibits ‘Decrease’ pattern in rumen during development. KRT13 and KRT4 exhibit a ‘Temporary Peak’ pattern in rumen at 10 days of post-delivery. TPM is indicated on the y-axis; scale is the same in each gene. Blue (KRT4) and orange (GSR, PTN and KRT13) RefSeq transcripts are transcribed from the plus and the minus strand reference genome, respectively (A color version of this figure appears in the online version of this article).

3.4. Expression profiles of genes related to genetic disorders and economic traits across tissues and developmental stages

The developmental gene expression profile provides valuable insights into the molecular mechanisms underlying the phenotypic expression of genetic disorders and economic traits. We observed the expression profiles of 162 causative genes for genetic disorders as listed in the Online Mendelian Inheritance in Animals (OMIA), and 12 causative genes for meat quantity and quality in the Japanese Black cattle (Supplementary Table S8 and S9). For examples of genetic disorders, we show the expression profiles of three causative genes for genetic disorders (SLC12A1, ANXA10, and MYH6) (Fig. 4, Supplementary File S2), because the developmental gene expression profile of these genes provided insight into the onset, as explained below. Recessive missense mutations in SLC12A1 are associated with hydrallantois in the Japanese Black cattle, which is the excessive accumulation of fluid in the allantoic cavity of a pregnant animal from Day 158 (mid-gestation) due to impaired Na–K–Cl cotransporter activity in the fetal kidney. However, the period of onset of this fetal disorder remains unknown. SLC12A1 was expressed in the kidney (Fig. 4A1) between 60 and 90 days of post-fertilization (Fig. 4A2), implying that renal dysfunction with hydrallantois may have already occurred at around 60 days of fetal age. We previously reported that the loss of maternal ANXA10 is associated with embryonic lethality in the Japanese Black cattle and mice.ANXA10 is only expressed in the abomasum, which is glandular stomach, in cattle from fetus to adulthood (Fig. 4B1). Thus, the mechanisms by which ANXA10 is involved in the maintenance of pregnancy are not clear. In this study, we found that ANXA10 was expressed in the uterus of adult Japanese Black cattle, albeit at a low level (Fig. 4B2), suggesting that maternal ANXA10 may be involved in fetal support. Another causative gene, MYH6, causes embryonic lethality in the Belgian Blue cattle through a recessive missense mutation. MYH6 is a myosin heavy chain protein, which is the major component of the thick filament of cardiac muscle and is necessary for cardiac muscle contraction.MYH6 is expressed in both striated muscles (i.e. heart and skeletal muscles; Fig. 4C1); however, MYH6 is expressed only in the heart during early development (Fig. 4C2). These results suggest that MYH6 deficiency may lead to embryonic lethality approximately 20 days of post-fertilization due to heart failure.,

Figure 4

Causative gene expression profiles across tissues and development for genetic disorders in the Japanese Black cattle. (A1) Expression profile of SLC12A1, a causative gene of hydrallantois in the Japanese Black cattle. (A2) The expression of SLC12A1 in the kidney during the developmental stages. (B1) Expression profile of ANXA10, a causative gene for embryonic lethality in the Japanese Black cattle. (B2) The expression of ANXA10 in the uterus during developmental stages. (C1) Expression profile of MYH6, the causative gene for embryonic lethality in the Belgian Blue cattle. (C2) The expression of MYH6 in skeletal muscles and heart during developmental stages. The order of the samples along the x-axis for A1, B1, and C1 is listed in Supplementary Table S2. TPM is indicated on the y-axis.

For example of economic traits, we showed the expression profiles of two causative genes for carcass weight in cattle (PLAG1 and NCAPG) (Fig. 4, Supplementary File S2). Researchers have repeatedly identified PLAG1 and NCAPG as major causative genes for carcass weight in the Japanese Black cattle and other cattle breeds. Notable, the allele substitution in PLAG1 and NCAPG led to an increase in the carcass weights, bringing up the weights to 28.4 and 35.3 kg, respectively. These large economic benefits have attracted the interest of researchers; however, the mechanism by which these genes produce this phenotype remains unknown. Our data showed that PLAG1 was widely expressed between 60 and 90 days of post-fertilization, but its expression decreased after birth (Supplementary File S2A). PLAG1 is a transcription factor that is essential for normal embryonic development. It regulates a subset of critical zygotic genes by binding to their promoters.NCAPG was also widely expressed between 60 and 90 days of post-fertilization (Supplementary File S2B). NCAPG is highly expressed in the small intestine, lymph nodes, thymus, ovary, and caruncle of both calves and adults. NCAPG is a subunit of the condensin complex, which is involved in chromosome condensation. Unlike the expression patterns of the causative genes for genetic disorders, the mechanisms of these adult carcass weight-related genes cannot be inferred from their expression patterns and will require additional research.

3.5. Japanese black cattle gene expression catalogue development using long-read sequencing

To analyse gene expression by long-read sequencing, we used 21 samples from a calf 10 days of post-delivery and a kidney sample from a fetus 90 days of post-fertilization. Sequencing yielded 44,620,920 reads with an average sequence length of 1,385 bp and an average mapping rate of ∼100% against the reference genome (Supplementary Table S10). We summarized the statistics describing the transcripts and compared them to the reference annotations (Supplementary Table S11). Of the 84,066 transcripts, we detected 35,979 genes that corresponded to the annotation reference file (Supplementary Table S11, 5_long_read_gene.txt), of which 21,224 had a coding sequence (Supplementary Table S11).

To compare gene expression across tissues, we determined the TPM for 21,224 coding genes and 14,755 non-coding genes. Among the 21,224 coding genes, we found an average of 16,217 genes expressed at TPM > 1, and 13,449 genes expressed at TPM > 10 across 22 tissues (Table 2). Similarly, an average of 13,923 genes were expressed at TPM > 1 and 9,896 genes were expressed at TMP > 10 across 21 tissues (Table 2). Figure 5A shows the number of genes represented in each tissue at each developmental stage (Fig. 5A, Supplementary Table S12). We then calculated the Z-score to measure the tissue specificity of gene expression in the calf at 10 days of post-delivery (Fig. 5B, Supplementary Table S13, 6_long_read_tissue_specifity_Z_score.txt). These genes were specifically expressed in the nervous, immune, and digestive systems and in the oviduct compared with other biological systems and tissues (Fig. 5B, Supplementary Table S13). For example, we found high expression of tissue-specific genes (Z-score ≥ 4) KRT4 and CXCL9 in the rumen and lymph nodes, respectively (Fig. 5C), which was consistent with the short-read sequencing results (Fig. 2C).

Table 2

Summary of represented genes with/without coding sequence across developmental stages and tissues in Japanese Black cattle using long-read sequencing

Developmental stage	Stages	Sex	No. of sample	Ave. no. of reads	No. of represented protein-coding genes (>1 TPM)	No. of represented non-protein-coding genes (>1 TPM)	No. of represented protein-coding genes (>10 TPM)	No. of represented non-protein-coding genes (>10 TPM)
Fetus	03_90d_fetal_stage	Female	1	2,317,965	14,482	2,207	11,026	370
Calf	04_10d_calf_stage	Female	21	2,014,426	17,952	5,569	15,872	1,411
Average					16,217	3,888	13,449	891

Figure 5

Number of represented coding genes and tissue-specific gene expression in 22 transcriptomes in the Japanese Black cattle using long-read sequencing. (A) Number of represented coding genes across developmental stages and tissues. The order of the samples on the x-axis is described in Supplementary Table S12. (B) The tissue-specific expression of a coding gene in each individual. To evaluate the tissue specificity of gene expression and calculate the Z-score, we scaled the TPM of each gene to have a mean of zero and a variance of one within an individual animal. Number of represented coding genes with Z-score ≥2 and ≥4 are shown. The order of the samples on the x-axis is described in Supplementary Table S13. (C) Highly tissue-specific genes (Z-score ≥ 4) of the rumen (KRT4) and lymph node (CXCL9). TPM is indicated on the y-axis; scale is the same in each gene. Blue (KRT4) and orange (CXCL9) RefSeq transcripts are transcribed from the plus and the minus strand reference genome, respectively.

3.6. Characterization of the Japanese Black cattle transcript isoforms using long-read sequencing

To detect full-length transcript isoforms in the 22 long-read sequencing samples, we classified isoform patterns into nine types: alternative first exon, alternative last exon, alternative 5′ splice site, alternative 3′ splice site, exon skipping, exon shuffling, intron retention, unannotated exon, and combination patterns. All isoforms and splicing junctions were validated by short-read sequencing. We identified 232 novel non-RefSeq isoforms (Fig. 6A and Supplementary Table S14). The isoforms have open reading frames in the range of 93–2,292 bp (Supplementary Table S14). The alternative first exon, intron retention, unannotated exon, and combination comprised the largest proportion (80%) of isoforms in the 22 samples (Fig. 6A, Supplementary Table S14). For example, three novel non-RefSeq isoforms are shown in Fig. 6B. TTC14 showed an intron-retention isoform with an intron present between exons 10 and 11 (Fig. 6B). TSEN34 had an alternative first-exon isoform (Fig. 6B). PRDX2 showed an exon 2 skipping isoform (Fig. 6B). These results indicate that long-read sequencing enables the detection of splice isoforms in the Japanese Black cattle. As it is unknown how these isoforms are expressed in different developmental stages of the Japanese Black cattle, further long-read sequencing analyses will require complete construction of transcript isoforms.

Figure 6.

Identification of novel non-RefSeq isoforms in 22 Japanese Black cattle transcriptomes using long-read sequencing. (A) The proportion of splicing events in 232 novel non-RefSeq isoforms. (B) The full-length structure of novel non-RefSeq isoforms. TTC14 (intron retention isoform), TSEN34 (alternative first exon), and PRDX2 (exon skipping). Blue (PRDX2) and orange (TTC14 and TSEN34) RefSeq transcripts are transcribed from the plus and the minus strand reference genome, respectively. Vertical magenta arrows indicate novel splice events (A color version of this figure appears in the online version of this article).

Click here for additional data file.