Extensive Variation in Gene Expression is Revealed in 13 Fertility-Related Genes Using RNA-Seq, ISO-Seq, and CAGE-Seq From Brahman Cattle.

Fertility is a key driver of economic profitability in cattle production. A number of studies have identified genes associated with fertility using genome wide association studies and differential gene expression analysis; however, the genes themselves are poorly characterized in cattle. Here, we selected 13 genes from the literature which have previously been shown to have strong evidence for an association with fertility in Brahman cattle (Bos taurus indicus) or closely related breeds. We examine the expression variation of the 13 genes that are associated with cattle fertility using RNA-seq, CAGE-seq, and ISO-seq data from 11 different tissue samples from an adult Brahman cow and a Brahman fetus. Tissues examined include blood, liver, lung, kidney, muscle, spleen, ovary, and uterus from the cow and liver and lung from the fetus. The analysis revealed several novel isoforms, including seven from SERPINA7. The use of three expression characterization methodologies (5' cap selected ISO-seq, CAGE-seq, and RNA-seq) allowed the identification of isoforms that varied in their length of 5' and 3' untranslated regions, variation otherwise undetectable (collapsed as degraded RNA) in generic isoform identification pipelines. The combinations of different sequencing technologies allowed us to overcome the limitations of relatively low sequence depth in the ISO-seq data. The lower sequence depth of the ISO-seq data was also reflected in the lack of observed expression of some genes that were observed in the CAGE-seq and RNA-seq data from the same tissue. We identified allele specific expression that was tissue-specific in AR, IGF1, SOX9, STAT3, and TAF9B. Finally, we characterized an exon of TAF9B as partially nested within the neighboring gene phosphoglycerate kinase 1. As this study only examined two animals, even more transcriptional variation may be present in a genetically diverse population. This analysis reveals the large amount of transcriptional variation within mammalian fertility genes and illuminates the fact that the transcriptional landscape cannot be fully characterized using a single technology alone.

Introduction

In tropical regions, Bos taurus indicus and crosses between Bos taurus indicus and Bos taurus taurus are extensively used as they are more resistant to heat stress, diseases, and ticks. Brahman cattle are a Bos taurus indicus breed extensively raised in tropical regions, including Northern Australia, Brazil, South Asia, and North America. In tropical beef production, fertility is a major driver of profitability; fertility levels can make the difference between a profitable and non-profitable enterprise. Despite the importance of fertility in tropical beef cattle, relatively little is known about the actual genes which contribute toward the genetic variation of the trait.

The FAANG (Functional Annotation of Animal Genomes) data types aim to characterize genes within the genomes of important animal species and breeds, eventually leading to an understanding of how genetic variation translates to phenotypic variation. FAANG data include but are not limited to 1) RNA-seq (Ozsolak and Milos, 2011), which uses short-read sequencing to quantify gene expression, as well as providing some information on gene structure through the presence of intron spanning reads (Wang Z. et al., 2009); 2) ISO-seq (Gonzalez-Garay, 2016), which uses long-read sequencing technology to identify isoforms and provide some information on gene abundance; and 3) CAGE-seq (cap analysis gene expression sequencing), which uses short-read sequencing and captures the 5′ end guanosine caps of eukaryotic mRNAs to identify transcription start sites (Takahashi et al., 2012a; Takahashi et al., 2012b; Salavati et al., 2020). Through the use of these technologies in combination, it is possible to characterize the structure and abundance of genes and identify potential mechanisms, such as allele-specific expression, that can lead to biological diversity.

A number of genes relevant to Brahman fertility have been identified using genome-wide association studies (Fortes et al., 2012a; Minten et al., 2013; Mota et al., 2017; Müller et al., 2017) and expression studies (Beerda et al., 2008; Moore et al., 2016; Dias et al., 2017; Nguyen et al., 2017; Nguyen et al., 2018; Nguyen et al., 2019). But within the Brahman breed, very little work has been conducted to characterize these critical genes. Here, we apply data from three gene expression datasets generated from a Brahman cow and fetus obtained from the same animal to characterize the transcriptional variation present within these genes within the Brahman breed of cattle.

Materials and Methods

Overview

In this study, genes associated with fertility traits in the literature were characterized in two Brahman cattle. Gene sequences from Bos taurus annotated genome (Rosen et al., 2020) were located within the Brahman genome (Ross et al., 2022) using BLASTn. CAGE-seq, ISO-seq, and RNA-seq data from 10 different tissues were mapped to these genes. The tissues used were from the same animal as was used to generate the Brahman genome assembly, and tissues from a female fetus she was carrying at the time of slaughter were also taken. The expression characteristics from each of these datatypes were then examined to characterize the expression variation within each gene.

Sample Collection

The tissues were obtained post commercial slaughter of an Australian Brahman cow from a commercial abattoir. Spleen, longissiumus dorsi muscle, thyroid, ovary, kidney, uterus, lungs, blood, and liver tissues were obtained from the Brahman cow, and lung and liver tissues were obtained from its developing fetus. The tissue samples were collected post commercial slaughter and immediately snap-frozen in liquid nitrogen. The samples were transferred on dry ice and then stored at −80°C until processing.

Locating Genes in Brahman Genome

Genes that have been previously associated with fertility traits in cattle were identified in the literature. Evidence for an association between the gene and variation in fertility traits included close proximity to a genome-wide association peak or a significant result in a differential expression RNA analysis and other public information, such as a fertility-related phenotype in other species (Table 1).

TABLE 1

Genes (identified from literature) as involved in Brahman fertility.

Full Gene Name (Abbreviate gene name)	Function	Evidence for association with fertility in Brahman
Androgen receptor (AR)	Member of nuclear receptor superfamily (Fortes et al., 2012b); regulates male fertility via regulation of androgen (Brinkmann et al., 1999; Wang R.-S. et al., 2009); interruption of gene function may cause spermatogenesis impairment; feminine character development in males (humans); cancers of prostate, ovary, or breast (MacLean et al., 1995; Dowsing et al., 1999).	Positional candidate in GWAS affecting male and female traits (Lyons et al., 2014).
Insulin-like growth factor 1 (IGF1)	Aids in cell growth, differentiation, and transformation (Le Reith, 1997; Yilmaz et al., 2006), regulates release of GnRH which affects age of onset of puberty, conception rate, and maintenance of pregnancy in mammals (Simmen et al., 1993; Wilson, 1998; Velazquez et al., 2008); IGF1 concentration positively regulates scrotal circumference, motile sperm quantity, and calving rate (Yilmaz et al., 2004; Fortes et al., 2012b).	Positional candidate affecting male and female traits in GWAS including related cattle breeds (Fortes et al., 2013a; Mota et al., 2017).
Inhibin subunit alpha (INHA)	Codes for α subunit of inhibin protein (Stelzer et al., 2016); protein complex with α and β subunit negatively regulates secretion of FSH (Richards and Pangas, 2010; Kaneko, 2016; Ulloa-Aguirre et al., 2018); regulates variability of inhibin before puberty (Fortes et al., 2012b).	Candidate gene for GWAS (Fortes et al., 2012b) in males and secretion related to ovulation in females (Burns et al., 2005).
Proenkephalin (PENK)	Codes for the neurotransmitters methionine enkephalin and leucine enkephalin via proteolytic cleavage (Takahashi, 2016); negatively regulates GnRH via modulation of progesterone release (Taylor et al., 2007); negatively regulates LH expression (Malven, 1995)	Biological candidate for GWAS affecting male and female traits (Fortes et al., 2012b) and SNP detection through RNA-seq study (Dias et al., 2017).
Pleiomorphic adenoma gene 1 (PLAG1)	Family of zinc finger transcription factors (Juma et al., 2016) influences the stature of cattle which negatively affects fertility (Karim et al., 2011), height of hip, food intake, calving ease, weight at birth, further weight gain, and body size (Snelling et al., 2010; Karim et al., 2011; Pausch et al., 2011; Littlejohn et al., 2012; Fortes et al., 2013b; Juma et al., 2016; Bouwman et al., 2018).	Haplotype analysis (Utsunomiya et al., 2017) and candidate gene for GWAS affecting male and female traits (Littlejohn et al., 2012; Fortes et al., 2013b; Mota et al., 2017).
Ribosomal protein S20 (RPS20)	Family of S10P ribosomal proteins (O’Leary et al., 2016) pleiotropically affects body size (McGowan et al., 2008); candidate gene for calving ease and puberty (Pausch et al., 2011; Fortes et al., 2012a).	Functional candidate for GWAS of female traits (Fortes et al., 2012a; Fortes et al., 2013b).
Rhotekin (RTKN2)	Family of rhotekin (Collier et al., 2004); codes for rhotekin protein, a part of Rho-binding domain group of Rho-GTPase effectors (Fu et al., 2000), influences exocytosis of the pituitary gland that produces GnRH which regulated the release of LH and FSH (Marques et al., 2000; Momboisse et al., 2011); GWAS study suggests that genetic variation in GnRH stimulation might be influenced by the RTKN2 gene (Fortes et al., 2012b)	Functional candidate for GWAS of male traits (Fortes et al., 2012b).
Serine peptidase inhibitor, clade A member 7 (SERPINA7)	Family of SERPIN (O’Leary et al., 2016) codes for serum thyroid hormone precursor, thyroglobulin (Nonneman et al., 2005) indirectly influences steroid hormone production and spermatogenesis (Nonneman et al., 2005), indirectly affects sex hormone binding globulin level that affects ovarian function (Poppe et al., 2008), and influences the size of testis and sperm count and volume in boars (Dongren et al., 2006; Wagner et al., 2008; Fortes et al., 2012b).	Candidate gene for GWAS in male traits (Fortes et al., 2012b).
SRY-transcription factor 9 (SOX9)	Family of the SOX gene (Thomsen et al., 2008) affects male fertility (Thomsen et al., 2008), influences differentiation of sertoli cells in testes, neural crest cells, and chondrocytes (Thomsen et al., 2008), and is associated with testis development. Prostate development and sex determination and sex reversal of embryo (Kent et al., 1996; Vidal et al., 2001; Thomsen et al., 2008; Alankarage et al., 2016).	Candidate gene for GWAS in male traits (Soares et al., 2017).
Signal transducer and activator of transcription 3 (STAT3)	Family of STAT protein involved in the growth hormone receptor (GHR) signaling pathway regulating growth hormone or somatotropin, involved in estrogen receptor pathway, and functional disruption causes infertility, obesity, hyperphagia, and thermal dysregulation (Gao et al., 2004).	Candidate gene for GWAS (Mota et al., 2017) and differential gene expression study (Nguyen et al., 2017) for female traits.
Serine/threonine kinase 11 interacting protein (STK11IP)	Affects male fertility (Fortes et al., 2012b), affects spermatogenesis in humans, is associated with inhibin expression in Brahman bulls (Fortes et al., 2012b), and affects spermatozoa development related to histone binding (Steilmann et al., 2011).	Functional candidate for GWAS of male traits (Fortes et al., 2012b).
TATA-box binding protein associated factor 1 (TAF1)	Family of TBP-associated factors, TAFs (Freiman et al., 2001), part of TFIID, a basal transcription factor (Freiman et al., 2001), involved in development of normal follicles in the ovary (Freiman et al., 2001), influences regulation of cell cycle and differentiation of spermatids in Drosophila (Metcalf and Wassarman, 2007), and affects cattle puberty, time of fertilization, and embryo development (Hecht et al., 2011).	Positional and functional candidate for GWAS of male traits (Fortes et al., 2012a).
TATA-box binding protein associated factor 9b (TAF9B)	Family of TBP-associated factors, TAFs (Freiman et al., 2001), functions similar to TAF1 and affects scrotal development in Brahman bulls (Fortes et al., 2012a).	Positional and functional candidate for GWAS of male traits (Fortes et al., 2012a).

The coding sequence (CDS) of each target gene was downloaded from the Bos taurus gene database of National Centre for Biotechnology Information (NCBI) (Rosen et al., 2020). Bos taurus gene sequences were aligned to the genome assembly of the Brahman animal (Ross et al., 2022) using BLASTn (Zhang et al., 2000; National Center for Biotechnology Information, 2008), with an e-value cutoff of 10−10.

Short-Read Whole Genome Sequencing

To identify and confirm the genomic sequence of the genes and to identify heterozygous loci, short-read data from both the Brahman cow and the fetus that were sequenced on the Novaseq6000 S4 flow-cell on a 2 × 150 bp paired-end run were aligned to the Brahman genome assembly (Ross et al., 2022).

To remove low-quality data and adapters, the sequences were quality-trimmed before analysis. The raw sequence data were quality-trimmed using the program QUADTrim (Robinson and Ross, 2019) using the options “-m 10” to direct QUADTrim to perform quality trim and N-base filter; “-g” to remove the guanosine tail, an error that often results from the optics of the NovaSeq6000; and “-d bulls” to specify the pre-set trimming filters specified in the 1000 Bull Genome Project (Hayes and Daetwyler, 2019).

Additionally, to identify genome-wide SNP, the 1000 Bull genomes pipeline (23) was followed. Briefly, reads were aligned to the ARS1.2 Bos taurus genome assembly. Alignments were quality-filtered, and SNPs were called. Loci that were called as the homozygous reference (Bos taurus) were removed. SNPs which were homozygous alternative (Ho) and heterozygous (He) were counted. Using the assumption that the whole genome contained 2.7 Gbp, the percentage divergence (D) between the whole genome at the haplotype level was calculated as:

RNA Extraction

Total RNA was isolated using mirVana miRNA Isolation Kit (Ambion) following the manufacturer’s instruction. RNA purity was evaluated with a Nanodrop ND-1000 spectrophotometer (v.3.5.2, Thermo Fisher Scientific). QubitTM 4.0 Fluorometer with the Qubit RNA BR (broad-range) assay kit (Thermo Fisher Scientific) was used to quantify RNA concentration. The assessment of RNA integrity was performed using Agilent 2100 Bioanalyser (Agilent Technologies). Only RNA with integrity number greater than 8.0 was used for library preparation for RNA-seq, CAGE-seq, and ISO-seq.

RNA Sequencing

All RNA samples were sent to the Ramaciotti Centre for Genomics (UNSW Sydney, Australia) for library preparation and sequencing. Stranded paired-end RNA-seq libraries were sequenced on a 2 × 100 bp paired-end NovaSeq6000 run with an S4 flowcell. The resulting reads were trimmed in the same way as the whole genome sequencing Illumina data.

CAGE-seq Sequencing

The RNA was sequenced on Illumina single-read flow cells utilizing the 27-nt-long tags prepared corresponding to the 5′-end of the capped RNAs as per Salavati et al., (2020). The libraries were sequenced on an Illumina HiSeq 2500 platform (50 nt single-read) at the Centre for Genomic Research, University of Liverpool, Liverpool. After sequencing, read quality was assessed using FastQC (Andrews, 2010), and quality trimming was administered using Trimmomatic, version 0.35 (Bolger et al., 2014) using the settings “CROP:9” to trim the last nine bases and “HEADCROP:14” to trim the initial 14 base pairs (Forutan et al., 2021).

ISO-seq Sequencing

First-strand cDNA synthesis was conducted using the TeloPrime Full-Length cDNA Amplification kit (Lexogen, Australia) from 1 μg of total RNA input according to the manufacturer’s guideline, with an exemption of the uterus and ovary samples. Due to low extracted RNA concentration, only 500 ng of total RNA from these two tissues was used at this step. Additionally, among these tissue samples, the full-length double-stranded cDNAs from the fetal liver, thyroid, and spleen were prepared using the TeloPrime Full-Length cDNA Amplification version 1 kit (Lexogen, Australia). The TeloPrime Full-Length cDNA Amplification version 2 kit was used for all other samples.

To determine the optimal PCR cycle number for the large-scale PCR, the full-length double-stranded cDNAs were fist amplified in a qPCR reaction using 3′ and 5′ end-specific primers from the TeloPrime Full-Length cDNA Amplification kit (Lexogen, Australia) combined with PCR master mix reagents from PrimerSTAR GXL DNA Polymerase (Takara, Australia). SYBR Green I (Invitrogen, United States) was added to a final concentration of ×0.1 in the qPCR reaction with a total of 40 cycles. QPCR results were evaluated based on the fluorescence value, electrophoresis images, and bioanalyzer results using the Agilent DNA 12000 kit (Agilent Technologies, Germany). Large-scale PCR was then performed using optimal PCR cycles determined during the optimization step for each sample.

PCR products for all samples were sent to Ramaciotti Centre for Genomics (UNSW Sydney, Australia) for library preparation and sequencing. Briefly, PacBio IsoSeq libraries were prepared using the PacBio SMRTBell template prep kit 1.0 SPv3 for sequel protocol. Aliquots of the cDNA products underwent a ×1 and ×0.4 Ampure PB clean-up (Beckman Coulter, Australia). The aliquots were combined post clean-up using different ratios with preference to the aliquots enriched for transcripts above 4 kb. The libraries were sequenced on the PacBio Sequel system using 10-hour movies and v3.0 sequencing chemistry. A total of 11 SMRT cells were sequenced on a PacBio Sequel system (Ramaciotti Centre for Genomics, UNSW Sydney, Australia) for 10 samples. The fetal lung sample was sequenced twice as the first run was overloaded.

Demultiplexing, filtering, and quality control were performed using SMRT Link version 6.0.0 (Pacific Biosciences). The raw reads (subreads) generated by PacBio were used for calling circular consensus sequence (CCS) using the CCS tool (version 4.2.2) with parameters “—skip-polish –min-passess=3.” Adapter sequences from these CCS reads were removed using Lima tool (version 1.11.0) with parameters “lima –isoseq –dumclips”. The polyA tails and artificial concatemers were trimmed and removed using the refine tool (isoseq refine –require-polya –min-length-polya 8).

Mapping of Expression Data

All three datasets (CAGE-seq, RNA-seq, and ISO-seq) used in this study and the gene sequences were mapped to the Brahman genome that was assembled from the same animal (cow) as all the nonfetal expression data were generated from.

The CAGE-seq data were mapped using BWA-mem (Li and Durbin, 2009) optimized as per Forutan et al. (2021) with the options -M (to mark the shorter split hits as secondary hits), -k 10 (to specify that the sequences with seed length below 10 were skipped), -T 10 (to filter out alignments with a mapping score less than 10), -L 4 (to specify the clipping penalty), and -B 5 (to specify the mismatch penalty).

RNA-seq data were mapped to the Brahman genome using STAR (Dobin et al., 2013). The reference genome was indexed using the “genomeGenerate” option. The alignments were output in sorted BAM format.

Minimap2 (Li, 2018; 2019) was used to align the gene sequences and ISO-seq data. The output was in “SAM” format, and the options “-uf” was used to find canonical splicing sites GT-AG on the transcript strand, “--secondary =no” to skip the output of secondary alignments, and “ -C5 -O6,24 -B4” which is the pre-set filter for long-read splice alignment of PacBio circular consensus sequencing reads. The alignments were converted to a sorted bam file using Samtools (Li et al., 2009).

Integrative Genome Viewer (IGV) (Robinson et al., 2011) was used to visualize alignments to interpret transcription start sites, isoforms, and gene expression. The minimum mapping quality was set to 10.

Identifying Heterozygous Loci

The heterozygous loci were identified using whole genome sequencing of the same animals. Single-nucleotide polymorphisms (SNPs) were observed manually in IGV. SNPs were only considered where both alleles were observed in at least two reads.

Statistical Analysis

To calculate the statistical significance of allele-specific expression for each tissue in the RNA-seq data, Pearson’s chi-squared test (LaMorte, 2016) was administered by comparing the difference between the expected and observed value with df = 1. The formula is as follows:

For each allele i, where Oi is the observed value and Ei is the expected value for each of the two alleles.

Only loci that were heterozygous, based on the observation of both alleles in the whole genome sequence data, were tested. The ratio of the two alleles at each of the heterozygous loci tested was assumed to be 1:1. Hence, the E-value is calculated by taking the average of the observed value. Therefore, where E is the expected value for the allele i and R is the total number of reads at that locus in the RNA-seq data (sequencing depth).

To control for any sequencing bias between the two alleles, instead of assuming equal allelic ratios, the ratio of alleles observed in the whole genome sequencing data from each of the two animals was also used to determine the expected values. Hence, the expected value for each allele i was

where E is the expected value for the allele i, R is the total observations at that locus in the RNA-seq data (sequencing depth), G is the observations of allele i in the whole genome sequencing (WGS) data, and ∑G is the total number of observations in the WGS data at that locus (sequencing depth).

To compare gene expression between tissues, the number of reads from each dataset (RNA-seq, ISO-seq, and CAGE-seq) was used. The relative expression in reads per million for each gene in each tissue (RPM ) was calculated as where T is the total number of read pairs mapped to gene i in tissue j and A is the total read pairs mapped in tissue j.

Results

Identification of Important Fertility Genes in the Brahman Genome

Thirteen genes important for Brahman fertility were identified in the literature (Table 1). The CDS sequence was extracted from the Bos taurus genome and aligned to the Brahman genome assembly using BLASTn to obtain the positions of those genes within the Brahman genome (Table 2). The mean homology of the 13 genes between the Bos taurus and Brahman genome was 99.76%. The mean homology across the entire genome of the two individuals was 99.50% (Table 3). Five genes had at least one isoform that shared a 100% homology between Bos taurus and Brahman genomes. All isoforms of proenkephalin (PENK), pleiomorphic adenoma gene 1 (PLAG1), SRY-transcription factor 9 (SOX9) and TATA-box binding protein-associated factor 9b (TAF9B), and three out of nine isoforms in insulin-like growth factor 1 (IGF1) had 100% homology. The least homology was observed in the three isoforms of rhotekin (RTKN2) ranging from 99.26 to 99.33% (Table 2).

TABLE 2

List of genes with their respective positions within Bos taurus genome and Brahman genome with their alignment length and identity.

Gene	Isoforms	Bos taurus		Bos indicus			Length of gene a	Chromosome	Strand direction b	Homology (%)
		Start Position	Stop Position	Start Position	Stop Position
AR	X1	51674157	51881942	84786064	84957236	171172		X	Negative	99.96
IGF	X1	66206081	66263849	66145894	66203733	57839		5	Negative	100
	X2	66206081	66263849	66145894	66203733	57839				99.65
	X3	66206081	66261980	66145894	66199198	53304				100
	X4	66206081	66261980	66145894	66199195	53301				100
	X5	66192424	66263849	66132237	66203733	71496				99.37
	X6	66192595	66263849	66132408	66147706	15298				99.33
	X7	66192424	66261980	66132237	66147706	15469				99.76
	X8	66192424	66261980	66132237	66147706	15469				99.75
	Preprotein	66192424	66263849	66132237	66147706	15469				99.79
INHA	X1	107501844	107504762	107614722	107617642	2920		2	Positive	99.82
PENK	X1	23542677	23546157	23420291	23423774	3483		14	Negative	100
PLAG1	X1	23330541	23332546	23192090	23194095	2005		14	Negative	100
	X2	23330541	23331794	23192090	23193343	1253
RTKN2	X1	18284769	18393694	17858619	17967072	108453		28	Negative	99.3
	X2	18284769	18385527	17858619	17958855	100236				99.26
	X3	18284769	18348138	17858619	17921964	63345				99.33
SERPINA7	X1	54824445	54829800	53237946	53241468	3522		X	Positive	99.69
	Precursor	54824445	54827949	53237963	53241468	3505				99.68
SOX9	X1	58919579	58922699	60166722	60169843	3121		19	Negative	100
STAT3	X1	42419849	42450618	43645371	43676141	30770		19	Negative	99.66
	X2	42419849	42450618	43645371	43676141	30770				99.66
	X3	42421282	42450618	43647176	43676141	28965				99.64
STK11IP	X1	107521189	107535476	107634090	107648385	14295		2	Positive	99.45
	X2	107521189	107535476	107634090	107648385	11551				99.5
	X3	107521189	107533477	107634090	107646987	12297				99.46
TAF1	X1	79206804	79276760	80853773	80923837	70064		X	Negative	99.86
	X2	79206804	79276760	80853773	80923837	70064				99.86
	X3	79206804	79276760	80853773	80923837	70064				99.79
	X4	79206804	79276760	80853773	80923837	70064				99.84
	X5	79206804	79276760	80853773	80923837	70064				99.84
	X6	79197983	79276760	80857229	80923837	66608				99.87
	X7	79225602	79276760	80872738	80923837	51099				99.9
	X8	79229398	79276760	80885675	80923837	38162				99.9
	X9	79227390	79276760	80885675	80923837	38162				99.9
TAF9B	X1	74232003	74255254	74070058	74093312	23254		X	Positive	100
	X2	74232003	74255924	74070058	74093982	23924				100
	Subunit	74232003	74239263	74070058	74077335	7277				100

Length from the first base of the first exon to the last base of the last exon.

Negative = 3′ to 5′ direction; positive = 5′ to 3′ direction.

TABLE 3

Genome-wide SNPs compared to Bos taurus genome.

	Mother	Foetus
Homozygous SNP a	7762707	7734830
Heterozygous SNP	11408849	11537221
Haplotype level Homology	99.501%	99.500%

Homozygous for alternate alleles from reference.

Within the 13 genes, heterozygous sites were identified in both the Brahman cow and fetus. There were 270 and 95 SNPs identified within the coding regions of the fertility genes for the cow and fetus, respectively. 55.07% of these were located within the untranslated regions (3.01% in 5′ UTR and 52.05% in the 3′ UTR). The highest number of heterozygous sites was observed in STK11IP; no significant heterozygous sites were observed in TAF1.

Transcription Start Sites

The transcription start sites (TSSs) were identified using CAGE-seq data (Supplementary Table S1) from Forutan et al. (2021), which were remapped to the Brahman genome (Ross et al., 2022). No CAGE-seq peak was identified in the PLAG1 region. In the 12 out of 13 genes that had a TSS identified, the CAGE-seq peaks were on an average 562 base pairs upstream of the start of the coding region in the first exon in the 5′ direction (Table 4). The largest 5′ UTR was found on SERPINA7, and the smallest was in TAF9B.

TABLE 4

Position and length of 5′ untranslated region.

Gene	5′ gene position	CAGE-seq TSS peak		UTR length a
		Start position	Stop position
AR	84957236	84958298	84958477	1,241
IGF1	66203733	66203715	66204034	301
INHA	107614722	107614985	107614491	231
PENK	23423774	23424308	23424459	685
RPS20	23140852	23141036	23141219	367
RTKN2	17967072	17966877	17967512	440
SERPINA7	53237946	53236117	53236027	1,919
SOX9	60169843	60169524	60170244	401
STAT3	43676141	43676023	43676320	179
STK11IP	107634090	107634194	107633817	273
TAF1	80923837	80923816	80923981	144
TAF9B	74070058	74070080	74069977	81

Distance between the start of the coding region within the first exon and the center of the TSS peak from the CAGE-seq data.

Tissue-Specific Expression

Tissue-specific expression of the 13 fertility-related genes was examined by comparing reads per million in all three datasets. In RNA-seq data (Figure 1; Supplementary Table S1) in blood, RPS20 has the most expression followed by STAT3, TAF9B, and TAF1. IGF1, SERPINA7, and SOX9 were not expressed at detectable levels in the blood sample, and the expression of the rest of the genes was very low. RPS20 was the most expressed gene in fetal liver followed by SERPINA7 and the least being PENK, RTKN2, and INHA with rest of the gene expression being insignificant. In adult liver tissue, STAT3 was the most prominently expressed gene followed by RPS20. While other genes showed a relatively similar expression profile between the two liver life stages, the expression level of AR was much higher in liver tissue than that in fetal liver. STAT3 was the most expressed gene in thyroid followed by RPS20. The expression of genes RTKN2, STAT3, STK11IP, TAF1, and TAF9B was higher in adult lung tissue than that in fetal lung, whereas for genes IGF1, INHA, PENK, PLAG1, RPS20, and SOX9, fetal lung showed higher expression. AR and SERPINA7 had no significant expression in both these tissues. TAF9B was the most expressed in kidney tissue, while TAF1 expression is considerably lower. STAT3 and RPS20 also had significant expression in kidney tissue. Most genes were lowly expressed in muscle tissue with an exception of RPS20, STAT3, and TAF9B. The spleen and fetal lung showed the highest expression level of RPS20, followed by STAT3 and TAF9B, when other genes have considerably lower expression levels. The fetal lung expresses SOX9 relatively highly, while its expression in the spleen was very low. Overall, fetal lung and spleen had the highest expression across all of the tested genes in the RNA-seq data.

FIGURE 1

Tissue-specific expression of genes in RNA-seq, ISO-seq, and CAGE-seq data. RPM, reads per million.

In the CAGE-seq data (Figure 1; Supplementary Table S1), RPS20 was the most expressed gene with its highest expression in the fetal lung and ovary. STK11IP and TAF9B also shows significant expression in all tissues. PLAG1 had no CAGE-seq expression data, and STAT3 shows no significant expression in any of the tissues. SERPINA7 and INHA showed significant expression in only the fetal liver and ovary, respectively. IGF1 was most expressed in the uterus along with significant expression in the liver and lung tissues, ovary, spleen, and thyroid. AR showed high expression levels in the ovary and uterus as well as the kidney, liver, and spleen. SOX9 has a fairly high expression in all tissues except blood and muscle, while TAF1 is well expressed in all tissues except the adult liver. PENK is highly expressed in the uterus in addition with notable expression in the ovary, spleen, and fetal lung.

In the ISO-seq data (Figure 1; Supplementary Tables S1, S2), RPS20 was the only gene with detected expression in all tissues, while PLAG1 showed no expression in any tissue. The ovary and uterus were the tissues where the most genes were expressed despite not having the highest sequencing depth. The ISO-seq data had the lowest sequencing depth of the three technologies used. Given the lower sequencing depth of the ISO-seq data, it is very likely that there are many more uncharacterized isoforms in the data; this is a limitation of this study.

In the ISO-seq data, genes including PLAG1, SERPINA7, STAT3, SOX9, and STK11IP were highly expressed in the ovary, while AR, IGF1, PENK, RPS20, STAT3, SOX9, and TAF1 were highly expressed in the uterus. Most genes were not expressed in the muscle and kidney. Only RPS20 was expressed in muscle, and only RPS20 and TAF9B were expressed in the kidney. TAF9B was expressed in the blood, fetal lung, kidney, ovary, spleen, and thyroid. RTKN2 shows expression in the fetal liver but not in the adult liver. It is also expressed in the fetal lung, ovary, spleen, and thyroid. INHA is only expressed in the ovary, STAT3 was only observed in the uterus and STK11IP only in blood. Despite SERPINA7 being expressed only in the liver tissues, it was the most highly expressed gene in the ISO-seq dataset with its highest expression level in the fetal liver. AR was expressed in the adult liver, ovary, and uterus and SOX9 is expressed in the fetal lung, thyroid, and uterus. IGF1 was expressed in the ovary, spleen, thyroid, and uterus. The highest expression of PENK is found in the uterus, while it is also expressed in the fetal lung, ovary, and spleen at lower levels.

Overall, RPS20 was the most highly expressed gene in any tissue in the RNA-seq, ISO-seq and CAGE-seq data. When the expression level was corrected by length to reads per kilobase million (by dividing by the length of the coding sequence for each gene in kilobasepairs), RPS20 was still the most highly expressed gene in the RNA-seq data, consistent with result in the CAGE-seq data.

Allele-Specific Expression

Allele-specific expression (Figure 2) was tested where there were expression data which overlapped a heterozygous SNP in the relevant animals (the cow or the fetal sample). A total of 117 SNPs were observed in the all genes. TAF1 had no SNP identified in these two animals. Due to the large number of loci being tested, a p-value cutoff of .0001 was used. Out of the 117 SNPs tested, nine had allele-specific expression in the cow and five in the fetus (chi-squared test; p < .0001). To determine the significance of these SNPs was not an effect of the allelic bias in the RNA-seq data; the expected ratios were adjusted to reflect the allelic ratios in the WGS data. In the cow, significant allele-specific expression was still observed in AR, RPS20, SERPINA7, and TAF9B (chi-squared test; p < .01). The SNP in SOX9 was not significant when using the new WGS-based ratio. Within the fetal tissues, RPS20 and SOX9 showed allele-specific expression when compared to both the 50:50 ratio and WGS-based expected ratios. Of the cow tissues, only the liver displayed allele-specific expression of AR, while all tested tissues showed highly significant allele-specific expression in the 5′ UTR of RPS20, which were the fetal tissues. SERPINA7 showed allele-specific expression in exons 1 and 2 in the liver tissue of the cow, which was the only sample with sufficient coverage to test for ASE. The 3′ UTR of SOX9 and TAF9B had allele-specific expression in both the fetal lung and adult lung. Overall, more SNPs associated with ASE were identified in the 5′ and 3′ UTR of the investigated genes than the exons.

FIGURE 2

Allele-specific expression for each tissue with significance. Significance at p < .0001 and p < .01 indicated for reference (dashed lines). The ratio of alleles observed in whole genome sequencing (WGS) of the same animals is indicated as the background of each panel. Only loci that were heterozygous in the WGS data were tested for allele-specific expression. In the position, 1, 2, and 3 refer to exons 1, 2, and 3, respectively.

Isoform Discovery

The isoforms present in the genes were observed in the ISO-seq data with supporting evidence from RNA-seq data (Figure 3; Supplementary Tables S3–S6). For the scope of this study, the isoforms absent in the NCBI database are considered to be novel.

FIGURE 3

Representation of isoforms with length of exons (red boxes) and introns (blue lines). Positions are displayed relative to the first expressed base pair (x axis) of the gene.

Out of all examined genes, the following had no ISO-seq data mapping to the gene region: AR, PLAG1, SOX9, and STAT3. The genes PENK and INHA had only one isoform each, which were not novel. STK11IP had one isoform expressed in blood, but it was not well supported with evidence from RNA-seq.

The highest number of isoforms was observed in SERPINA7 with seven isoforms expressed in adult and fetal liver tissues. Out of these isoforms, isoforms I, II, IV, V, VI, and VII were novel (Figure 3). Genes RPS20 and IGF1 (Figure 3) have four isoforms each. Isoform 1 of RPS20 is present in the fetal liver, kidney, liver, muscle, ovary, spleen, and uterus, while isoform III is present in the kidney and IV in the liver. Isoform II is present in all tissues. Out of these isoforms I, II, and II are found to be novel. In IGF1, isoform I is found in the ovary and uterus, while isoform II is there in the spleen and thyroid. Isoforms III and IV, which are both expressed in the uterus, are novel.

The gene RTKN2 was expressed in the fetal lung, ovary, spleen, and thyroid and has three isoforms out of which all three were novel. Isoform I, with a longer first exon, was found only in the fetal lung and thyroid. Isoform II was observed in the spleen, whereas III was found in the ovary and had the first 2 exons missing (Figure 3).

A novel isoform, which missed the first seven exons when compared to the NCBI database of TAF1 isoforms, was identified for the TAF1 gene. It is expressed only in the spleen (Figure 3).

TAF9B had two isoforms. The absence of the first exon in 5′ direction when compared to the Bos taurus gene sequence made both these isoforms novel. Isoform I was present in the blood and ovary and II was observed in the fetal lung, kidney, spleen, and thyroid (Figure 3).

Interestingly, TAF9B was a partially nested gene. The last exon of TAF9B was nested in phosphoglycerate kinase 1 in blood, kidney, lung, muscle, thyroid, fetal liver, and fetal lung tissues.

Discussion

Here, we investigated expression variation of 13 fertility-associated genes in Brahman cattle using RNA-seq, ISO-seq, and CAGE-seq. Within these few genes, the variety of data available allowed us to identify previously unknown levels of variation at the genome and transcriptome level.

The homology between the Bos taurus genome annotation and Brahman genome of 99.8% reveals that on average, these genes have a variable site within their expressed regions roughly every 500 bp. This is a slightly higher level of conservation than that of the whole genome, which has 99.5% homology in both the mother and fetus based on identified SNP loci. A higher level of conservation within these genes could be expected given their role in important reproductive traits and the selective pressure in coding or regulatory regions. However, given the extremely small number of animals used in this study, this finding needs to be validated across a large cohort of genetically diverse animals.

Identification of transcription start sites (TSSs) in these genes revealed that (where detected) the TSSs were on an average 562 basepairs upstream of the start of the coding region in the first exon. 5′ UTRs are known to regulate the posttranscriptional modification (Kim et al., 1992; Lawless et al., 2009; Araujo et al., 2012). These play an important role during embryonic development (Van Der Velden and Thomas, 1999; Leppek et al., 2018). 5′ UTRs also regulate the translation of mRNAs (Van Der Velden and Thomas, 1999). Studies suggest that 5′UTR may be utilized to control the expression of genes (Halder et al., 2009). The TSSs within our investigated genes revealed variable 5′ UTR lengths (81–1,919 bp). This is within the observed ranges of 5′ UTR across mammals (Pesole et al., 2001); however, 12 out of the 13 genes had 5′ UTR longer than the reported average for other mammals (nonhuman and nonmouse), while nine out of 13 genes had longer 5′ UTRs than the average reported for humans. Compared to UTR lengths in humans observed from the study by Davuluri et al., (2000), the 5′ UTRs of these genes are substantially longer in most cases (averaging from 47 to 250 bp depending on the expression class). This suggests that these genes related to fertility may have longer than average 5′ UTRs, which may have a role in their regulation (Pesole et al., 2001). Davuluri et al. (2000) found that genes which were thought to be poorly translated tended to have longer 5′ UTRs, suggesting that the translation of these fertility-related genes may be low, especially for SERINPINA7 and AR, both of which had 5′ UTRs longer than 1KB.

The long 5′ UTR observed in AR may lead to low translation. Lyons et al. (2014) suggested that this low level of translation may be overcome by a shift in the transcription start site or alternative splicing. AR has a cluster of SNPs and indels in the 5′UTR region in cattle. Within this polymorphism cluster, there are four putative SRY-binding sites in a perfect LD with an SNP associated with scrotal circumference in Brahman bulls (Lyons et al., 2014). Additionally, a knockout study in mice found that AR also affects fertility in females (Yeh et al., 2002), suggesting that even though AR is primarily involved in male development, it may also contribute to female fertility variation.

This is only the second report of CAGE-seq data analysis in cattle. Here, we used the same dataset as used by Forutan et al. (2021) who reported the structure of TSS positioning within the bovine genome across subspecies and developmental stages. Additionally, a CAGE-seq analysis mapped, identified, and predicted novel and previously unannotated transcription start sites (TSSs) and TSS enhancer cluster (Salavati et al., 2020) in sheep. Here, we took a more specific approach to investigate a subset of genes.

The RNA-seq results revealed that the highest relative expression of these genes is observed in the spleen and fetal lung. This is consistent with the fact that most of these genes considered for this study are related to hormonal regulation and reproductive development (Kent et al., 1996; Le Reith, 1997; Brinkmann et al., 1999; Vidal et al., 2001; Yilmaz et al., 2006; Wang R.-S. et al., 2009). Expression of genes observed in ISO-seq and CAGE-seq data was not as prominent as that in RNA-seq. The sequencing depth within these datasets is much lower than that in RNA-seq, limiting their power to gene expression. It is important to note that a cDNA size selection step was included in the ISO-seq library preparation, which could skew the quantification dramatically. For the CAGE-seq results, the short-read length could limit the ability of the data to map specifically to the genome as there are regions of high homology, such as recent duplications. Therefore, RNA-seq is likely the most accurate representation of relative abundance between the tissues of the three expression datasets used in this study.

The absence of any data for PLAG1 in ISO-seq and CAGE-seq datasets could be a direct effect of the lack of deep sequencing. Genes such as AR, SERPINA7, SOX9, and STK11IP have shown very low expression in all the datasets despite being found to be expressed well in some of these tissues by previous studies (O’Leary et al., 2016). These genes are mainly associated with male fertility, while the available datasets were compiled from a female specimen (both the cow and the fetus were female); this is the most likely cause of the lower expression levels.

Allele-specific expression was identified in five out of the 13 genes. Liver and fetal liver tissues showed the highest proportion of SNPs with allele-specific expression. Allele-specific expression has been found to be abundant in dairy cattle transcriptome data (Chamberlain et al., 2015), and this analysis confirms that at least for these genes, it is likely the same in Brahman cattle also. Allele-specific expression in a direct mechanism where genetic variation can be linked to phenotypic variation for an important trait as the amount of a gene being expressed can be directly affected by which alleles are present in the animal.

A total of 24 isoforms were found within all the observed genes, out of which 16 are novel. Isoform III of RPS20 shows an unspliced intron in the ovary and uterus. It is also notable that when intron 3 is retained, intron 2 is always retained, but when intron 2 is retained, intron 3 is not necessarily always retained. Unspliced introns are thought to have a major role in gene expression regulation in plants (Pleiss et al., 2007; Syed et al., 2012), but the part they play in mammals is not fully known. Studies suggest that intron retention may regulate the production of isoforms, stability, and efficiency of translation of RNA and rapid gene expression through posttranscriptional splicing of these introns (Jacob and Smith, 2017). Multiple ISO-seq studies in different species (Xie et al., 2018; Beiki et al., 2019; Feng et al., 2019) have identified that intron retention is quite common as a source of isoform variation.

An important limitation of this study is that only two related animals were used. Given the effect that genetic variation has been observed to have of isoform expression (Garrido-Martín et al., 2021), it could be hypothesized that in a larger population of genetically diverse animals, more variation in isoforms would be observed. The economics of data generation using current technology mean that whole population isoform characterization using long-read sequencing is not currently feasible. However, as sequencing technology continues to improve and costs decrease, population scale isoform discover is likely to be a reality in the near future. Large-scale population datasets will likely show dramatically more variation that has been observed here.

An important aspect of isoform discovery en masse is that stringent filters must be placed on the data interpretation to ensure that reported isoforms are robust. This generally means that isoform calling pipelines will collapse down isoforms where the only difference is in the first or last exon or the untranslated regions. The isoform variation in TAF9B, where only the length of the untranslated region differs between the two isoforms, is one such example. Usually isoforms such as these would raise suspicion of the shorter isoform being 5′ degraded; however, as there are multiple forms of evidence used to identify these isoforms (ISO-seq and RNA-seq) we can be confident that there are, indeed, two gene forms being expressed. This particular example highlights the usefulness of using multiple data types to identify expression variation.

Despite the absence of the last exon in TAF9B in the ISO-seq data, RNA-seq data suggest the presence of an exon downstream in 3′ direction in blood, fetal lung, kidney, ovary, spleen, and thyroid tissues. This is confirmed by a large number of intron spanning RNA-seq reads. This suggests the presence of an exon within another gene (phosphoglycerate kinase 1) present right next to TAF9B on the reverse strand. The relevance of this exon on the regulation of expression of both TAF9B and the other gene needs to be considered. Studies of nested genes, where genes were fully located within other genes (Yu et al., 2005) showed that the nested genes were under strong selection and displayed reciprocal expression with each other as well as strong tissue-specific expression. We did not observe strong tissue-specific expression to TAF9B, and the partially-nested gene it is associated with was also expressed in six of the tissues where TAF9B was observed, suggesting that TAF9B does not follow the same pattern of expression.

Genes relevant to fertility in Brahman were identified and shown to have tissue-specific expression, allele-specific expression, variation in transcription start sites, untranslated regions, and novel isoforms. This case study is an example of the detailed information that can be obtained from combining information from multiple expression datasets. It is clear that no one datatype is able to fully characterize the transcriptome, and efforts to strategically align data generation efforts will be most beneficial.