Characterization and functional roles of paternal RNAs in 2–4 cell bovine embryos

Embryos utilize oocyte-donated RNAs until they become capable of producing RNAs through embryonic genome activation (EGA). The sperm’s influence over pre-EGA RNA content of embryos remains unknown. Recent studies have revealed that sperm donate non-genomic components upon fertilization. Thus, sperm may also contribute to RNA presence in pre-EGA embryos. The first objective of this study was to investigate whether male fertility status is associated with the RNAs present in the bovine embryo prior to EGA. A total of 65 RNAs were found to be differentially expressed between 2–4 cell bovine embryos derived from high and low fertility sires. Expression patterns were confirmed for protein phosphatase 1 regulatory subunit 36 (PPP1R36) and ataxin 2 like (ATXN2L) in three new biological replicates. The knockdown of ATXN2L led to a 22.9% increase in blastocyst development. The second objective of this study was to characterize the parental origin of RNAs present in pre-EGA embryos. Results revealed 472 sperm-derived RNAs, 2575 oocyte-derived RNAs, 2675 RNAs derived from both sperm and oocytes, and 663 embryo-exclusive RNAs. This study uncovers an association of male fertility with developmentally impactful RNAs in 2–4 cell embryos. This study also provides an initial characterization of paternally-contributed RNAs to pre-EGA embryos. Furthermore, a subset of 2–4 cell embryo-specific RNAs was identified.


Gene expression analysis of 2-4 cell embryos.
A total of 65 genes were differentially expressed between embryos derived from high and low fertility sires. There were 35 genes with decreased expression and 30 genes with increased expression in the embryos from high fertility sires compared to the embryos from low fertility sires (FDR < 0.1). The subset of genes selected for validation were ATXN2L, PPP1R36, amyloid beta precursor protein binding family A member 1 (APBA1), plakophilin 2 (PKP2), cadherin EGF LAG seven-pass G-type receptor 3 (CELSR3), ELAV like RNA binding protein 4 (ELAVL4), WD repeat domain 93 (WDR93), Drebrin Like (DBNL), and ENSBTAG00000046943 (LOC112444303 oncomodulin). The genes PPP1R36, PKP2, CELSR3, and ELAVL4 were chosen because they had previously-identified associations with reproductive development and fertility. On the other hand, ATXN2L, APBA1, WDR93, DBNL, and LOC112444303 oncomodulin were selected because their roles in fertility were still unclear. In total, nine genes were chosen for qRT-PCR validation ( Table 1). All three validation replicates for ATXN2L and PPP1R36 reproduced the same trend found in RNA-Seq results. The gene ATXN2L had an average fold difference of 2.45 (P = 0.052), and PPP1R36 had an average fold difference of 4.08 (P = 0.065). The genes WDR93, ELAVL4, and LOC112444303 oncomodulin showed replication of RNA-Seq results for two out of three biological validation replicates. WDR93 had an average fold difference of 0.65 (P = 0.39), ELAVL4 had an average fold difference of 0.43 (P = 0.13), and LOC112444303 oncomodulin had an average fold difference of 0.70  Table 1. qRT-PCR results for nine genes identified as differentially-expressed between high and low fertility sires.
(P = 0.40). The genes DBNL, CELSR3, APBA1, and PKP2 were not consistent across biological validation replicates. The 65 differentially expressed genes are found in Supplementary Table S1. Additionally, there was no significant difference in the cleavage rate between high and low fertility sire embryos (P = 0.70) (Supplementary Table S2).

Knockdown of ATXN2L with gapmer treatment.
To test the possible role of ATXN2L in embryo development, knockdown experiments were performed using gapmer technology. Embryos which were treated with ATXN2L gapmer at the zygote stage showed a 22.9% increase in blastocyst rate compared to control embryos (P < 0.05). The fold change in ATXN2L was 0.44 in treated vs. control embryos (P < 0.05).
origin of RnAs. The parental origin of 6385 RNAs from 2-4 cell embryos was characterized (Fig. 1).
There were 2675 RNAs derived from both the sperm and the oocyte. An additional 2575 RNAs were exclusively oocyte-derived, while 472 RNAs were exclusively sperm-derived (Supplementary Table S3). Finally, there were also 663 RNAs that were found only in 2-4 cell embryos.
The expression levels of 24 of the differentially expressed RNAs (FDR < 0.1) in the 2-4 cell embryos was further determined. Of these, eight RNAs were downregulated and 16 RNAs were upregulated in embryos from high fertility sires versus low fertility sires. A total of four differentially expressed genes were sperm-derived. There were also two differentially expressed genes found exclusively in 2-4 cell embryos and nine that were oocyte-derived. Additionally, the sperm and oocyte co-contributed nine differentially expressed genes (Fig. 2). The parental origins of the remaining 41 differentially expressed genes were considered provisional and can be found in Supplementary Table S2.

Discussion
The RNAs donated by the oocyte at fertilization are known to affect the developmental competence of embryos [4][5][6][7][8] . It is also recognized that the sperm donates a variety of non-genetic elements [12][13][14][15][16][17][18][19][20][21] during fertilization that can alter RNA expression [22][23][24][25][26] . This study demonstrates that male fertility is associated with RNA content in the pre-EGA embryo. The functionality of one of these RNAs was interrogated through knockdown experiments. Additionally,   Supplementary  Table S3. this study reveals a substantial number of RNAs that are directly contributed from the sperm to the embryo at fertilization. Furthermore, a suite of transcripts that are specific to the 2-4 cell embryo have also been identified.
Several of the differentially expressed genes identified through this study have been previously associated with fertility. For example, a genome-wide association study found that PKP2 and CTTNBP2NL N-terminal like (CTTNBP2NL) are associated with conception rates in Holstein cows 40 . Knockdown of PKP2 and overexpression of CTTNBP2NL each resulted in reduced embryo implantation rates in mice 40 . In this study, PKP2 was upregulated and CTTNBP2NL was downregulated in embryos from high fertility sires. Additionally, both genes were expressed in sperm. Therefore, this study demonstrates that there is paternal control over the presence of these RNAs in pre-EGA embryos. Further, the direction of differential expression for these RNAs in the pre-EGA embryos reinforces previous findings about their functions in fertility.
Additional differentially expressed genes have been associated with the regulation of spermatogenesis. For instance, the gene PPP1R36 drives autophagy during spermatogenesis 41 . Inhibiting autophagy in mice leads to infertile spermatozoa 42 . PPP1R36 was found in spermatozoa in this study. Additionally, PPP1R36 was upregulated in embryos from high fertility sires. Sperm-specific knockout mice for the autophagy regulator autophagy related 5 (ATG5) generate embryos that fail to develop beyond the 4-8 cell stage 43 . Thus, future studies should evaluate whether PPP1R36 is transferred to the zygote to help regulate autophagy. The mechanistic target of rapamycin kinase (mTOR) is another spermatogenesis-regulating gene that was differentially expressed in pre-EGA embryos. mTOR can shape male reproductive potential by affecting spermatogonial stem cell maintenance, Sertoli cell physiology, and blood-testis barrier architecture 44 . Blocking mTOR activity in mouse Sertoli cells is detrimental to sperm quality and causes infertility 45 . mTOR was upregulated in embryos of high fertility sires, and was also found in sperm. Future studies should evaluate whether these genes serve a dual purpose in regulating both spermatogenesis and early embryonic development.
Although many of the differentially expressed genes in this study were sperm-derived, there was a subset of these genes that were either oocyte-derived or embryo-specific. Therefore, alternative mechanisms for paternal regulation of pre-EGA RNA may be possible. Previously, our lab showed that differential DNA methylation patterns in sperm from high and low fertility sires was correlated with altered expression of the same genes in blastocysts from high and low fertility sires 27 . In this study, the paternally-expressed imprinted genes APBA1 and GATA binding protein 3 (GATA3) 46 were both upregulated in embryos from high fertility sires. Altered DNA methylation in sperm of these genes is associated with motility 47 and male fertility status 48 , respectively. APBA1 transcripts were found in sperm, but GATA3 was categorized as oocyte-derived in this study. Future studies should evaluate alternative sperm-derived contributions that may regulate pre-EGA embryo RNA, such as DNA methylation.
The knockdown of ATXN2L led to an increase in blastocyst development. Intriguingly, the expression of this gene was elevated in 2-4 cell embryos from high fertility sires. However, ATXN2L can have equally dynamic effects when its expression is increased or decreased 49 . ATXN2L is considered a regulator of stress granule and processing body (P-body) formation 49 . When ATXN2L is upregulated in HeLa cells, P-body numbers are reduced and stress granule numbers are increased 49 . Conversely, downregulation of ATXN2L in HeLa cells causes a reduction of stress granules 49 . P-bodies are present in non-stressed cells and contain non-translating mRNA and proteins which facilitate mRNA decay and translational repression 50 . Stress granules are formed when blocks of translation initiation occur in cells, and they contain various translation initiation factors and regulators of mRNA stability 50 . Control over RNA presence and activity is crucial to EGA 51 . Therefore, ATXN2L may require a careful balance during early embryonic development. Previously, a "maternal mode" and a "zygotic mode" have been proposed to initiate mRNA decay during EGA 51 . In the future, the use of ATXN2L as a "paternal mode" to regulate the mRNA presence should be explored. It is unclear whether the knockdown of ATXN2L is truly beneficial for embryonic development. However, our results indicate that reducing the expression of this gene can upregulate blastocyst development. Future studies should investigate the role of ATXN2L when it is overexpressed at the zygote stage. Moreover, it will be important to elucidate its possible role in regulating maternal RNAs. ATXN2L was not proven to be directly contributed by sperm in this study. Therefore, the paternal mechanism of control over ATXN2L expression should also be further assessed.
Delivery of individual sperm RNAs has been demonstrated in a variety of species [14][15][16][17][18] . The CLU RNA was shown to be transferred from sperm to oocytes in multiple studies 14,16 . CLU was found in sperm, oocytes, and 2-4 cell embryos in our dataset. This indicates it is both maternally and paternally-contributed in bovine. Additionally, calcium binding tyrosine phosphorylation regulated (CABYR) was previously identified as a sperm-specific RNA delivered to the zygote in mice 52 . Our study also confirmed this finding. The full subset of sperm-derived RNAs will provide a useful reference to identify fertility-related interactions between paternally-derived RNAs. For instance, this study revealed that a binding partner of CABYR called A-kinase anchoring protein 3 (AKAP3) was also sperm-derived. CABYR, AKAP3, and Ropporin form a complex in the fibrous sheath of sperm, which has been implied to function as signaling for capacitation 53 . Additional research should be carried out to determine whether these RNAs produce proteins that form functional complexes during early embryonic development.
Ultimately, these data also allowed additional insight into the minor EGA of the 2-4 cell embryo. Minor EGA is a pre-emptive wave of transcriptional activity which occurs at the 1 or 2-4 cell stage of bovine embryo development 54 . This study uniquely compiles information from both the maternally and paternally contributed RNAs. Therefore, we identified a novel subset of RNAs that are exclusively transcribed by the 2-4 cell bovine embryo during minor EGA. Conversely, a subset of the gene transcripts we identified as paternally-derived in this study is actively transcribed the 2-4 cell embryo. For example, we classified cyclin A2 (CCNA2) and cyclin dependent kinase 2 (CDK2) as maternally-derived RNAs in our study. These genes cooperate during activation of the embryonic genome 55 . In another study, Graf et al. 49 found that CCNA2 and CDK2 belong to the category of RNAs that are maternally-derived but also transcribed by the embryo. However, several RNAs that were identified as embryo-specific through this study were also found to be present at the 4-cell stage in the study by Graf www.nature.com/scientificreports www.nature.com/scientificreports/ et al. 49 . These were deoxyribonuclease 1 (DNASE1), membrane associated ring-CH-type finger 3 (MARCH3), and ras-related C3 botulinum toxin substrate 3 (RAC3). DNASE1 can regulate DNA turnover and affects the stability of promoters 56 . MARCH3 mediates protein ubiquitination of Fc gamma receptor (FcγR), which affects responses to antibody-coated tumor cells 57 . RAC3 is a GTPase that alters cell growth by promoting cell adhesion and spreading 58 . These genes represent typical functions for early EGA RNAs transcribed by the embryo. However, this study also uncovered less traditional functions in the novel subset of early EGA genes.
Many of the unique embryo-specific RNAs in this study have roles in cell-cell signaling. These RNAs may serve as important regulators of embryo-mother communication during preimplantation development. For example, this study identified that sex hormone binding globulin (SHBG) expression was embryo-specific. SHBG is an androgen and estrogen transport protein which regulates plasma concentration of steroid hormones 59 . A deleterious mutation in SHBG is associated with prenatal death in dairy cattle 60 . This protein can be internalized by neurons or prostate cancer cells 61,62 . Further, SHBG can transport sex steroids into these target cells 61,62 . Another transcript found only in 2-4 cell embryos was tumor necrosis factor (TNF), a known mediator of embryo-maternal communication in mammals 63 . This gene modulates prostaglandin F2alpha (PGF2α) and is thought to aid in transferring the embryo from the oviduct to the uterus 64 . Additionally, secreted seminal-vesicle Ly-6 protein 1 (SOLD1) expression was also found to be embryo-specific in this study. SOLD1 is a signal at the fetomaternal interface that can regulate trophoblast invasiveness 65 . Future studies should interrogate whether these 2-4 cell embryo-specific RNAs can regulate signaling in the reproductive tract.
In conclusion, the present study demonstrated that sperm assists in regulating a substantial amount of pre-EGA RNA content in the early embryo and that male fertility is associated with RNA patterns in pre-EGA embryos. The knockdown of one of these genes, ATXN2L, resulted in substantially increased blastocyst development. This study also provided a whole-transcriptome characterization for oocyte-and sperm-contributed RNAs prior to EGA. A novel group of RNAs that are produced specifically by the 2-4 cell embryo was also identified. Future studies should evaluate whether non-RNA components contributed by the sperm to the zygote are capable of regulating pre-EGA RNA content in the embryo. Additionally, the roles of paternally-influenced RNAs in pre-EGA embryos should be further explored.

Methods ethics statement.
This study did not require approval from the Animal Care and Use Committee. Cows used for oocyte aspirations were not cared for at the University of Wisconsin -Madison facilities. Ovaries used for oocyte aspirations were purchased from Applied Reproductive Technology, LLC (Monona, WI, USA), and we were permitted by the company to perform in vitro fertilization (IVF) using the ovaries.
Bull semen selection. Semen samples were donated by Semex, Canada and the fertility status of sires was designated by the company as either high or low fertility based on the Repromax ™ system, which simultaneously accounts for Sire Conception Rate (SCR), Agri-Tech Analysis (ATA), and Canada's Non-Return Rate (NRR) data. Permission was granted by Semex to perform IVF. A total of 12 bulls were used for embryo production, which included six high fertility and six low fertility bulls.
In vitro embryo production. In vitro production of embryos was performed as described previously 66,67 .
Ovaries were supplied by Applied Reproductive Technology, LLC (Monona, WI, USA), from which follicles were aspirated for oocyte collection. Oocytes were washed in a TL-Hepes solution which contained 3% bovine serum albumin, sodium pyruvate, and gentamicin. Cohorts of 10 oocytes were matured for 24 h in 50 µl of M-199 maturation media which contained gonadotropins (FSH and LH), estradiol, sodium pyruvate, 10% fetal bovine serum, and gentamicin. In each of the six IVF replicates, matured oocytes were washed in the supplemented TL-Hepes solution and then split evenly into two groups for fertilization (200-380 oocytes/sire). One oocyte group was fertilized using a high fertility sire, and the second group was fertilized using a low fertility sire (Fig. 3). Overall, a total of 1737 oocytes were allocated to high fertility bulls and 1606 oocytes were allocated to low fertility bulls across six IVF replicates. New bulls were used in each replicate. In total, six high fertility bulls and six low fertility bulls were used for 2-4 cell embryo generation. Age differences of the bulls paired in each IVF replicate were checked using a paired t-test, and bull pairs used for IVF showed no substantial differences in age. Cohorts of 10 oocytes were transferred to 44 µl drops of fertilization medium, which was supplemented with fatty acid-free bovine serum albumin (FAF-BSA), sodium pyruvate, and gentamicin. Sperm preparation was performed with a Percoll gradient as previously described 68 . The final sperm concentration was adjusted to 1 million per mL, using 2 µl of sperm per drop. Penicillamine-hypotaurine-epinephrine and heparin (2 µl each) were added to fertilization droplets. Oocytes and sperm were co-cultured for 20 h. Then, presumptive zygotes were stripped of cumulus cells, washed with the supplemented TL-Hepes solution, and placed in 50 µl drops of CR1aa culture media 69,70 which were supplemented with FAF-BSA, sodium pyruvate, amino acids, and gentamicin. Presumptive zygotes were cultured in cohorts of 20-25. At 40 h post-fertilization, embryos that had reached the 2-4 cell stage were counted and collected into 100 µl of lysis solution for RNA-Sequencing (RNA-Seq). Statistical analysis of cleavage rates was performed using a paired t-test. Three biological replicates of embryo groups containing 37-70 cleavage-stage embryos for each high and low fertility bull were used for RNA-Seq. www.nature.com/scientificreports www.nature.com/scientificreports/ USA). Then, 150 bp paired-end sequencing was carried out on the HiSeq X Ten (Illumina). The quality of fastq files containing the raw sequence reads was evaluated using FastQC software (v 0.11.8). Adapters were trimmed from all sample fastq files using Trimmomatic software (v. 0.38) and the reads were simultaneously quality trimmed by removing reads shorter than 35 bp. All the trimmed files were re-checked for quality using FastQC software. Paired-end read files were aligned to the UMD3.1 bovine genome assembly (Bos taurusUMD3.1.dna. toplevel.fa) using Tophat (v 2.1.1) 71 . The number of reads per gene in each sample was calculated by HTSeq-count (v 0.6.1) 72 using UMD3.1 as reference genome (Bos_taurus.UMD3.1.94.gtf). Differential expression analysis was performed on 2-4 cell samples using EdgeR (v 3.20.7). The statistical tests were corrected for multiple testing using the Benjamini-Hochberg method as implemented in EdgeR 73 . A false discovery rate (FDR) cutoff of 10% was used to identify significant genes.

Evaluation of gene expression in biological replicates. Three additional biological replicates of in
vitro-produced 2-4 cell embryos were generated for validation purposes. New high-and low-fertility sires were used for each replicate. RNA was extracted from groups of 180-280 2-4 cell embryos for each sire. Whole transcriptome amplification was performed as described above. The iTaq Universal SYBR Green Supermix (Bio-Rad) was used for quantitative real-time PCR (qRT-PCR). Cycling was performed with a CFX-Connect Real-Time PCR Detection System (Bio-Rad), under the following conditions: 95 °C for 30 seconds, followed by 40 cycles of 95 °C for 5 seconds, and 60 °C for 30 seconds. To select a reference gene, the genes beta-actin (ACTB), ATP synthase F1 subunit beta (ATP5B), beta-2-microglobulin (B2M), and calnexin (CANX) were compared for expression stability across 2-4 cell embryo samples. The most stable gene, ACTB, was selected as the internal control 74 . All primers were designed to span an intron, in order to avoid DNA amplification (Supplementary Table S4). Fold change in gene expression was determined using the 2 −ΔΔ CT method as described 75 . If a sample CT was greater than 33, it was considered beyond the threshold for quantification. Validation within individual biological replicates was considered, and those replicates with >1.5-fold difference in accordance with sequencing results were categorized as individually validated. For determination of significance across all three replicates, a paired t-test of normalized gene expression values (ΔCT) was used.
Knockdown of ATXN2L using gapmer supplementation. ATXN2L was selected for further functional analysis based on its significance in validation, and because its role in fertility had not yet been explored. Knockdown was performed with a locked nucleic acid gapmer specific to ATXN2L, which was custom designed by Qiagen. Both 1 µM and 3 µM concentrations of the gapmer were evaluated, and the optimal concentration of 1 µM was selected. Each gapmer was supplemented to culture media 24 hours post-fertilization. Simultaneously, a control group of embryos was cultured, where an equivalent volume of water (vehicle of gapmer) was supplemented. Embryos remained in culture until blastocysts were graded at eight days post-fertilization. For grading, embryos were categorized as either blastocyst or degenerate. The blastocyst rate was calculated as the percentage of cleaved embryos that developed to the blastocyst stage. Three IVF replicates of ATXN2L gapmer knockdowns were carried out. Additionally, a cohort of embryos was supplemented with a bovine negative control gapmer, which is a scramble sequence designed so that it does not align to any known bovine sequences. Immediately after blastocyst grading, embryos were collected in lysis buffer, and RNA was extracted with the same method described above for the 2-4 cell embryos. Gene expression analysis was performed in order to confirm gapmer activity, and fold change in gene expression was calculated using the 2 −ΔΔCT method 75 . The significance of both blastocyst rate and gene expression was determined using a paired t-test.

Identification of RNAs expressed in sperm.
Sperm RNA was extracted from the two sires (one high fertility, one low) from one of the sequenced IVF replicates. One straw of cryopreserved semen per sire was www.nature.com/scientificreports www.nature.com/scientificreports/ subjected to RNA extraction. Each straw was thawed for one minute at 35-37 °C. Thawed semen was transferred to a 1.5 mL tube and centrifuged for four minutes at 4000 RPM. The supernatant was removed, and cells were suspended in a somatic cell lysis buffer 76 for four minutes on ice. Sperm was then centrifuged four minutes at 4000 RPM and lysis supernatant was aspirated. Somatic cell lysis buffer suspension was repeated once. Following somatic cell lysis, sperm were suspended in 0.75 mL of TRIzol (Invitrogen, Carlsbad, CA, USA), and RNA extraction was carried out according to manufacturer instructions. Then, RNA was amplified through WTA, purified, and sequenced. Reads were processed as described for 2-4 cell embryo samples. RNAs with a read count greater than five in either sample were considered to be expressed in sperm.

Identification of RNAs expressed in oocytes.
Studies that evaluated bovine oocyte RNA expression in multiple replicates were identified using NCBI Gene Expression Omnibus (GEO) (http://www.ncbi.nlm. nih.gov/geo). Two studies provided freely available data on bovine oocyte RNA expression in multiple replicates 77,78 , and were thus used for comparison to our sequencing results. Data from Xie et al. 77 is under accession number GSE18290 and data from Graf et al. 49 is under accession number GSE52415. The averages across oocyte replicates was used for comparison. The data from Xie et al. 77 was generated using the Affymetrix GeneChip ® Bovine Genome Array, which includes all high quality and publicly identified bovine transcripts. The Affymetrix-gene-code was converted to gene ID through DAVID 79,80 , and a minimum probe threshold of 40 raw intensity units 81 was used as a cutoff to classify whether or not a gene was expressed. The data from Graf et al. 49 were generated through RNA-Seq and aligned to UMD 3.1, and a cutoff of greater than five reads 78 was used to classify whether a gene was expressed in oocytes for this dataset. characterization of RnA origin. To classify the parental origin of RNAs, the reads from 2-4 cell embryo samples were merged with sequencing results from sperm samples and literature-derived oocyte RNA expression data. Any gene with greater than five reads in at least one of the six replicates of 2-4 cell embryos was considered expressed, and the thresholds stated above were used to determine sperm and oocyte expression. Genes that did not meet the threshold for any of the 2-4 cell embryo samples were omitted from the analysis. The parental origin of RNAs was determined based on RNA expression in different sample types. Ultimately, the RNAs were sorted into four categories: (1) both oocyte-and sperm-derived (2) oocyte-derived; (3) sperm-derived; (4) only 2-4 cell embryo. If a gene was only found in one gamete-specific dataset (rather than both) it was labeled as provisional. Provisionally-labeled genes were not included in final counts of parental origin, but all genes and parental origin labels can be found in Supplementary Table S3.
Accession codes. Oocyte RNA expression data from Xie et al. 77 is under accession number GSE18290, and oocyte RNA expression data from Graf et al. 49 is under accession number GSE52415.

Data availability
The majority of data generated is published in the current study (and its Supplementary Information files). The datasets generated during the current study are available from the corresponding author on reasonable request.