Anna K. Goldkamp
Yahan Li
Rocio M. Rivera
Darren E. Hagen- 1Department of Animal and Food Sciences, Oklahoma State University, Stillwater, OK, United States
- 2Division of Animal Sciences, University of Missouri, Columbia, MO, United States
Background: As couples struggle with infertility and livestock producers wish to rapidly improve genetic merit in their herd, assisted reproductive technologies (ART) have become increasingly popular in human medicine as well as the livestock industry. Utilizing ART can cause an increased risk of congenital overgrowth syndromes, such as Large Offspring Syndrome (LOS) in ruminants. A dysregulation of transcripts has been observed in bovine fetuses with LOS, which is suggested to be a cause of the phenotype. Our recent study identified variations in tRNA expression in LOS individuals, leading us to hypothesize that variations in tRNA expression can influence the availability of their processed regulatory products, tRNA-derived fragments (tRFs). Due to their resemblance in size to microRNAs, studies suggest that tRFs target mRNA transcripts and regulate gene expression. Thus, we have sequenced small RNA isolated from skeletal muscle and liver of day 105 bovine fetuses to elucidate the mechanisms contributing to LOS. Moreover, we have utilized our previously generated tRNA sequencing data to analyze the contribution of tRNA availability to tRF abundance.
Results: 22,289 and 7,737 unique tRFs were predicted in the liver and muscle tissue respectively. The greatest number of reads originated from 5′ tRFs in muscle and 5′ halves in liver. In addition, mitochondrial (MT) and nuclear derived tRF expression was tissue-specific with most MT-tRFs and nuclear tRFs derived from LysUUU and iMetCAU in muscle, and AsnGUU and GlyGCC in liver. Despite variation in tRF abundance within treatment groups, we identified differentially expressed (DE) tRFs across Control-AI, ART-Normal, and ART-LOS groups with the most DE tRFs between ART-Normal and ART-LOS groups. Many DE tRFs target transcripts enriched in pathways related to growth and development in the muscle and tumor development in the liver. Finally, we found positive correlation coefficients between tRNA availability and tRF expression in muscle (R = 0.47) and liver (0.6).
Conclusion: Our results highlight the dysregulation of tRF expression and its regulatory roles in LOS. These tRFs were found to target both imprinted and non-imprinted genes in muscle as well as genes linked to tumor development in the liver. Furthermore, we found that tRNA transcription is a highly modulated event that plays a part in the biogenesis of tRFs. This study is the first to investigate the relationship between tRNA and tRF expression in combination with ART-induced LOS.
Introduction
Assisted reproductive technologies (ART) are described as treatments that manipulate reproduction to increase chances of conception and encompasses a wide array of procedures such as in vitro fertilization (IVF), intracytoplasmic sperm injection (ICSI), and embryo transfer (Huang and Rosenwaks, 2014). ART is often used to increase genetic gain and advance reproductive potential, and its use has rapidly increased in beef and dairy cattle populations (Hansen, 2006). In fact, a 2017 report showed a dramatic shift in worldwide embryo production, in which significantly higher numbers of bovine embryos are now produced in vitro compared to in vivo (Viana and J., 2017). According to the International embryo technology society (IETS) newsletter, a record of more than 1.5 million in vitro-conceived bovine embryos were produced or collected in 2020 alone (Viana, 2021). Although in vitro embryo production has quickly become the preferred technique globally, it is important to consider the effects of in vitro procedures on genomic output.
Several studies have investigated the link between ART use and the increased occurrence of congenital overgrowth syndromes, such as Beckwith-Wiedemann syndrome (BWS) in humans and Large Offspring Syndrome (LOS) in ruminants (McEvoy et al., 2000; Young et al., 2001; Butler, 2009; Vermeiden and Bernardus, 2013; Mussa et al., 2017). LOS is often characterized by overgrowth, tongue enlargement, and abdominal wall defects (Young et al., 1998; McEvoy et al., 2000; Kohler et al., 2019). BWS shares clinical features with LOS and is also associated with an increased risk of liver tumors (hepatoblastoma) (Rump et al., 2005). Livestock are often bred for economically beneficial characteristics related to production, making LOS an issue for breeders and a source of economic loss for producers. Due to their large size, LOS offspring have an increased chance of dystocia (difficult birth) which can result in death of the calf and/or dam (Sinclair et al., 2000). In addition to cow and calf mortality, dystocia can result in financial losses associated with decreased milk production and fertility, and an increased likelihood of health issues (e.g., respiratory and digestive disorders, uterine disease, mastitis) (Dematawewa and Berger, 1997; Lombard et al., 2007; Mee, 2008; Atashi et al., 2012). However, the mechanism of ART-induced fetal overgrowth remains poorly understood.
Our previous work has detected dysregulation of transcripts and differentially methylated regions (DMRs) in LOS, and some of these regions resulted in dysregulation of imprinted loci (genes expressed in a parent-specific fashion) (Chen et al., 2013; Chen et al., 2015). We also have shown that DNA methylation is associated with a very small percent of gene misregulation in LOS individuals, suggesting other factors may be influencing gene regulation (Chen et al., 2017). Therefore, there is still a lack of clarity in diagnosis due to the variation in molecular basis and presence of major clinical symptoms. Due to their crucial role in protein synthesis, our recent study investigated tRNA expression within skeletal muscle and liver in LOS. This study revealed differential expression of tRNA genes as well as tissue- and treatment- specific tRNA transcripts with unique sequence variations (Goldkamp et al., 2022). These findings as well as the discovery of small non-coding RNAs derived from tRNAs, led us to consider the role of tRNA-derived fragments (tRFs) in LOS. To date, no study has examined the relationship between bovine tRNA expression and their processed regulatory products.
During tRNA maturation, the 5′ leader and 3′ trailer sequence of precursor tRNAs (pre-tRNAs) is cleaved by RNase Z and RNase P (Lee et al., 2009; Jarrous et al., 2022) (Figure 1). Following the addition of a 3′ CCA tail and enzymatic splicing, the mature tRNA is actively transported through the nuclear pore complex. Mature tRNAs may be cleaved through a Dicer-dependent or Dicer-independent pathway and several classes of tRFs are produced based on the tRNA cleavage position: 5′ tRFs, 3′ tRFs, internal tRFs (i-tRFs; internal fragments spanning anywhere within the tRNA), 5′ halves, and 3′ halves (Figure 1). Generally, 5′ tRFs, 3′ tRFs, and i-tRFs are 16–26 nt, whereas 5′ and 3′ halves are 27–36 nt (Lee et al., 2009). Initially considered to be random tRNA degradation products, growing evidence indicates that tRFs are an emerging class of non-coding RNAs with implications in multiple biological processes, namely regulation of protein translation (Ivanov et al., 2011; Huang et al., 2021). There are several suggested mechanisms of translational inhibition, such as disrupted ribosomal interactions through mRNA competition (Sobala and Hutvagner, 2013), displacement of initiation factors necessary for translation (Kapur et al., 2017), and recruitment of RNase Z to cleave target mRNAs (Elbarbary et al., 2009). Other studies suggest that stress granules may be formed in response to tRF-mediated inhibition of protein synthesis, which can reduce apoptosis in cancer cells (Decker and Parker, 2012; Olvedy et al., 2016). The primary pathway that is frequently suggested is similar to microRNAs, in which the tRF is loaded into an RNA-induced silencing complex (RISC) to target partially complementary mRNA (Shigematsu et al., 2014; Shigematsu and Kirino, 2015; Venkatesh et al., 2016). Various tRF types have been identified in plants, humans, and cattle with some acting as promoters of metastasis or aiding in homeostasis in humans (Green et al., 2016; Torres et al., 2019), and others for example have been reported to respond to nutritional deficiency in Arabidopsis (Hsieh et al., 2009) or Bovine Leukemia Virus in cattle (Taxis et al., 2018). Considering the alterations in tRNA expression and the dysregulation of mRNA transcripts in LOS, these tRNA genes may be selectively transcribed to give rise to unique tRF subtypes capable of targeting transcripts related to growth and/or liver tumor development. Furthermore, the diverse functions of tRFs across health states indicates a possible role in syndrome development.
FIGURE 1. tRNA and tRF biogenesis pathway and silencing mechanisms in eukaryotes. tRNAs and tRNA-derived fragments (tRFs) are produced through this pathway in eukaryotes. Arrows indicate each step and help visualize how mRNA transcripts are targeted for translational repression or degradation, as well as how they influence stress response. There are differences in cellular ribonucleases used for cleavage in bacteria and yeast that are not shown. RNA POL III, RNA Polymerase III; TSEN Complex, tRNA Splicing Endonuclease Complex; AGO1-4, Argonaute 1–4; XPO5, Exportin-5; RanGTP, GTP-bound Ras-related nuclear protein; Nsun2, NOP2/Sun RNA methyltransferase family member 2; Dnmt2, DNA methyltransferase 2; NPC, Nuclear Pore Complex; ARS, aminoacyl-tRNA synthetase; ANG, Angiogenin; Mt-DNA, Mitochondrial DNA; RISC, RNA-induced silencing complex; eIF4A, eukaryotic translation initiation factor 4A; eIF4G, eukaryotic translation initiation factor 4G. Figure created with BioRender.com.
In this study, we performed small RNA sequencing on skeletal muscle and liver samples collected from day 105 artificial insemination-conceived fetuses (AI-Control), ART-conceived bovine fetuses with a body weight above the 97th percentile relative to Control-AI (ART-LOS), and ART-conceived bovine fetuses with a body weight below the 97th percentile (ART-Normal). In addition, previously generated tRNA sequencing data was used to compare the expression of mature tRNAs and their processed regulatory products (Goldkamp et al., 2022). We detected differentially expressed tRFs due to method of conception (AI vs. ART) as well as syndrome development (ART-Normal vs. ART-LOS). Our results indicate that tRNA expression is highly dynamic based on tissue type and syndrome development. This brings the possibility that some tRNA expression can act as a means of tRF production in order to regulate gene expression. This study contributes insights on the mechanisms of tRF biogenesis and their role in targeting transcripts related to growth and development.
Materials and methods
Animals and RNA isolation
Day 105 Bos taurus indicus (B. t. indicus; Nelore breed) × Bos taurus taurus (B. t. taurus; Holstein breed) F1 fetal conceptuses were previously generated by us (Li et al., 2019a). Tissues were flash frozen in liquid nitrogen and stored at −80° until RNA extraction. Total RNA was extracted from skeletal muscle and liver tissues of F1 hybrid controls (artificial insemination; Control-AI), in vitro produced ART-Normal (similar weight as controls), and in vitro produced ART-LOS (body weight greater than 97th percentile relative to controls) using TRIzol Reagent (Invitrogen, Carlsbad, CA, United States) following the manufacturer’s instructions. Quality and concentration of the RNA samples was assessed using the Agilent TapeStation RNA ScreenTape (Agilent, Santa Clara, CA, United States) and RNA integrity numbers (RIN) for all samples were ≥ 7.4. The weights of individuals within each treatment group were compared in our previous study to classify the 97th percentile (Chen et al., 2013). The weight (in grams) and identifier of the fetuses in each treatment group are as follows: 1) Control-AI fetuses: 392 g (CON#1), 404 g (CON#2), 416 g (CON#3), and 360 g (CON #4), 2) ART-Normal fetuses: 360 g (ART#1), 376 g (ART#2), and 390 g (ART#3), 3) ART-LOS fetuses: 514 g (LOS#1), 518 g (LOS#2), and 620 g (LOS#4). Control-AI and ART-Normal body weights are very similar. The 97th percentile of Control-AI weight was selected as the threshold to characterize LOS because it has been previously used to define the equivalent overgrowth syndrome in humans, Beckwith-Wiedemann Syndrome (BWS) (Weksberg et al., 2010).
Library preparation and sequencing
Small RNA library preparation was done using the TruSeq® Small RNA Library Preparation Kit (Illumina, Inc., San Diego, CA, United States) and following the manufacturer’s instructions. 1 μg of total RNA was briefly ligated to 3′ RNA adapters with ligation buffer, RNase Inhibitor, and T4 RNA Ligase 2. After the addition of stop solution, 5′ RNA adapters were also ligated with T4 RNA Ligase 2. Reverse transcription was performed with each adapter-ligated RNA library to produce cDNA constructs. Each resulting cDNA library was amplified via Polymerase Chain Reaction (PCR). A unique RPIX was used for each library sample for multiplexed sequencing and analysis. Following PCR and before cDNA construct purification, each library was run on a High Sensitivity DNA chip (Agilent, Santa Clara, CA, United States) with expected peaks of approximately 140–160 bp. The pooled libraries were resolved on a 6% Novex TBE PAGE gel (polyacrylamide gel) and a size selection of 140–180 bp (predicted size of tRNA fragments and adapters) was performed on the gel. The purified and pooled libraries were sequenced using Illumina NextSeq 500 System High-Output Kit (Illumina, Inc., San Diego, CA, United States) and conducted by the OSU Microarray Core Facility. All samples were sequenced in one lane at the same time to prevent a batch effect. The liver small RNA-seq data was provided from our previous study (GEO database accession # GSE117015) (Li et al., 2019a) and was sequenced using Illumina NextSeq 500 System High-Output Kit, the same library preparation kit and sequencing platform as the skeletal muscle tRFs, by the University of Missouri-Columbia DNA core facility.
Processing and alignment of small RNAseq data
The raw sequence reads were filtered using the fastq-mcf command from ea-utils (version 0.148d4) in order to remove the TruSeq Small RNA adapter sequence (TGGAATTCTCGGGTGCCAAGG) (Aronesty, 2011). The adapter trimmed reads were then quality trimmed using SolexaQA++ (version 3.1.6) dynamictrim utility with a Phred cut off score of 19 (Cox et al., 2010). The quality trimmed reads were kept if they had a length of at least 13 bp or greater and were sorted using the SolexaQA++ lengthsort utility. The resulting reads were then mapped against the bovine genome ARS-UCD1.2 using the MINTmap pipeline in order to predict tRNA fragments from the small RNA-seq data (Elsik et al., 2016; Loher et al., 2017). MINTmap aligned reads to a look up table that contains unique tRF sequences ranging from 16 to 50 bp that are exclusively located in regions associated with annotated tRNA genes. The reads that mapped to bovine tRFs were measured with the default setting of MINTmap, which allows no mismatches, no insertions, and no deletions and also analyzes the whole genome to retrieve all possible alignments (Loher et al., 2017). Additionally, MINTmap outputs the parental tRNA source(s) that the tRF is potentially derived from, the tRF sequence, tRF subtype and the unique MINTplate associated with the tRF. Only the exclusive tRF expression output of unnormalized reads was used for data analysis.
Differential expression analysis
Non-linear full quantile normalization was used with the betweenLaneNormalization function on CPM transformed read counts using EDAseq v2.24.0 in order to produce PCA and RLE plots (Risso et al., 2011). Principal component analyses (PCA) and relative log expression (RLE) plots were created with the plotPCA and plotRLE function of the EDASeq package, respectively. Only tRFs that had at least 5 counts per million in all of the control, or all of the ART-normal, or at least 2 ART-LOS were considered moderately expressed and kept for DE analysis. EdgeR v3.24.3 was used to conduct a differential expression analysis and the trimmed mean of M values method (TMM) of EdgeR was used for normalization (Lun et al., 2016). A likelihood ratio test was conducted using the glmLRT function of edgeR in order to identify differential expression in skeletal muscle and liver (Control-AI vs. ART-normal, ART-normal vs. ART-LOS, and Control-AI vs. ART-LOS). Differentially expressed tRFs were defined as those with a false discovery rate (FDR) of ≤ 0.05. Heatmaps of differentially expressed tRFs were created for skeletal muscle and liver tissue to graphically represent gene expression. The normalized read counts were transformed into moderated log-counts per million and heatmaps were produced using RColorBrewer v1.1-2 and the heatmap.2 function of the gplots package v3.0.1.1.
Target prediction
RNA-seq data for Control and ART-LOS individuals from our previous study was retrieved from NCBI Gene Expression Omnibus (GEO) (accession # GSE63509) (Chen et al., 2016). RNA-seq data was used to predict potential gene candidates targeted by DE tRFs. Differential expression analysis of the RNA-seq data was done using the same method previously described for the tRF analysis.
All DE tRFs that were identified between Control vs. ART-Normal, ART-Normal vs. ART-LOS, and Control vs. ART-LOS were analyzed for target prediction. The 3′ UTR sequences of all expressed protein-coding genes in the ARS-UCD1.2 bovine genome were obtained from Ensembl Release 98 (Cunningham et al., 2019). miRanda v3.3a is a program commonly used for miRNA target prediction and was used for DE tRF target prediction in this study (Enright et al., 2003; Riffo-Campos et al., 2016). The 3′ UTR sequences of the protein coding genes were used as a reference for alignment of the DE tRF sequences with a binding score cutoff of ≥ 150 and an energy cutoff of ≤ −20 (Enright et al., 2003). Since it has been proposed that tRFs may target transcripts that are only partially complementary, unlike miRNAs, the strict parameter was not used and a partially complementary seed sequence was allowed (Martinez et al., 2017; Xiong et al., 2019). Predicted tRF targets were compared with the DE transcripts obtained from analysis of RNAseq mentioned previously. Downregulated mRNA targets were overlapped with targets of upregulated tRFs and vice versa for each treatment comparison.
Functional enrichment analysis
Enrichment analysis was done using the list of candidate gene targets of DE tRFs for each treatment comparisons and all expressed protein coding genes used as the background gene set. Gene set enrichment analysis (GSEA) of the GO terms was performed using Fisher’s exact test as implemented in R package topGO v2.42.0 (Alexa A, 2022). KEGG pathway enrichment analysis was performed using the Wilcoxon rank-sum tests via the R package KEGGREST v1.30.1 (Tenenbaum, 2021). Human Phenotype Ontology (HPO) enrichment analysis was done using the g:Profiler web server with p-values corrected by the g:SCS threshold significance criterion (Raudvere et al., 2019). For the identified candidate target genes, we used mouse mutant phenotype information and performed a mammalian phenotype enrichment analysis with the Fisher’s Exact test implemented by MamPhEA (Weng and Liao, 2010). Enrichment results with a p-value of ≤ 0.05 were classified as significant. Dot plots depicting enrichment results were created with ggplot2 package v3.2.1.
YAMAT-seq data
Mature tRNA sequencing data from our previous study was used to compare tRNA and tRF levels (Goldkamp et al., 2022). MINTmap was used to provide all possible parental tRNA sources of each tRF (Elsik et al., 2016; Loher et al., 2017). In order to evaluate the relationship between parent tRNA expression and tRF abundance, both YAMAT-seq and small RNAseq data were then merged. In an effort to not exclude any parent tRNA predicted by MINTmap, the counts for each tRF were divided by the number of parental tRNAs it was predicted to be derived from and were then log transformed. If there was no detected expression in both the parental tRNA and the tRF in any treatment group, the tRNA species was not included in the scatter plot. Scatter plots were made to show tRNA and tRF expression relative to the tRNA species with ggplot2 package v3.2.1 and by tissue type to calculate Pearson’s correlation coefficient with ggpubr package v0.2.4.
Results and discussion
Small RNA sequencing
In order to understand tRF expression in bovine fetuses with congenital overgrowth syndromes, we performed small RNA sequencing to generate tRF expression profiles in skeletal muscle and liver. This resulted in an average of 10,750,864 (76.7%) and 9,074,956 (86%) reads retained per sample for muscle and liver respectively (Supplementary Table S1). Adapter and quality trimmed reads were aligned to the ARS-UCD1.2 bovine reference genome using the MINTmap pipeline in order to predict tRNA fragments from the small RNA-seq data. A total of 936,898 and 2,854,063 reads exclusively mapped to tRFs in the skeletal muscle and liver. A lower proportion of retained reads were mapped due to MINTmap’s strategy: only exact matches are allowed, one sequence is counted once no matter how often it appears within the genome, and only tRFs that map exclusively to genomic tRNA locations are counted. Therefore, we excluded ambiguous reads that mapped to locations both within and outside of tRNA loci in order to prevent false positives.
Detection of tRNA-derived fragments
Five subtypes of mapped tRFs were predicted in muscle and liver datasets: 5′-tRF, 3′-tRF, i-tRF, 5′ half, and 3′ half. A total of 22,289 unique tRFs were predicted in the liver tissue and 7,737 unique tRFs were predicted in the muscle tissue (Supplementary Table S1). The larger number of predicted tRFs in liver could be a result of high transcriptional activity in the liver tissue. Our recent tRNA study detected a greater number of tRNA genes expressed in liver compared to muscle (487 vs. 474), which could contribute to changes in the tRF profile of each tissue (Goldkamp et al., 2022). The liver acts as a key player in nutrient metabolism and detoxification, which could result in excess transcripts in an effort to effectively regulate metabolic homeostasis. In fact, the human liver transcriptome has been described to have increased complexity and significant variability in transcript expression (Shackel et al., 2006; Bahar Halpern et al., 2015). We included a filtering step, in which tRFs with counts present in any two individuals within a tissue (n = 10) were classified as expressed and kept for analysis. This filtration step yielded a total of 13,231 tRFs in the liver and 3,508 tRFs in the muscle. Out of the 13,231 expressed tRFs in the liver, the distribution of tRFs by subtype are as follows: 11,102 i-tRFs, 1,492 5′-tRFs, 305 5′ halves, 294 3′-tRFs, and 38 3′ halves. Out of the 3,508 expressed tRFs in the muscle, there were 2,748 i-tRFs, 687 5′-tRFs, 48 5′ halves, and 25 3′-tRFs (Supplementary Table S1). i-tRFs were the most common of the list of predicted tRFs in both tissues. i-tRFs arise from a variety of positions and may be derived upstream, within, or downstream of the anticodon loop (Loher et al., 2017). This creates an opportunity for more than one i-tRF to be processed from a single tRNA molecule. Despite most of the predicted tRF species being of the i-tRF subtype in both tissues, the distribution of reads derived from a particular subtype was tissue-specific. The largest portion of reads were derived from 5′ tRFs in the muscle and 5′ halves in the liver (Figure 2A). These results are consistent with previous studies in human and mouse which reported that hematopoietic tissues, such as the liver, have greater expression of 5′ tRNA halves compared to non-hematopoietic tissues and are suggested to function as immune signaling molecules (Fu et al., 2009; Dhahbi, 2015). Consistent with these observations, we found that most transcripts ranged from 22 to 24 nt in the skeletal muscle and 33–36 nt in the liver. This represents the expected size of 5′ tRFs and 5′ halves respectively (Figure 2B). Since tRFs of this size (22–24 nt) resemble miRNAs, this could indicate a higher likelihood of association with AGO proteins for gene silencing in the skeletal muscle (Stavast and Erkeland, 2019). Because the average size of a miRNA is ∼22 nt, we were curious if any of the predicted tRF sequences aligned to known miRNAs. We retrieved the mature sequences of all annotated bovine miRNAs from miRbase and aligned the tRF sequences using blast + v2.10.1 (Camacho et al., 2009; Kozomara et al., 2019). We found that none of the tRF sequences perfectly aligned to any of the bovine miRNAs.
FIGURE 2. Quantitative analysis of tRF subtypes and size distribution. (A) Predicted tRFs were classified based on the region of the mature tRNA molecule that they are derived from across all samples (n = 10). (B) Predicted tRFs were also classified based on size across all samples (n = 10). All reads were categorized based on size or subtype and the y axis represents the percent of total CPM-normalized tRF transcript counts. The y axis sums to 100% for each tissue. Summary statistics were computed with the SummarySE function of the Rmisc package and standard error bars are shown in black.
We observed approximately 2.97% and 6.05% of the expressed transcripts in muscle and liver were derived from mitochondrial (MT) tRNAs. These observations suggest that MT-derived tRFs make a minor contribution to the tRFome. Consistent with a previous tRF study, the parental tRNA from which mitochondrial and nuclear tRFs originated, varied between muscle and liver (Telonis et al., 2019). For example, most MT-tRFs were derived from LysUUU in muscle and AsnGUU in liver, while most nuclear tRFs were derived from initiator MetCAU (iMetCAU) in muscle and GlyGCC in liver (Figures 3A,B). This finding demonstrates the unique expression profiles for nuclear- and MT-derived tRFs in the muscle and liver, which could underlie tissue-specific biological processes. In addition, we observed certain tRNAs did not produce tRFs in any of the treatment groups in the muscle or liver (AlaGGC, ArgGCG, AspAUC, CysACA, GlyACC, HisAUG, SerACU, SeCeUCA, ThrGGU, TyrAUA, and MT-IleGAU). We previously found 10 of these tRNAs (excluding MT-IleGAU) to be transcriptionally silent across all treatment groups in muscle and liver (Goldkamp et al., 2022). Despite the annotation of these silent tRNAs in the bovine assembly, a previous report has illustrated that 9 of these genes (excluding SeCeUCA and MT-IleGAU) are reportedly missing from eukaryotic, bacterial, and/or archaeal species (Ehrlich et al., 2021). This could suggest a selective pressure on anticodon bias across species. As far as we know, there are no reports of any tRNA isodecoders that do not participate in tRF biogenesis. This may indicate that some tRNA species are more resistant to processing events, and is perhaps linked to tRNA modifications offering protection from cleavage (Goll et al., 2006; Schaefer et al., 2010; Guzzi and Bellodi, 2020). Finally, PheAAA is an example of a previously identified silent isodecoder with detected tRF expression, suggesting that some tRNAs are transcribed and cleaved to solely give rise to unique tRF species (Goldkamp et al., 2022).
FIGURE 3. MT and nuclear tRF distribution. A bar graph depicting the levels of tRFs derived from (A) nuclear and (B) mitochondrial parental tRNAs across Control-AI, ART-Normal, and ART-LOS groups in muscle and liver tissue. Log transformed CPM values were used and each parental tRNA was grouped at the level of the anticodon. The SummarySE function was implemented to calculate the statistics of continuous variables by treatment group and standard error bars are shown for each anticodon in each treatment group.
Data visualization by relative log expression and principal component analysis
Relative log expression (RLE) plots were used to visualize the normalized tRF expression data across and within treatment groups (Supplementary Figure S1). Most samples were constant although there was increased variation in ART-LOS #2 in the muscle (Supplementary Figure S1A) and ART-LOS #1 in the liver (Supplementary Figure S1B). This is consistent with our previous work using tissue samples from the same ART-LOS individuals, in which genes in ART-LOS #1 in liver and ART-LOS #2 in muscle were expressed differently from other LOS individuals (Chen et al., 2015). Principal component analysis (PCA) plots show the clustering of individuals based on the normalized tRFs in muscle and liver (Supplementary Figure S2). In the muscle, the Control-AI vs. ART-Normal and ART-LOS vs. ART-Normal cluster together, yet there is no clustering in the Control-AI vs. ART-LOS groups (Supplementary Figure S2A). Similarly, Control-AI vs. ART-Normal and ART-Normal vs. ART-LOS comparisons show clustering in the liver (Supplementary Figure S2B). However, the PCA displaying all three treatment groups indicates that ART-LOS #2 clusters with the Control-AI group in the liver and away from other treatment groups in the muscle. Overall, we observed variation in tRF expression within treatment groups, particularly in ART-LOS individuals. This might be due to the nature of the syndrome, as certain LOS phenotypes differ in severity and are not always present (Chen et al., 2013). Previous reports have suggested that tRFs may be less tightly regulated than other small RNAs, due to their larger abundance and the ability of each tRF to originate from several tRNA genes (Umu et al., 2018; Veneziano et al., 2019). However, several studies demonstrate that certain tRFs describe underlying mechanisms in cellular states and disease progression (Olvedy et al., 2016; Krishna et al., 2019). Although overgrowth is one of the most common characteristics of LOS, liver tumor predisposition is variable and the classification of LOS fetuses based on body weight alone likely introduced a preference for tRF dysregulation in the muscle.
Identification of differentially expressed tRNA-derived fragments
We used EdgeR v 3.24.3 to conduct differential expression analyses of the normalized tRF read counts in muscle and liver tissue (Lun et al., 2016). We conducted DE analysis across three comparisons: Control-AI vs. ART-LOS, Control-AI vs. ART-Normal, and ART-Normal vs. ART-LOS. For Control-AI vs. ART-LOS, we identified 24 DE tRFs in muscle and detected no DE tRFs in liver (Supplementary Table S2A). For ART-Normal vs. ART-LOS, we identified 764 DE tRFs in muscle and 43 DE tRFs in liver (Supplementary Table S2B). For Control-AI vs. ART-Normal, we detected 196 DE tRFs in muscle and 44 DE tRFs in liver (Supplementary Table S2C). Few studies have investigated tRF expression in muscle and liver. While we have not fully elucidated the biological processes underlying tissue-specific tRF expression, studies have found that the most critical stage of fetal skeletal muscle development occurs during early to mid-gestation in cattle and sheep, whereas liver tumors may not have formed by day 105 of fetal development (Yan et al., 2013). More specifically, the size of the liver increases throughout development but the disease of the liver may not yet be present. These differences in tRF expression could therefore be related to the time of tissue collection, where dynamic changes in expression are occurring in the muscle but not in the liver. We consistently saw higher numbers of DE tRFs in the muscle, which could suggest that the potential for gene targeting is higher in muscle tissue due to the recruitment of small RNAs that are similar in size to miRNAs (Figure 2B). DE tRFs in Control-AI vs. ART-Normal and ART-Normal vs. ART-LOS could suggest that tRF expression can be influenced by method of conception (AI vs. ART) as well as syndrome development (ART-Normal vs. ART-LOS). Additionally, heatmaps of DE tRFs in muscle and liver were produced to visualize the degree of up and down regulation across all individuals (Supplementary Figures S3, S4). The heatmap for muscle showed consistent DE expression across all treatment groups. There was much variation within Control and ART-LOS groups, which also could be explained by fewer DE tRFs detected in that comparison. The heatmap for liver displayed little consistency in expression between treatment groups, which could be due to the assignment of treatment group based on weight.
Mature tRNAs are tightly regulated for non-canonical functions
Previously generated data from our study characterizing tRNA expression profiles in Control-AI, ART-Normal, and ART-LOS individuals was used in order to better understand the relationship between mature tRNA expression and tRF abundance (Goldkamp et al., 2022). Due to the high levels of sequence conservation across tRNA species, MINTmap can identify numerous parental tRNA sources for a single tRF. In an effort to not exclude any parent tRNA source, the counts for each tRF were divided by the number of tRNAs it was predicted to be derived from. In order to determine if there was an association between tRNA and tRF abundance, we performed a Pearson correlation analysis between tRNA and tRF expression. We found the Pearson correlation coefficients were 0.47 and 0.6 (p-value ≤ 0.05) for the muscle and liver respectively (Figures 4A,B). One explanation for these moderately positive correlation coefficients is that selective transcription of tRNA genes can bias the availability of certain mature tRNAs and ultimately the population of tRFs. These findings agree with a previous study, which reported tissue-specific modulation of tRNA transcription to support its dual function in translation as well as gene regulation by tRFs (Torres et al., 2019). This data demonstrates that tRF expression is non-random and dependent on the availability of highly regulated tRNA molecules. We acknowledge that the redundancy of tRNA genes and difficulties in efficient sequencing remains a major challenge in tRNA studies.
FIGURE 4. The relationship between muscle mature tRNA expression and tRF abundance. A scatterplot showing the correlation between tRNA (x-axis) and tRF expression (y-axis) in (A) muscle and (B) liver. Counts for the tRNA and tRF dataset were CPM normalized and log-transformed. tRNA species with no detected expression in any treatment group from either dataset were not included. The Pearson’s product moment correlation coefficient was used to determine the relationship between tRNA and tRF expression. Linear regression lines were added using the geom_smooth function of ggplot2 with the linear model argument.
Differentially expressed tRNA-derived fragments target transcripts in large offspring syndrome individuals
Target prediction was done via miRanda with the sequences of the DE tRFs as well as the 3′ UTR sequences of all expressed protein-coding genes. The RNAseq datasets from a previous LOS study were used in order to identify DE mRNA transcripts and 3′ UTR sequences were retrieved from Ensembl Release 98 for the ARS-UCD1.2 reference genome (Chen et al., 2015; Cunningham et al., 2019). DE mRNA targets with an inverse relationship to DE tRFs were overlapped for each treatment group comparison and combined in order to generate candidate gene lists for each pairwise comparison. Pairwise comparisons were used for target prediction and enrichment. However, Control-AI vs. ART-LOS in the liver was not used for further analysis because there were no statistically significant DE tRFs identified. R packages topGO v2.42.0 and KEGGREST 1.30.1 were used in order to identify functionally enriched biological processes, molecular functions, and pathways of all candidate target genes. In addition, g:Profiler was used to perform an analysis of human phenotype ontology (HPO) in order to identify enriched genes that are associated with phenotypic abnormalities in human disease (Raudvere et al., 2019). There was no significant KEGG pathway or HPO enrichment in Control-AI vs. ART-LOS in the muscle tissue. This is likely due to the low number of differentially expressed tRFs (24 DE tRFs predicted). Our enrichment analysis identified several affected biological processes, molecular functions, and signaling pathways between ART-Normal and ART-LOS groups in the muscle (Figure 5A) as well as abnormalities related to the targeted genes (Figure 5B). Certain enriched HPO terms in the muscle were related to phenotypes often observed in LOS, such as hernia of the abdominal wall and abnormality of limb bone/skeletal morphology (Figure 5B) (Li et al., 2019b). The full outputs for all performed enrichment analyses can be found in Supplementary Table S3.
FIGURE 5. Biological processes, molecular functions and KEGG pathway enriched for target genes. (A) Subset of the affected biological processes, molecular functions, and KEGG pathways associated (p-value ≤ 0.05) with gene targets in ART-Normal vs. ART-LOS muscle. Gene ontology (GO) enrichment tests were performed with TopGO using the classic algorithm with the fisher test. The pathway function of the KEGGREST package was used applying the Wilcoxon rank-sum test to calculate significantly enriched pathways. (B) The top 25 enriched Human Phenotype Ontology (HPO) terms for gene targets in ART-Normal vs. ART-LOS muscle. HPO enrichment was determined using g:Profiler with the g:SCS threshold significance criterion and HPO terms with an adjusted p-value ≤ 0.05 were considered significant. The number of genes enriched in an HPO term is shown on the x axis.
In the liver tissue, several GO terms associated with metabolic processes were enriched in Control-AI vs. ART-Normal and ART-Normal vs. ART-LOS (e.g., carbohydrate derivative metabolic process and glycoprotein metabolic process). In addition, there was an enrichment of genes related to immune response in both comparisons, such as regulation of phagocytosis, regulation of lymphocyte differentiation, and regulation of T cell differentiation (Supplementary Table S3). Previous reports in other species suggest immune cells and inflammatory responses, such as phagocytosis, have a role in the progression of tumor development (Grivennikov et al., 2010; Yang et al., 2017; Lecoultre et al., 2020). We also found that both the Wnt and cGMP-PKG signaling pathways were targeted in the liver of ART-LOS individuals. Of the enriched genes, RACK1 and MAPK3 were both upregulated in ART-LOS liver tissue and enriched in the Wnt signaling pathway and the cGMP-PKG signaling pathway respectively. RACK1 is known to negatively regulate the Wnt signaling pathway, yet activate Sonic hedgehog (Shh) signaling (Che et al., 2012; Yang et al., 2019). Activation of Shh signaling has been implicated as a potential prognosis predictor in human hepatocellular carcinoma and upregulation of RACK1 in lung cancer correlates with metastasis and tumor differentiation (Shi et al., 2012; Li and Xie, 2015). Our previous LOS study reported microRNAs targeting genes in the Wnt Signaling pathway as well, suggesting complementary mechanisms affecting gross regulators of LOS development (Li et al., 2019a). The expression of MAPK3 has been implicated in several cancer types, in which upregulation of MAPK3 correlates with tumor recurrence and poor prognosis (Du et al., 2020; Yuan et al., 2020; Xiao et al., 2021). As previously mentioned, both ART-Normal and ART-LOS groups had enrichment of processes related to tumor formation. This could be due to the variability in the presence of LOS phenotypes and the assignment of individuals to a treatment group based on weight, suggesting both ART-Normal and ART-LOS could have increased chances of tumor development in liver.
In the muscle, gene targets were enriched in GO terms related to the regulation of biological process, cell cycle regulation, and tissue-specific developmental processes (Figure 5A; Supplementary Table S3). We found SMAD1 was enriched in the regulation of biological and cellular processes and was downregulated in ART-LOS individuals. SMAD1 belongs to a family of anti-differentiation transcription factors that are critical to the bone morphogenetic protein pathway, which regulates muscle mass and regeneration (Saad et al., 2021). The inhibition of SMAD1 by microRNAs results in the promotion of skeletal muscle differentiation and regeneration (Dey et al., 2012; Saad et al., 2021). Furthermore, BMI1 was upregulated in the muscle of ART-LOS individuals. Overexpression of BMI1 in mouse mesenchymal stem cells causes an increase in body size, weight, length of tibiae, and width of the cartilaginous growth plate (Chen et al., 2019). In addition, RAI1 was downregulated in ART-LOS individuals. Changes in RAI1 dosage can have significant impacts on growth and development. For example, overexpression of RAI1 can result in extreme growth retardation, whereas haploinsufficiency of RAI1 causes increased weight and fat deposition (Girirajan et al., 2008; Alaimo et al., 2014; Falco et al., 2017). RAI1 was enriched in the HPO term, abnormal appendicular skeleton morphology, suggesting that it could be responsible for phenotypic abnormalities in bovine (Figure 5B; Supplementary Table S3). According to our analysis and previous result, we confirmed the dysregulation of genes known to be associated with overgrowth: IGF2R, GNAS, DNMT3A, and CDKN1C (Chen et al., 2015). IGF2R was downregulated in ART-LOS muscle and low levels of IGF2R has been identified in both bovine and ovine fetal overgrowth due to IVF (Young et al., 2001; Li et al., 2022). GNAS was upregulated in the ART-LOS muscle, which is consistent with hypomethylation of the GNAS loci that has been observed in Beckwith-Wiedemann syndrome (Bliek et al., 2009). Although mutations in DNMT3A are typically associated with overgrowth, we observed downregulation of this gene (O’Doherty et al., 2012). This could indicate that tRFs are capable of targeting DNA methyltransferases and modulating DNA methylation imprinting. In addition, a study using a mouse model found that embryos with CDKN1C deficiency can mimic phenotypes of BWS, such as overgrowth and abdominal wall defects (Tunster et al., 2011). Similar to this report, we observed downregulation of CDKN1C in ART-LOS individuals (Tunster et al., 2011; Robbins et al., 2012). Finally, we identified that IGF1 was upregulated in ART-LOS. Shi and colleagues determined that the inhibition of a FBXO40, a negative regulator of IGF1 signaling, resulted in elevated IGF1 levels as well as increased body size and muscle mass in mice (Shi et al., 2011). A list of tRFs that targeted the described genes above is shown in Supplementary Table S3. Although our analysis was limited to protein-coding genes for target prediction, we must recognize that long non-coding RNAs could also be regulated by small regulatory RNAs, such as miRNAs, and should be considered for future investigation (Fatica and Bozzoni, 2014).
While we were unable to determine if differences in tRF and tRF-targeted gene expression exist at distinct stages of development, we expect they do because previous work has revealed dynamic changes in the transcriptome throughout development and has provided information about the control of normal development. For example, a study quantified the mouse developmental transcriptome by applying polyA-RNAseq to tissues sampled from day 10.5 of embryogenesis to birth, revealing that transcriptomes clustered by tissue type and developmental stage (He et al., 2020). Therefore, stage-specific molecular alterations are associated with normal phenotypes. Furthermore, variations in tRF expression have been observed at different time points during mouse fetal development (Su et al., 2020). Together, this information suggests the abundance of tRFs and their gene targets may change as development progresses in Control-AI and ART-conceived individuals.
In order to address the relationship between tRF subtypes and pathways, we first investigated the subtype percentages of DE tRFs in each pairwise comparison (Supplementary Figure S5). We then evaluated tRF subtype distribution within significant pathways (Supplementary Figure S6). For example, enriched pathways in the muscle included regulation of developmental process and tissue morphogenesis. In these pathways, 5′ tRFs and i-tRFs are present in nearly equal proportions. We also looked at specific examples of gene targets, such as IGF1 and IGF2R. Finally, we evaluated tRF subtypes in immune response pathways (Regulation of phagocytosis and positive regulation of T cell differentiation), where we observe the majority of tRFs are of the 5′ tRF subtype. Although it is not clear yet, these findings suggest that 5′ tRFs and i-tRFs are the major subtypes responsible for targeting in our dataset.
Mammalian phenotype enrichment analysis (MamPhEA) based on mutant mouse phenotypes further revealed that tRF-regulated genes in muscle and liver tissue were associated with traits observed in LOS (Figure 6). In muscle, genes were enriched for abnormal birth body size and abnormal skeleton morphology (Figure 6A). In liver, genes were enriched for increased tumor growth/size and altered tumor pathology (Figure 6B). Full results can be found in Supplementary Table S3.
FIGURE 6. Mammalian Phenotype Enrichment Analysis. Enriched MGI phenotypes in ART-Normal vs. ART-LOS (A) Muscle and (B) Liver. MamPhEA outputs shows enriched phenotypes that are hierarchically structured, representing mutant phenotypes enriched with more DEGs than by chance. Level 3 shows a general phenotype class and Levels 4/5 show detailed phenotype terms nested within the Level 3 class. Enrichment was determined using the Fisher’s exact test and phenotypes with a p-value ≤ 0.05 were considered significant. Enriched terms are shown in shades of orange.
Conclusion
Overall, these data sets demonstrate that tRFs are commonly found in the muscle and liver tissue of Control-AI and ART-conceived individuals. Despite a moderate amount of variation in expression, we detected DE tRFs that may target pathways related to tumor progression or overgrowth. These outcomes provide deeper insight into the epitranscriptomic alterations that occur in ART-LOS individuals. This study is the first to examine the effect of altered tRNA availability on the differential expression of tRFs and its relationship to overgrowth syndrome.
Data availability statement
The original contributions presented in the study are publicly available. This data can be found at https://www.ncbi.nlm.nih.gov/geo/ (Accession Numbers PRJNA480853, PRJNA876238, GSE213525).
Ethics statement
The animal study was reviewed and approved by All animal procedures were performed at TransOva Genetics by veterinarians, and all procedures were approved by their animal care and use committee (Protocol number—MRP2010–001) and were conducted in a manner conforming to Trans Ova Genetics policies and procedures and the Guide for the Care and Use of Laboratory Animals.
Author contributions
AG carried out the experiment, sequencing library preparation, statistical data analysis, data interpretation, and the writing/revision of this manuscript. YL contributed to the data analysis and interpretation and revision of the manuscript. RR provided tissue samples, assisted in conceptual design, contributed to data interpretation, and revision of the manuscript. DH conceived the study, designed the experiment, was involved with the interpretation of the results and in the writing/revision of the manuscript. All authors read, revised, edited, and approved the final manuscript.
Funding
This project is supported by the Agriculture and Food Research Initiative competitive Grant Nos. 2021-67016-33417 and 2018-67015-27598 from the United States Department of Agriculture National Institute of Food and Agriculture.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.1055343/full#supplementary-material
SUPPLEMENTARY FIGURE S1 | RLE plots to visualize variation of muscle tRFs across all treatment groups before and after betweenLaneNormalization in (A) muscle and (B) liver.
SUPPLEMENTARY FIGURE S2 | PCA plots of predicted tRFs before and after betweenLaneNormalization in (A) muscle and (B) liver in different pairwise comparisons: Control vs. ART-LOS, Control vs. ART-Normal, ART-Normal vs. ART-LOS, and all three treatment groups.
SUPPLEMENTARY FIGURE S3 | Heatmaps of differentially expressed tRFs in muscle. (A) ART-Normal vs. ART-LOS, (B) Control-AI vs. ART-Normal, and (C) Control-AI vs. ART-LOS.
SUPPLEMENTARY FIGURE S4 | Heatmaps of differentially expressed tRFs in liver. (A) ART-Normal vs. ART-LOS and (B) Control-AI vs. ART-Normal.
SUPPLEMENTARY FIGURE S5 | Pie charts showing the percentage of differentially expressed tRFs belonging to each subtype between pairwise comparisons in muscle and liver.
SUPPLEMENTARY FIGURE S6 | Pie charts depicting the specific tRF subtypes targeting different pathways in (A) muscle and (B) liver.
References
Alaimo, J. T., Hahn, N. C., Hahn, N. H., Mullegama, S. V., and Elsea, S. H. (2014). Dietary regimens modify early onset of obesity in mice haploinsufficient for Rai1. PLoS One 9 (8), e105077. doi:10.1371/journal.pone.0105077
Aronesty, E. (2011). ea-utils: Command-line tools for processing biological sequencing data. Durham, NC.
Atashi, H., Abdolmohammadi, A., Dadpasand, M., and Asaadi, A. (2012). Prevalence, risk factors and consequent effect of dystocia in holstein dairy cows in Iran. Asian-Australas. J. Anim. Sci. 25 (4), 447–451. doi:10.5713/ajas.2011.11303
Bahar Halpern, K., Tanami, S., Landen, S., Chapal, M., Szlak, L., Hutzler, A., et al. (2015). Bursty gene expression in the intact mammalian liver. Mol. Cell. 58 (1), 147–156. doi:10.1016/j.molcel.2015.01.027
Bliek, J., Verde, G., Callaway, J., Maas, S. M., De Crescenzo, A., Sparago, A., et al. (2009). Hypomethylation at multiple maternally methylated imprinted regions including PLAGL1 and GNAS loci in Beckwith-Wiedemann syndrome. Eur. J. Hum. Genet. 17 (5), 611–619. doi:10.1038/ejhg.2008.233
Butler, M. G. (2009). Genomic imprinting disorders in humans: A mini-review. J. Assist. Reprod. Genet. 26 (9-10), 477–486. doi:10.1007/s10815-009-9353-3
Camacho, C., Coulouris, G., Avagyan, V., Ma, N., Papadopoulos, J., Bealer, K., et al. (2009). BLAST+: Architecture and applications. BMC Bioinforma. 10, 421. doi:10.1186/1471-2105-10-421
Che, L., Yuan, Y. H., Jia, J., and Ren, J. (2012). Activation of sonic hedgehog signaling pathway is an independent potential prognosis predictor in human hepatocellular carcinoma patients. Chin. J. Cancer Res. 24 (4), 323–331. doi:10.3978/j.issn.1000-9604.2012.10.10
Chen, G., Zhang, Y., Yu, S., Sun, W., and Miao, D. (2019). Bmi1 overexpression in mesenchymal stem cells exerts antiaging and antiosteoporosis effects by inactivating p16/p19 signaling and inhibiting oxidative stress. Stem Cells 37 (9), 1200–1211. doi:10.1002/stem.3007
Chen, Z., Hagen, D. E., Elsik, C. G., Ji, T., Morris, C. J., Moon, L. E., et al. (2015). Characterization of global loss of imprinting in fetal overgrowth syndrome induced by assisted reproduction. Proc. Natl. Acad. Sci. U. S. A. 112 (15), 4618–4623. doi:10.1073/pnas.1422088112
Chen, Z., Hagen, D. E., Ji, T., Elsik, C. G., and Rivera, R. M. (2017). Global misregulation of genes largely uncoupled to DNA methylome epimutations characterizes a congenital overgrowth syndrome. Sci. Rep. 7 (1), 12667. doi:10.1038/s41598-017-13012-z
Chen, Z., Hagen, D. E., Wang, J., Elsik, C. G., Ji, T., Siqueira, L. G., et al. (2016). Global assessment of imprinted gene expression in the bovine conceptus by next generation sequencing. Epigenetics 11 (7), 501–516. doi:10.1080/15592294.2016.1184805
Chen, Z., Robbins, K. M., Wells, K. D., and Rivera, R. M. (2013). Large offspring syndrome: A bovine model for the human loss-of-imprinting overgrowth syndrome beckwith-wiedemann. Epigenetics 8 (6), 591–601. doi:10.4161/epi.24655
Cox, M. P., Peterson, D. A., and Biggs, P. J. (2010). SolexaQA: At-a-glance quality assessment of Illumina second-generation sequencing data. BMC Bioinforma. 11 (1), 485. doi:10.1186/1471-2105-11-485
Cunningham, F., Achuthan, P, Akanni, W., Allen, J., Amode, M. R., Armean, I. M., et al. (2019). Ensembl 2019. Nucleic Acids Res. 47, D745–D751. doi:10.1093/nar/gky1113
Decker, C. J., and Parker, R. (2012). P-Bodies and stress granules: Possible roles in the control of translation and mRNA degradation. Cold Spring Harb. Perspect. Biol. 4 (9), a012286. doi:10.1101/cshperspect.a012286
Dematawewa, C. M., and Berger, P. J. (1997). Effect of dystocia on yield, fertility, and cow losses and an economic evaluation of dystocia scores for Holsteins. J. Dairy Sci. 80 (4), 754–761. doi:10.3168/jds.s0022-0302(97)75995-2
Dey, B. K., Gagan, J., Yan, Z., and Dutta, A. (2012). miR-26a is required for skeletal muscle differentiation and regeneration in mice. Genes. Dev. 26 (19), 2180–2191. doi:10.1101/gad.198085.112
Dhahbi, J. M. (2015). 5' tRNA halves: The next generation of immune signaling molecules. Front. Immunol. 6, 74. doi:10.3389/fimmu.2015.00074
Du, Y., Zhang, J., Meng, Y., Huang, M., Yan, W., and Wu, Z. (2020). MicroRNA-143 targets MAPK3 to regulate the proliferation and bone metastasis of human breast cancer cells. Amb. Express 10 (1), 134. doi:10.1186/s13568-020-01072-w
Ehrlich, R., Davyt, M., Lopez, I., Chalar, C., and Marin, M. (2021). On the track of the missing tRNA genes: A source of non-canonical functions? Front. Mol. Biosci. 8, 643701. doi:10.3389/fmolb.2021.643701
Elbarbary, R. A., Takaku, H., Uchiumi, N., Tamiya, H., Abe, M., Takahashi, M., et al. (2009). Modulation of gene expression by human cytosolic tRNase Z(L) through 5'-half-tRNA. PLoS One 4 (6), e5908. doi:10.1371/journal.pone.0005908
Elsik, C. G., Unni, D. R., Diesh, C. M., Tayal, A., Emery, M. L., Nguyen, H. N., et al. (2016). Bovine genome database: New tools for gleaning function from the Bos taurus genome. Nucleic Acids Res. 44, D834–D839. doi:10.1093/nar/gkv1077
Enright, A. J., John, B., Gaul, U., Tuschl, T., Sander, C., and Marks, D. S. (2003). MicroRNA targets in Drosophila. Genome Biol. 5 (1), R1. doi:10.1186/gb-2003-5-1-r1
Falco, M., Amabile, S., and Acquaviva, F. (2017). RAI1 gene mutations: Mechanisms of smith-magenis syndrome. Appl. Clin. Genet. 10, 85–94. doi:10.2147/TACG.S128455
Fatica, A., and Bozzoni, I. (2014). Long non-coding RNAs: New players in cell differentiation and development. Nat. Rev. Genet. 15 (1), 7–21. doi:10.1038/nrg3606
Fu, H., Feng, J., Liu, Q., Sun, F., Tie, Y., Zhu, J., et al. (2009). Stress induces tRNA cleavage by angiogenin in mammalian cells. FEBS Lett. 583 (2), 437–442. doi:10.1016/j.febslet.2008.12.043
Girirajan, S., Patel, N., Slager, R. E., Tokarz, M. E., Bucan, M., Wiley, J. L., et al. (2008). How much is too much? Phenotypic consequences of Rai1 overexpression in mice. Eur. J. Hum. Genet. 16 (8), 941–954. doi:10.1038/ejhg.2008.21
Goldkamp, A., Li, Y., Rivera, R. M., and Hagen, D. E. (2022). Characterization of tRNA expression profiles in large offspring syndrome. BMC Genomics 23 (273), 273. doi:10.1186/s12864-022-08496-7
Goll, M. G., Kirpekar, F., Maggert, K. A., Yoder, J. A., Hsieh, C. L., Zhang, X., et al. (2006). Methylation of tRNAAsp by the DNA methyltransferase homolog Dnmt2. Science 311 (5759), 395–398. doi:10.1126/science.1120976
Green, D., Fraser, W. D., and Dalmay, T. (2016). Transfer RNA-derived small RNAs in the cancer transcriptome. Pflugers Arch. 468 (6), 1041–1047. doi:10.1007/s00424-016-1822-9
Grivennikov, S. I., Greten, F. R., and Karin, M. (2010). Immunity, inflammation, and cancer. Cell. 140 (6), 883–899. doi:10.1016/j.cell.2010.01.025
Guzzi, N., and Bellodi, C. (2020). Novel insights into the emerging roles of tRNA-derived fragments in mammalian development. RNA Biol. 17 (8), 1214–1222. doi:10.1080/15476286.2020.1732694
Hansen, P. J. (2006). Realizing the promise of IVF in cattle-an overview. Theriogenology 65 (1), 119–125. doi:10.1016/j.theriogenology.2005.09.019
He, P., Williams, B. A., Trout, D., Marinov, G. K., Amrhein, H., Berghella, L., et al. (2020). The changing mouse embryo transcriptome at whole tissue and single-cell resolution. Nature 583 (7818), 760–767. doi:10.1038/s41586-020-2536-x
Hsieh, L. C., Lin, S. I., Shih, A. C. C., Chen, J. W., Lin, W. Y., Tseng, C. Y., et al. (2009). Uncovering small RNA-mediated responses to phosphate deficiency in Arabidopsis by deep sequencing. Plant Physiol. 151 (4), 2120–2132. doi:10.1104/pp.109.147280
Huang, J. Y., and Rosenwaks, Z. (2014). Assisted reproductive techniques. Methods Mol. Biol. 1154, 171–231. doi:10.1007/978-1-4939-0659-8_8
Huang, P., Tu, B., Liao, H. J., Huang, F. Z., Li, Z. Z., Zhu, K. Y., et al. (2021). Elevation of plasma tRNA fragments as a promising biomarker for liver fibrosis in nonalcoholic fatty liver disease. Sci. Rep. 11 (1), 5886. doi:10.1038/s41598-021-85421-0
Ivanov, P., Emara, M. M., Villen, J., Gygi, S. P., and Anderson, P. (2011). Angiogenin-induced tRNA fragments inhibit translation initiation. Mol. Cell. 43 (4), 613–623. doi:10.1016/j.molcel.2011.06.022
Jarrous, N., Mani, D., and Ramanathan, A. (2022). Coordination of transcription and processing of tRNA. FEBS J. 289 (13), 3630–3641. doi:10.1111/febs.15904
Kapur, M., Monaghan, C. E., and Ackerman, S. L. (2017). Regulation of mRNA translation in neurons-A matter of life and death. Neuron 96 (3), 616–637. doi:10.1016/j.neuron.2017.09.057
Kohler, S., Carmody, L., Vasilevsky, N., Jacobsen, J. O. B., Danis, D., Gourdine, J. P., et al. (2019). Expansion of the human phenotype ontology (HPO) knowledge base and resources. Nucleic Acids Res. 47 (D1), D1018–D1027. doi:10.1093/nar/gky1105
Kozomara, A., Birgaoanu, M., and Griffiths-Jones, S. (2019). miRBase: from microRNA sequences to function. Nucleic Acids Res. 47 (D1), D155–D162. doi:10.1093/nar/gky1141
Krishna, S., Yim, D. G., Lakshmanan, V., Tirumalai, V., Koh, J. L., Park, J. E., et al. (2019). Dynamic expression of tRNA-derived small RNAs define cellular states. EMBO Rep. 20 (7), e47789. doi:10.15252/embr.201947789
Lecoultre, M., Dutoit, V., and Walker, P. R. (2020). Phagocytic function of tumor-associated macrophages as a key determinant of tumor progression control: A review. J. Immunother. Cancer 8 (2), e001408. doi:10.1136/jitc-2020-001408
Lee, Y. S., Shibata, Y., Malhotra, A., and Dutta, A. (2009). A novel class of small RNAs: tRNA-derived RNA fragments (tRFs). Genes. Dev. 23 (22), 2639–2649. doi:10.1101/gad.1837609
Li, J. J., and Xie, D. (2015). RACK1, a versatile hub in cancer. Oncogene 34 (15), 1890–1898. doi:10.1038/onc.2014.127
Li, Y., Boadu, F., Highsmith, M. R., Hagen, D. E., Cheng, J., and Rivera, R. M. (2022). Allele-specific aberration of imprinted domain chromosome architecture associates with large offspring syndrome. iScience 25 (5), 104269. doi:10.1016/j.isci.2022.104269
Li, Y., Donnelly, C. G., and Rivera, R. M. (2019). Overgrowth syndrome. Vet. Clin. North Am. Food Anim. Pract. 35 (2), 265–276. doi:10.1016/j.cvfa.2019.02.007
Li, Y., Hagen, D. E., Ji, T., Bakhtiarizadeh, M. R., Frederic, W. M., Traxler, E. M., et al. (2019). Altered microRNA expression profiles in large offspring syndrome and Beckwith-Wiedemann syndrome. Epigenetics 14 (9), 850–876. doi:10.1080/15592294.2019.1615357
Loher, P., Telonis, A. G., and Rigoutsos, I. (2017). MINTmap: Fast and exhaustive profiling of nuclear and mitochondrial tRNA fragments from short RNA-seq data. Sci. Rep. 7, 41184. doi:10.1038/srep41184
Lombard, J. E., Garry, F. B., Tomlinson, S. M., and Garber, L. P. (2007). Impacts of dystocia on health and survival of dairy calves. J. Dairy Sci. 90 (4), 1751–1760. doi:10.3168/jds.2006-295
Lun, A. T., Chen, Y., and Smyth, G. K. (2016). It's DE-licious: A recipe for differential expression analyses of RNA-seq experiments using quasi-likelihood methods in edgeR. Methods Mol. Biol. 1418, 391–416. doi:10.1007/978-1-4939-3578-9_19
Martinez, G., Choudury, S. G., and Slotkin, R. K. (2017). tRNA-derived small RNAs target transposable element transcripts. Nucleic Acids Res. 45 (9), 5142–5152. doi:10.1093/nar/gkx103
McEvoy, T. G., Sinclair, K. D., Young, L. E., Wilmut, I., and Robinson, J. J. (2000). Large offspring syndrome and other consequences of ruminant embryo culture in vitro: Relevance to blastocyst culture in human ART. Hum. Fertil. 3 (4), 238–246. doi:10.1080/1464727002000199061
Mee, J. F. (2008). Prevalence and risk factors for dystocia in dairy cattle: A review. Vet. J. 176 (1), 93–101. doi:10.1016/j.tvjl.2007.12.032
Mussa, A., Molinatto, C., Cerrato, F., Palumbo, O., Carella, M., Baldassarre, G., et al. (2017). Assisted reproductive techniques and risk of beckwith-wiedemann syndrome. Pediatrics 140 (1), e20164311. doi:10.1542/peds.2016-4311
O'Doherty, A. M., O'Shea, L. C., and Fair, T. (2012). Bovine DNA methylation imprints are established in an oocyte size-specific manner, which are coordinated with the expression of the DNMT3 family proteins. Biol. Reprod. 86 (3), 67. doi:10.1095/biolreprod.111.094946
Olvedy, M., Scaravilli, M., Hoogstrate, Y., Visakorpi, T., Jenster, G., and Martens-Uzunova, E. S. (2016). A comprehensive repertoire of tRNA-derived fragments in prostate cancer. Oncotarget 7 (17), 24766–24777. doi:10.18632/oncotarget.8293
Raudvere, U., Kolberg, L., Kuzmin, I., Arak, T., Adler, P., Peterson, H., et al. (2019). g:Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res. 4, W191–W198. doi:10.1093/nar/gkz369
Riffo-Campos, A. L., Riquelme, I., and Brebi-Mieville, P. (2016). Tools for sequence-based miRNA target prediction: What to choose? Int. J. Mol. Sci. 17 (12), E1987. doi:10.3390/ijms17121987
Risso, D., Schwartz, K., Sherlock, G., and Dudoit, S. (2011). GC-content normalization for RNA-Seq data. BMC Bioinforma. 12, 480. doi:10.1186/1471-2105-12-480
Robbins, K. M., Chen, Z., Wells, K. D., and Rivera, R. M. (2012). Expression of KCNQ1OT1, CDKN1C, H19, and PLAGL1 and the methylation patterns at the KvDMR1 and H19/IGF2 imprinting control regions is conserved between human and bovine. J. Biomed. Sci. 19, 95. doi:10.1186/1423-0127-19-95
Rump, P., Zeegers, M. P., and van Essen, A. J. (2005). Tumor risk in beckwith-wiedemann syndrome: A review and meta-analysis. Am. J. Med. Genet. A 136 (1), 95–104. doi:10.1002/ajmg.a.30729
Saad, N. Y., Al-Kharsan, M., Garwick-Coppens, S. E., Chermahini, G. A., Harper, M. A., Palo, A., et al. (2021). Human miRNA miR-675 inhibits DUX4 expression and may be exploited as a potential treatment for Facioscapulohumeral muscular dystrophy. Nat. Commun. 12 (1), 7128. doi:10.1038/s41467-021-27430-1
Schaefer, M., Pollex, T., Hanna, K., Tuorto, F., Meusburger, M., Helm, M., et al. (2010). RNA methylation by Dnmt2 protects transfer RNAs against stress-induced cleavage. Genes. Dev. 24 (15), 1590–1595. doi:10.1101/gad.586710
Shackel, N. A., Seth, D., Haber, P. S., Gorrell, M. D., and McCaughan, G. W. (2006). The hepatic transcriptome in human liver disease. Comp. Hepatol. 5, 6. doi:10.1186/1476-5926-5-6
Shi, J., Luo, L., Eash, J., Ibebunjo, C., and Glass, D. J. (2011). The SCF-Fbxo40 complex induces IRS1 ubiquitination in skeletal muscle, limiting IGF1 signaling. Dev. Cell. 21 (5), 835–847. doi:10.1016/j.devcel.2011.09.011
Shi, S., Deng, Y. Z., Zhao, J. S., Ji, X. D., Shi, J., Feng, Y. X., et al. (2012). RACK1 promotes non-small-cell lung cancer tumorigenicity through activating sonic hedgehog signaling pathway. J. Biol. Chem. 287 (11), 7845–7858. doi:10.1074/jbc.M111.315416
Shigematsu, M., Honda, S., and Kirino, Y. (2014). Transfer RNA as a source of small functional RNA. J. Mol. Biol. Mol. Imaging 1 (2), 8.
Shigematsu, M., and Kirino, Y. (2015). tRNA-derived short non-coding RNA as interacting partners of Argonaute proteins. Gene Regul. Syst. Bio. 9, 27–33. doi:10.4137/GRSB.S29411
Sinclair, K. D., Young, L. E., Wilmut, I., and McEvoy, T. G. (2000). In-utero overgrowth in ruminants following embryo culture: Lessons from mice and a warning to men. Hum. Reprod. 15, 68–86. doi:10.1093/humrep/15.suppl_5.68
Sobala, A., and Hutvagner, G. (2013). Small RNAs derived from the 5' end of tRNA can inhibit protein translation in human cells. RNA Biol. 10 (4), 553–563. doi:10.4161/rna.24285
Stavast, C. J., and Erkeland, S. J. (2019). The non-canonical aspects of MicroRNAs: Many roads to gene regulation. Cells 8 (11), E1465. doi:10.3390/cells8111465
Su, Z., Frost, E. L., Lammert, C. R., Przanowska, R. K., Lukens, J. R., and Dutta, A. (2020). tRNA-derived fragments and microRNAs in the maternal-fetal interface of a mouse maternal-immune-activation autism model. RNA Biol. 17 (8), 1183–1195. doi:10.1080/15476286.2020.1721047
Taxis, T. M., Kehrli, M. E., D'Orey-Branco, R., and Casas, E. (2018). Association of transfer RNA fragments in white blood cells with antibody response to bovine leukemia Virus in holstein cattle. Front. Genet. 9, 236. doi:10.3389/fgene.2018.00236
Telonis, A. G., Loher, P., Magee, R., Pliatsika, V., Londin, E., Kirino, Y., et al. (2019). tRNA fragments show intertwining with mRNAs of specific repeat content and have links to disparities. Cancer Res. 79 (12), 3034–3049. doi:10.1158/0008-5472.CAN-19-0789
Tenenbaum, D. M. B., Keggrest: Client-Side REST access to the kyoto encyclopedia of genes and genomes (KEGG). R package version 1.34.0., 2021.
Torres, A. G., Reina, O., Stephan-Otto Attolini, C., and Ribas de Pouplana, L. (2019). Differential expression of human tRNA genes drives the abundance of tRNA-derived fragments. Proc. Natl. Acad. Sci. U. S. A. 116 (17), 8451–8456. doi:10.1073/pnas.1821120116
Tunster, S. J., Van de Pette, M., and John, R. M. (2011). Fetal overgrowth in the Cdkn1c mouse model of Beckwith-Wiedemann syndrome. Dis. Model. Mech. 4 (6), 814–821. doi:10.1242/dmm.007328
Umu, S. U., Langseth, H., Bucher-Johannessen, C., Fromm, B., Keller, A., Meese, E., et al. (2018). A comprehensive profile of circulating RNAs in human serum. RNA Biol. 15 (2), 242–250. doi:10.1080/15476286.2017.1403003
Veneziano, D., Tomasello, L., Balatti, V., Palamarchuk, A., Rassenti, L. Z., Kipps, T. J., et al. (2019). Dysregulation of different classes of tRNA fragments in chronic lymphocytic leukemia. Proc. Natl. Acad. Sci. U. S. A. 116 (48), 24252–24258. doi:10.1073/pnas.1913695116
Venkatesh, T., Suresh, P. S., and Tsutsumi, R. (2016). tRFs: miRNAs in disguise. Gene 579 (2), 133–138. doi:10.1016/j.gene.2015.12.058
Vermeiden, J. P., and Bernardus, R. E. (2013). Are imprinting disorders more prevalent after human in vitro fertilization or intracytoplasmic sperm injection? Fertil. Steril. 99 (3), 642–651. doi:10.1016/j.fertnstert.2013.01.125
Viana, J. (2021). 2020 Statistics of embryo production and transfer in domestic farm animals. Embryo Technol. Newsl. 39.
Viana, J., (2017). Statistics of embryo collection and transfer in domestic farm animals. Embryo Transf. Newsl. 36, 08–25.
Weksberg, R., Shuman, C., and Beckwith, J. B. (2010). Beckwith-Wiedemann syndrome. Eur. J. Hum. Genet. 18 (1), 8–14. doi:10.1038/ejhg.2009.106
Weng, M. P., and Liao, B. Y. (2010). MamPhEA: A web tool for mammalian phenotype enrichment analysis. Bioinformatics 26 (17), 2212–2213. doi:10.1093/bioinformatics/btq359
Xiao, H., Wang, B., Xiong, H. X., Guan, J. F., Wang, J., Tan, T., et al. (2021). A novel prognostic index of hepatocellular carcinoma based on immunogenomic landscape analysis. J. Cell. Physiol. 236 (4), 2572–2591. doi:10.1002/jcp.30015
Xiong, W., Wang, X., Cai, X., Liu, Y., Li, C., Liu, Q., et al. (2019). Identification of tRNAderived fragments in colon cancer by comprehensive small RNA sequencing. Oncol. Rep. 42 (2), 735–744. doi:10.3892/or.2019.7178
Yan, X., Zhu, M. J., Dodson, M. V., and Du, M. (2013). Developmental programming of fetal skeletal muscle and adipose tissue development. J. Genomics 1, 29–38. doi:10.7150/jgen.3930
Yang, H., Zhu, Q., Cheng, J., Wu, Y., Fan, M., Zhang, J., et al. (2019). Opposite regulation of Wnt/β-catenin and Shh signaling pathways by Rack1 controls mammalian cerebellar development. Proc. Natl. Acad. Sci. U. S. A. 116 (10), 4661–4670. doi:10.1073/pnas.1813244116
Yang, Y., Chen, L., Gu, J., Zhang, H., Yuan, J., Lian, Q., et al. (2017). Recurrently deregulated lncRNAs in hepatocellular carcinoma. Nat. Commun. 8, 14421. doi:10.1038/ncomms14421
Young, L. E., Fernandes, K., McEvoy, T. G., Butterwith, S. C., Gutierrez, C. G., Carolan, C., et al. (2001). Epigenetic change in IGF2R is associated with fetal overgrowth after sheep embryo culture. Nat. Genet. 27 (2), 153–154. doi:10.1038/84769
Young, L. E., Sinclair, K. D., and Wilmut, I. (1998). Large offspring syndrome in cattle and sheep. Rev. Reprod. 3 (3), 155–163. doi:10.1530/ror.0.0030155
Yuan, J., Dong, X., Yap, J., and Hu, J. (2020). The MAPK and AMPK signalings: Interplay and implication in targeted cancer therapy. J. Hematol. Oncol. 13 (1), 113. doi:10.1186/s13045-020-00949-4
Glossary
ART Assisted reproductive technologies
LOS Large Offspring Syndrome
BWS Beckwith-Wiedemann Syndrome
tRFs tRNA-derived fragments
MT mitochondrial
DE Differentially expressed
DMRs Differentially methylated regions
Pre-tRNAs precursor tRNAs
i-tRFs Internal-tRFs
RISC RNA-induced silencing complex
AI artificial insemination
RIN RNA integrity number
PCR Polymerase Chain Reaction
PCA Principal component analyses
RLE Relative log expression
FDR False discovery rate
GSEA Gene set enrichment analysis
HPO Human Phenotype Ontology
MamPhEA Mammalian phenotype enrichment analysis
SHH Sonic hedgehog
IVF In vitro fertilization
CDC Center for Disease Control and Prevention
SRA Sequence Read Archive
RNA POL III RNA Polymerase III
TSEN Complex tRNA Splicing Endonuclease Complex
AGO1-4 Argonaute 1-4
XPO5 Exportin-5
RanGTP GTP-bound Ras-related nuclear protein
Nsun2 NOP2/Sun RNA methyltransferase family member 2
Dnmt2 DNA methyltransferase 2
NPC Nuclear Pore Complex
ARS aminoacyl-tRNA synthetase
ANG Angiogenin
Mt-DNA Mitochondrial DNA
RISC RNA-induced silencing complex
eIF4A eukaryotic translation initiation factor 4A
eIF4G eukaryotic translation initiation factor 4G
Keywords: tRNA, tRNA fragments, bovine, large offspring syndrome, regulation
Citation: Goldkamp AK, Li Y, Rivera RM and Hagen DE (2022) Differentially expressed tRNA-derived fragments in bovine fetuses with assisted reproduction induced congenital overgrowth syndrome. Front. Genet. 13:1055343. doi: 10.3389/fgene.2022.1055343
Received: 27 September 2022; Accepted: 28 October 2022;
Published: 15 November 2022.
Edited by:
Yuan Zhou, Peking University, ChinaReviewed by:
Renhua Li, Henry M Jackson Foundation for the Advancement of Military Medicine (HJF), United StatesYongsheng Yu, Jilin Academy of Agricultural Sciences (CAAS), China
Copyright © 2022 Goldkamp, Li, Rivera and Hagen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Darren E. Hagen, ZGFycmVuLmhhZ2VuQG9rc3RhdGUuZWR1