Next Article in Journal
Differences in Charge Distribution in Leishmania tarentolae Leishmanolysin Result in a Reduced Enzymatic Activity
Next Article in Special Issue
Age-Dependent Variations in Functional Quality and Proteomic Characteristics of Canine (Canis lupus familiaris) Epididymal Spermatozoa
Previous Article in Journal
Endothelium-Dependent Induction of Vasculogenic Mimicry in Human Triple-Negative Breast Cancer Cells Is Inhibited by Calcitriol and Curcumin
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
IJMS
Volume 23
Issue 14
10.3390/ijms23147642
ijms-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

In Silico Identification of lncRNAs Regulating Sperm Motility in the Turkey (Meleagris gallopavo L.)

by
Jan Pawel Jastrzebski
1,*,
Aleksandra Lipka
2,
Marta Majewska
3,
Karol G. Makowczenko
4,
Lukasz Paukszto
1,5,
Joanna Bukowska
6,
Slawomir Dorocki
7,
Krzysztof Kozlowski
8 and
Mariola Slowinska
9
1
Department of Plants Physiology, Genetics and Biotechnology, Faculty of Biology and Biotechnology, University of Warmia and Mazury in Olsztyn, 10-719 Olsztyn, Poland
2
Department of Gynecology and Obstetrics, School of Medicine, Collegium Medicum, University of Warmia and Mazury in Olsztyn, 10-561 Olsztyn, Poland
3
Department of Human Physiology and Pathophysiology, School of Medicine, Collegium Medicum, University of Warmia and Mazury in Olsztyn, 10-561 Olsztyn, Poland
4
Department of Animal Anatomy and Physiology, Faculty of Biology and Biotechnology, University of Warmia and Mazury in Olsztyn, 10-719 Olsztyn, Poland
5
Department of Botany and Nature Protection, Faculty of Biology and Biotechnology, University of Warmia and Mazury in Olsztyn, 10-719 Olsztyn, Poland
6
In Vitro and Cell Biotechnology Laboratory, Institute of Animal Reproduction and Food Research, Polish Academy of Sciences in Olsztyn, 10-748 Olsztyn, Poland
7
Department of Entrepreneurship and Spatial Management, Institute of Geography, Pedagogical University of Cracow, 30-084 Kraków, Poland
8
Department of Poultry Science, Faculty of Animal Bioengineering, University of Warmia and Mazury in Olsztyn, Oczapowskiego 5, 10-719 Olsztyn, Poland
9
Gamete and Embryo Biology, Institute of Animal Reproduction and Food Research, Polish Academy of Sciences in Olsztyn, 10-748 Olsztyn, Poland
*
Author to whom correspondence should be addressed.
Int. J. Mol. Sci. 2022, 23(14), 7642; https://doi.org/10.3390/ijms23147642
Submission received: 21 June 2022 / Revised: 6 July 2022 / Accepted: 9 July 2022 / Published: 11 July 2022
(This article belongs to the Special Issue Omics Strategies in Male Reproduction)
Browse Figures
Review Reports Versions Notes

Abstract

:
Long non-coding RNAs (lncRNAs) are transcripts not translated into proteins with a length of more than 200 bp. LncRNAs are considered an important factor in the regulation of countless biological processes, mainly through the regulation of gene expression and interactions with proteins. However, the detailed mechanism of interaction as well as functions of lncRNAs are still unclear and therefore constitute a serious research challenge. In this study, for the first time, potential mechanisms of lncRNA regulation of processes related to sperm motility in turkey were investigated and described. Customized bioinformatics analysis was used to detect and identify lncRNAs, and their correlations with differentially expressed genes and proteins were also investigated. Results revealed the expression of 863 new/unknown lncRNAs in ductus deferens, testes and epididymis of turkeys. Moreover, potential relationships of the lncRNAs with the coding mRNAs and their products were identified in turkey reproductive tissues. The results obtained from the OMICS study may be useful in describing and characterizing the way that lncRNAs regulate genes and proteins as well as signaling pathways related to sperm motility.

1. Introduction

The formation and maturation of the sperm in testis, known as spermatogenesis, is a series of complex molecular processes and morphological changes, including mitotic and meiotic divisions and differentiation leading from germ cells to mature spermatozoa [1,2]. These processes are under multiple regulators, such as hormones [3], genetic factors [4,5], long non-coding RNA (lncRNA) and epigenetics [6,7].
The molecular mechanism of sperm production from the primordial germ cells is well-documented in mammals, but there are only a few studies detailing this process in birds [5,8]. The final product of spermatogenesis is a spermatozoon that consists of the head, which is the carrier of the male genetic material, and the flagellum, which provides the motive force to move towards the oocyte. Spermatogenesis consists of several stages, including spermatocytogenesis, mitosis, meiosis and spermiogenesis [9]. Spermiogenesis is the transformation of spermatids into spermatozoa without further cell division. It includes the formation of an acrosome and an axoneme, loss of cytoplasm, and the replacement of nucleohistones with nucleoprotamine, which accompanies nuclei condensation [10,11]. In avians, spermiogenesis varies from species to species, e.g., in guineafowls, spermiogenesis involves ten distinct morphological steps [12]. However, in red junglefowl, spermiogenesis entails 8–10 steps [13,14,15] and 12 stages in Japanese quail [16,17].
The only function of the sperm is to fertilize the egg. For this, the sperm must be highly mobile, especially in the female reproductive system. Since sperm is transcriptionally inactive, mobility depends on flagellum and actin element formation in testes and post-testicular development of sperm motility [18]. Mobility dysfunction, which may be impaired at the molecular level, handicaps sperm and leads to increased infertility up to the complete loss of fertilization possibility [19,20].
The turkey (Meleagris gallopavo) is one of the most important livestock animals and the second largest contributor to the world’s poultry meat production [21]. The industrial turkey breeding is based on the artificial insemination because natural mating results in adequate fertility levels [22]. In comparison to other Galliformes species, turkey semen is characterized by some unique features, such as oxidative respiration as the main pathway in spermatozoa energy metabolism [23,24,25], distinctive proteome pattern of seminal plasma [26,27] and low semen sensitivity for cooling/freezing procedures [28]. The molecular mechanisms and key factors responsible for producing spermatozoa with good quality are still poorly understand for turkey spermatogenesis. So far, only the pathways associated with specificity of the reproductive tract have been described for turkey [5]. There is a lack of information describing molecular mechanism of spermatogenesis such as lncRNAs that can be involved in regulation of mature spermatozoa quality.
Semen quality is of economic importance and is studied to optimize reproductive performance [29] and scientific significance [30]. Investigating the genetic and epigenetic mechanisms underlying spermatogenesis allows us to better understand the entire process of sperm production and maturation on the molecular level. Sperm motility is the basic parameter of semen quality [4,31,32,33] in both mammals and birds [34], including turkeys [35].
Semen quality, as well as sperm motility, is influenced by gene expression and regulatory factors such as lncRNAs [7]. LncRNA is the most abundant group of transcripts apart from coding sequences and may be crucial in regulating gene expression and protein activity in maintaining basic cell functions as well as spermatogenesis [6,36]. In the current research, we focused on the identification of lncRNAs in the whole transcriptome of three turkey tissues: testis (T), epididymis (E), and ductus deferens (DD). Moreover, the deep in silico analyses revealed the potential interactions between differentially expressed lncRNAs (DELs) and differentially expressed genes (DEGs) engaged in the sperm motility within the male reproductive tract tissues.

2. Results

2.1. Sequencing, Mapping and Expression

As the result of sequencing, more than 100 million reads were obtained in each sample (Table 1). The highest number of raw reads was obtained for one of the individual T samples (165.8 raw reads), and the lowest was for the DD sample (102.6 of raw reads). The highest value of sequencing depth for both individual samples and mean was in the epididymis samples (176.9 million of reads in individual samples, mean = 152 million of reads). The trimming process reduced the quantity by about 15%, resulting in mean values of 83.22%, 83.47% and 84.33% of processed reads in ductus deferens, epididymis and testes, respectively. The mapping efficiency ranged from 96.6% to 98.9% with means of more than 98% in each group of samples. The count calling process resulted in the expression values for 56,037 genes and 103,847 transcripts in all samples. The number of expressed genes in individual samples ranged from 17,596 (“12_E”) to 43,510 (“14_T”), while the number of expressed genes in three analyzed tissues was 53,385 in testes, 35,313 in ductus deferens and 41,029 in epididymis. The number of expressed transcripts ranged from 22,109 (“12_E”) to 73,569 (“14_T”), and in summary for each tissue: testes—94,726, ductus deferens—64,905 and epididymis—75,765.

2.2. lncRNA Identification

The set of all pooled transcripts (103,847) from all samples was filtered following the pipeline presented in Figure 1. The whole assembled transcriptome contained 55,168 sequences of more than one exon, 86,326 sequences longer than 200 nucleotides, 87,374 non-protein coding and 100,248 transcripts with an expression of more than ten counts (Figure 2A).
The dataset of filtered sequences was then analyzed for coding potential using seven different tools. The results of various methods were different. The highest number of non-coding sequences was obtained from CPC2—20,344, PLEK—20,168 and CNCI—17,931, whereas FEELnc, CPAT, LncFinder and MDeep resulted in 1262, 1085, 922 and 2582 non-coding sequences, respectively (Figure 2B). Due to such differences, it was decided to apply the consensus that, if a given sequence, was confirmed by at least five methods (out of seven), it was qualified as non-coding. As a final result, 4391 lncRNA transcripts were obtained, where 863 were predicted and 3487 were annotated using the reference turkey ENSEMBL genome (version 5.01.104).

2.3. Differentially Expressed lncRNAs

The final lncRNA dataset was used to identify the differentially expressed lncRNAs (DELs) in each comparison. It resulted in 139 DELs in “epididymis/testis” comparison (93 up- and 46 downregulated), 182 DELs in “ductus deferens/testis” (132 up- and 50 downregulated) and 25 only upregulated DELs in “epididymis/ductus deferens” comparison (Figure 3, Table 2 and Table S1).
Potential interactions between all identified DEGs and DELs within each comparison were analyzed (Tables S2–S4), but only relations with DEGs involved in sperm motility processes, described in detail by Słowińska et al. [5], were further visualized, described and discussed (Table 3, Figure 4). Several CIS-relations were detected between DEGs and DELs (Table 4), but no such relation was found for genes involved in sperm motility. This is why these results were presented as supplementary data and are not discussed in the manuscript.
Ijms 23 07642 g004 550
Figure 4. The expression profiles of DELs in three comparisons based on trans-correlated DEGs: (A) DEGs of statistically significant expression difference between ductus deferens (DD) and other tissues-testis (T) and epididymis (E), (B) epididymis vs. testis and ductus deferens, (C) testis vs. epididymis and ductus deferens. The red triangle indicates the tissue against which the expression profiles of the other two tissues (green circles) were compared. The expression profiles of DEGs are presented as grey lines, blue lines are the profiles of positively correlated DELs and red are negatively correlated lncRNAs. Expression profiles (the Z-score values of FPKM-Fragments Per Kilo base of transcript per Million mapped fragments) of all identified DELs are presented as the heatmap in Figure 5.
Figure 4. The expression profiles of DELs in three comparisons based on trans-correlated DEGs: (A) DEGs of statistically significant expression difference between ductus deferens (DD) and other tissues-testis (T) and epididymis (E), (B) epididymis vs. testis and ductus deferens, (C) testis vs. epididymis and ductus deferens. The red triangle indicates the tissue against which the expression profiles of the other two tissues (green circles) were compared. The expression profiles of DEGs are presented as grey lines, blue lines are the profiles of positively correlated DELs and red are negatively correlated lncRNAs. Expression profiles (the Z-score values of FPKM-Fragments Per Kilo base of transcript per Million mapped fragments) of all identified DELs are presented as the heatmap in Figure 5.
Ijms 23 07642 g004
Ijms 23 07642 g005 550
Figure 5. The expression heatmap (the Z-score of the FPKM value) of the identified lncRNAs in each sample of each tissue grouped in the clusters (the dendrogram inside). The outer bars’ length corresponds to the chromosomal distribution of the DELs.
Figure 5. The expression heatmap (the Z-score of the FPKM value) of the identified lncRNAs in each sample of each tissue grouped in the clusters (the dendrogram inside). The outer bars’ length corresponds to the chromosomal distribution of the DELs.
Ijms 23 07642 g005

2.4. TRANS-Acting Genes

Analyses of the correlation of expression profiles resulted in 33,315 unique relations between DEGs and DELs. By limiting the results to the relationships between DELs and genes associated with the processes responsible for sperm motility 519 (32 DELs with 78 DEGs), 749 (72 DELs with 85 DEGs) and 99 (7 DELs with 28 DEGs) correlations were detected, respectively, in “epididymis/testis”, “ductus deferens/testis” and “epididymis/ductus deferens” comparison (Table 4). Detailed information about the types of relationships between individual protein-coding genes and DELs is presented in Tables S2–S4, and the summary of these relations for genes associated with sperm motility is available in Table S5.
The number of trans-correlations between DELs and DEGs involved in the DD/(T + E) group of sperm motility processes (the processes enriched for genes of significantly different expression between ductus deferens and both other tissues) was 229 (blue links—positive correlations and red—negative correlations in Figure S1A). Fourteen DEGs (the middle column in Figure 6 and grey lines in Figure 4A) were associated with six processes (microtubule cytoskeleton organization, microtubule-based movement, flagellated sperm motility, axonemal dynein complex assembly, outer and inner dynein arm assembly) showing statistically significant co-expression with 70 DELs (blue and red lines in Figure 4A and the third column in Figure S1A). Eleven genes of significantly different expression between the epididymis and other two tissues (group E/(T + DD)), associated with three processes (calcium ion binding, lipid metabolism and reproductive process), have significant co-expression with 29 DELs (Figure 4B and Figure 6A) by 88 relations. Genes that differ in expression between testis and other analyzed tissues (55 genes) were engaged in 12 processes related to sperm motility (six actin building processes, protein phosphorylation, apoptosis, response to stress, estrogen and oxygen-containing compounds and cation symporter activity) and showed 644 co-expression relations with 48 DELs (Figure 4C and Figure S1A). Almost all genes within each process showed a high and statistically significant correlation of expression with DELs. There were only a few genes without such co-expression. This is why the knowledge of each relationship should be extensively studied and enriched by additional analyses.

2.5. Direct Interactions

A total of 3101 RNA-RNA interactions were identified between 124 transcripts (50 DEGs) and 188 lncRNAs in the T/(E + DD) group (Figure S1B). Moreover, 796 potential interactions between 25 transcripts of 10 DEGs and 173 lncRNAs in the group E/(T + DD) (Figure 6B) were discovered. A total of 607 unique interactions between 26 transcripts of 14 DEGs and 163 lncRNAs were distinguished in the DD/(E + T) group (Figure S2B). The direct RNA–protein interactions’ analyses resulted in 5951 highly probable relations between 42 proteins (20 DEGs) and 338 lncRNAs in the T/(E + DD) group, 1089 results for 255 lncRNAs and nine targets (4 DEG) lncRNAs in the group E/(T + DD) and 757 interactions between 154 lncRNAs and eight protein sequences (3 DEGs). All detailed data on the interactions are stored in Table S5. The genes that showed the potential for RNA-RNA interactions (and co-expressing) but did not interact in protein type-RNA with DELs in the group DD/(T + E) were: DNAH7 and ZMYND10; in group E/(T + DD), these genes included DNAH5, CDHR2, CDH20; and, in the T/(E + DD) group, these genes were: SDCBP, CTB1, FAT1, TWF1, ACVR1, ARPC5, NEDD9, MET, SORBS1, ATP1A1, KRT19. The genes that showed the potential for RNA-protein but not RNA–RNA interactions in the DD/(T + E) group were: PLK4 and LRRC6 in group E/(T + DD) are CLGN, and in the group T/(E + DD) were: PACSIN2 and POF1B. The CCDC40 gene was detected in all analyzed types of correlations and interactions in the group DD/(T + E), LRP2 in group E/(T + DD), and the genes that showed all co-expression, RNA–RNA and RNA–protein interactions were: SLC9A3R1, SLIT2, GATA3, NET1, GSN, SPRY2, PDLIM3, BAG3, EVL, APP and EZR. The analyses of functional annotations using STRING resulted in 27 nodes and 21 edges with a p-value equal to 1.41 × 10−5 (see Figure 7).

2.6. Real-Time PCR-Validation of the Expression of Selected lncRNAs

The real-time PCR analysis, in all cases, confirmed the NGS expression levels of selected lncRNAs in turkey testis, epididymis and ductus deferens (Figure 8).

3. Discussion

In this study, lncRNAs were identified in the transcriptome of three turkey tissues, i.e., testes, epididymis and ductus deferens. In particular, the study focused on the analysis of lncRNAs, which can regulate the activity of genes and proteins involved in sperm motility. The present studies are a continuation of the Slowińska et al. [5] study. However, the current studies are the first attempt to identify and analyze the differences in the expression of lncRNAs that may have the potential to regulate sperm motility.
Among all 4391 identified lncRNAs, 139, 182 and 25 DELs were discovered in the following cases: “epididymis/testis”, “ductus deferens/testis” and “epididymis/ductus deferens”, respectively. The study analyzed the distance relations (CIS), co-expression (TRANS), RNA–RNA and protein–RNA direct interactions between DELs and DEGs for each comparison. Based on these results, lncRNAs were indicated which may regulate processes related to sperm motility. The results obtained from the OMICS study may be useful to describe and characterize the way that lncRNAs regulate genes and proteins as well as signaling pathways related to sperm motility.

3.1. The Regulations between Testes and Post-Testicular Areas of Spermatogenesis

Testes (Group T/(E + DD))

The evidence indicates that testicular and post-testicular development of sperm varies substantially [37]. Consistent with previous results [5], the differences in the expression of DELs between studied tissues were observed. The abundance and diversity of the revealed DELs with strong and statistically significant correlations as well as potential interactions suggest that lncRNAs may play an important role in the genetic background by determining differences in processes associated with sperm maturation and motility, as well as the entire spermatogenesis. The genes that had significant differences in expression between the testis and epididymis and ductus deferens revealed the potential for interaction with lncRNA and thus the ability to regulate spermatogenesis processes in turkey. The current study identified lncRNAs that may affect mentioned genes involved in the spermatogenesis, in particular those related to sperm motility. One such gene whose activity is differentially regulated between the testes and the ductus deferens that may be crucial in sperm motility is ARPC5 (short names of genes and macromolecules are explained in the Abbreviations section). ARPC5 takes part in actin polymerization, essential for sperm motility and capacitation processes. Moreover, it acts through the actin nucleation for the production of the actin filaments. The ARP2/3 complex consists of the two proteins ARP2 and ARP3 and five additional subunits ARPC1-5 [38,39]. ARP2 and ARP3 are the core of the mechanism that catalyzes the nucleation of actin fibers by ATP hydrolysis [40,41,42], while the remaining subunits stabilize the ARP2/3 complex [43]. The ARP2/3 complex plays a pivotal role in the regulation of sperm function and male fertility through actin polymerization [44]. The authors’ research of both co-expression analyses and direct RNA–RNA interactions revealed that the ARPC5 gene may be regulated by numerous lncRNAs. Both positive and negative co-expression was confirmed in the testis vs. epididymal (Table S3) and testis vs. ductus deferens comparisons (Table S2). However, in the case of direct interactions, a potential interplay was demonstrated only at the RNA-RNA level, and ARPC5 has a high probability of interacting with more than 60 DELs, and 18 of those DELs have the potential to interact only with the ARPC5 (Figure S1).
The current results identified two lncRNAs (MSTRG.13455 and MSTRG.34057) that had a negative correlation with the expression of ARPC5. This phenomenon indicated that the expression of MSTRG.13455 and MSTRG.34057 in the testicles was higher than in other tissues. In contrast, the direct interaction of those lncRNAs has been demonstrated with other genes related to sperm motility, i.e., SPTBN1, NEDD9, CGNL1 and MYH11.
SPTBN1 connects the cell membrane with the actin cytoskeleton and plays a role in determining the shape of the cells, transmembrane protein distribution and organization of organelles [45]. The current research indicates that the SPTBN1 protein can directly interact with a few identified lncRNAs (MSTRG.13455, MSTRG.34057, MSTRG.35249, MSTRG.45880). The expression of the aforementioned lncRNAs was elevated in the testes, but the RNA–protein interaction may take place in the area of further stages of spermatogenesis, which is suggested by the negative correlation of expression. It was also shown that MSTRG.13455, MSTRG.31149, MSTRG.35249 lncRNAs can directly interact with CGLN1 protein. The physiological expression profile of CGNL1 plays a key role in testes development through the differentiation and dynamics of the Sertoli cell cytoskeleton. Moreover, non-disturbed expression and full activity of the CGLN1 are required for the appropriate flow of spermatogenesis [46]. In turn, NEDD9 regulates actin dynamics through cortactin deacetylation [47].
In this research, it was hypothesized that LEF1 may be regulated by lncRNA on both RNA–RNA and RNA–protein levels. The LEF1 is involved in the apoptotic process and is an important member of the Wnt/β-catenin signaling pathway. Moreover, it was reported that LEF1 may contribute to the arrest of spermatogonial differentiation [48]. ANXA1 plays diversified roles in cell functioning. It affects cell proliferation and differentiation as well as plasma membrane repair and apoptosis [49]. The expression level of ANXA1 is related to semen DNA fragmentation, making it a possible new biomarker for sperm DNA quality [50]. Both NRP1 and its cooperation with VEGFA are critical for vascular development in various angiogenic processes. It is known that NRP1 enhances signal transduction of VEGFA isoforms. It is suggested that decreased expression of the NRP1 in Sertoli cells may inactivate VEGFA and thus ultimately hinder the establishment and proliferation of the spermatogonial stem cells (SSCs), leading to an imbalance between undifferentiated spermatogonia self-renewal, differentiation and impaired spermatogenesis [51]. The SSCs may also be regulated by ITGB1 as an indirect way of response to androgens that involves interaction with E-cadherin on SSCs to regulate its fates [52,53]. Furthermore, SORBS1 is associated with the proper function of microtubules and cytoskeleton structure. Moreover, it can also be involved in the formation of multinuclear round spermatids and the formation of mature sperm [54]. Stöckl et al. revealed that CNN3 regulates the actin cytoskeleton and is involved in cell fusion and myogenesis. However, its specific function in spermatogenesis is still unknown [55]. Nevertheless, it may be speculated that CNN3 may be engaged in cytoskeletal dynamics and the dramatic transformation that the SSCs undergo to ultimately produce highly specialized spermatozoa [56]. NET1 is also involved in the control of the cell cycle, and through regulation of rearrangements of actin, cytoskeleton can control cell motility [57,58]. The methylation and expression level of GATA3 can affect sperm motility and male fertility status [59]. MICAL2 covalently modifies and affects actin [60] and, due to its interactions with cas proteins, may regulate integrin-mediated cell adhesion and cytoskeletal organization [61].
Differentially expressed lncRNAs detected between analyzed turkey male reproductive tissues have the potential to regulate the aforementioned genes. Moreover, the current results indicated that lncRNAs may have the potential for multilevel regulation of actin and cytoskeleton modifications, as well as basic cell functions including proliferation and differentiation. Furthermore, these processes and functions are known to be crucial for appropriate sperm production; however, the current research gives an in-depth picture of how they may be regulated.

3.2. Post-Testicular Development

3.2.1. Epididymis (Group E/(T + DD))

In the epididymis, among all DELs, five lncRNAs were identified with significantly different expressions than in the other two analyzed tissues. Those lncRNAs showed a high potential for direct RNA–protein interactions with multiple proteins. The results indicate that MSTRG.45728, MSTRG.44914 and MSTRG.39415 may have the ability to directly regulate the activity of the PLA2G12B and CLGN proteins. Moreover, MSTRG.18092 is likely to interact with PLA2G12B, CLGN, HNF4A and LRP2. The HNF4A is a member of the nuclear transcription factors and contains a conserved nucleic-acid binding domain (Figure 9) [62].
Interestingly, HNF4A, CLGN and LRP2 genes can be regulated by lncRNAs at the level of RNA–RNA interactions. Moreover, another two lncRNAs were identified, i.e., MSTRG.45871 and MSTRG.47677, whose expression profiles revealed a high correlation (r > 0.9, p-value < 1 × 10−4) with the LRP2, as well as a high probability of direct interaction with its product. Possible direct interactions were also identified on the RNA–RNA level with genes such as DNAH5, HEXB, CDHR2, CDH20, HPGDS and RDH10. The DNAH5 gene product is involved in the assembly of the outer axon dynein arm of flagella or cilia, and a mutation within this gene contributes to a risk of male infertility [63,64]. It is highly possible that any changes in its expression may affect post-testicular sperm maturation. According to available data, HPGDS may play a key role in the regulation of the spermatogenesis mechanism [65]. On the other hand, the HEXB gene is involved in the binding of the sperm to zona pellucida [45,66], which is crucial for effective fertilization. However, there is no information on those genes’ function in the epididymis. However, taking into consideration that epididymis is important for avian sperm to achieve appropriate quality for fertilizing capacity and motility [5,67], it can be expected that presently identified lncRNAs with multivariate interactions on different levels may have huge potential to regulate the course of the spermatogenesis and sperm maturation.

3.2.2. Ductus Deferens (Group DD/(T + E))

Detailed analysis showed that the processes related to the formation of dynein arms can be regulated by the cross-talk between lncRNA with the following genes: ARMC4, CCDC40, DNAH7, LRRC6, PIH1D3, ZMYND10, DEUP1, NDC80, PLK4, RSPH9, TTC26, CFAP221, ROPN1L and CFAP69. DEGs associated with inner and outer dynein arm assembly in epididymis were also detected and described in the authors’ previous study [5]. The high value of the co-expression correlation coefficient between CCDC40 and lncRNA, which was identified both in the analyses of direct RNA–protein interactions and in the analysis of expression profiles: MSTRG.17922 (r = 0.97 p-value = 1.88 × 10−7) and CCDC40 and MSTRG.45871 (r = 0.92, p-value = 2.18 × 10−5), suggest that these two lncRNAs may play a supporting role in the formation of protein complexes of dynein arms (MSTRG.17922 log2(FC) = 2.35 and MSTRG.45871 log2(FC) = 2.65 in comparison DDE). The ARMC4 gene shows a high co-expression correlation (r = 0.94, p-value = 6.26 × 10−6) with MSTRG.1945 lncRNA. This lncRNA is highly overexpressed in the epididymis in relation to the ductus deferens (log2(FC) = 10.2). The DNAH7 reflects the potential for direct interaction and co-expression with MSTRG.47677 lncRNA, which is expressed in the epididymis but not identified in the ductus deferens. It has also been revealed that two lncRNAs characterized by high overexpression in the testes relative to the ductus deferens (MSTRG.34479: log2(FC) = 3.59 and MSTRG.40956: log2(FC) − 6.39) depicted very similar expression profile with DNAH7 (in both cases r = 1 and p-value = 9.46 × 10−7 and 7.45 × 10−9 respectively) and a high potential for direct RNA–RNA interaction. ARMC4 may have an evolutionarily conserved function in spermatogenesis by being involved in cilia/flagella assembly, structure or function [68,69], and DNAH7 is a component of the inner dynein arm of ciliary axonemes [70]. Interestingly, the CCDC40 plays a key role in cilia formation and function and is an element of the dynein regulatory complex [71]. Mutations within the CCDC40 are associated with cilia and sperm dysmotility [72], and are suspected to be involved in idiopathic male infertility [64]. Moreover, co-expression and interactions between lncRNAs and the CCDC40 and ZMYND10 are involved in bovine endometritis [73]. These distant examples demonstrate that those genes are prone to the regulatory function of the lncRNAs. Hence, the proper formation of the cilia may be one of the processes regulated by such an interaction.
All DELs detected in the current study were novel lncRNAs, not yet annotated; therefore, their potential function was discussed based on postulated functions of genes with which they interact or correlate. However, it should be mentioned that the impact of lncRNAs not only manifests through gene expression regulation. Besides transcriptional level regulation, at the post-transcriptional level, lncRNAs may be involved in mRNA processing, as well as direct interactions with proteins to regulate translation and post-translation modifications [74,75]. Within the authors’ research, numerous lncRNAs interactions as distance relations (CIS), co-expression (TRANS), RNA–RNA and protein–RNA direct interactions were detected. CIS-regulation can act by two main, independent mechanisms that do not exclude each other, the lncRNA can regulate neighboring loci (within the same chromosome), and/or can affect chromatin state or steric impediment that influences the expression of nearby genes [76]. Conversely, lncRNAs act as trans when they affect genes located on other chromosomes [36]. RNA–RNA interactions may also refer to interactions between two RNAs (inter-molecular) or between different regions of one RNA molecule (intra-molecular) and both mechanisms are significant for RNA function and regulation [77]. RNA–protein interactions act through the formation of dynamic ribonucleoprotein complexes [78]. Although the nature and dynamics of the aforementioned interactions still need to be investigated [79], nevertheless, multiple features of lncRNAs can determine their functionality. These features include sequence, expression level, processing, cellular localization, structural organization and interactions with other molecules. Moreover, the regulatory effect is mediated by the lncRNA transcript or by its transcription [75,79].
The current research focused on the comparison of gene and lncRNA expression in different reproductive tissues regardless of sperm quality. The study aimed to establish a general pattern of lncRNA expression in each part of the turkey’s reproductive tracts. Since spermatogenesis in birds is very diverse, it should be helpful to characterize this process in the turkey. The introduced in silico approach indicated target genes and regulatory elements that should be further investigated to provide mechanistic insight. It would be very interesting to check how identified DEGs and DELs are reflected in the proteome and how this correlates with semen quality parameters. Further research including single-cell sequencing should also be considered to distinguish cell clusters within reproductive tracts and their role during spermatogenesis [80]. The results are a step forward and potentially may indicate directions for further research.

4. Materials and Methods

4.1. Tissue Collection

The turkeys came from an experimental farm of the Poultry Department (University of Warmia and Mazury in Olsztyn, Poland). As described previously [5], male turkeys were photostimulated at 26 wks of age (14 h light to 10 h darkness) and produced semen by 30 wks of age. Semen was collected at one-week intervals, by abdominal massage [67]. The sperm concentration, viability, motility and volume were 9.5 ± 1.5 × 109 sperm/mL, 97.9 ± 0.3%, 85.2 ± 9.6% and 0.6 ± 0.05 mL, respectively. Male reproductive tract tissue samples were collected from six, 38-week-old turkeys that were slaughtered at a local slaughterhouse. Tissues were immediately frozen in liquid nitrogen and used for total RNA isolation. All experiments were conducted in accordance with the guidelines of the International Guiding Principles for Biomedical Research Involving Animals proposed by the Society for Reproductive Research and the Polish Act on Animal Welfare. These experiments did not require ethical consent. The sampling process was described in detail previously [5].

4.2. RNA Isolation, Library Preparation and Sequencing

All procedures from RNA isolation up to sequencing procedures were performed as described previously [81]. The key procedures in this process included the following steps. The total RNA was extracted using a Total RNA Mini kit, and the concentration was determined spectrometrically (NanoDrop-1000 Technologies LLC, Wilmington, DE, USA). The RNA integrity was checked using a 2100 Bioanalyzer, with an RNA 6000 Nano LabChip Kit and the samples with RIN above 8.0 and rRNA ratios above 1.0 were used for NGS. The final strand-specific RNA sequencing libraries were prepared by synthesizing the sequences with PCR and using the TruSeq Stranded mRNA Library Prep Kit (Illumina, San Diego, CA, USA) used according to the manufacturer’s protocol. The cDNA libraries were sequenced using the NovaSeq 6000 platform (Illumina, San Diego, CA, USA) as the paired-end, 2 × 100-bp sequences. Library preparation and RNA-Seq were outsourced to Macrogen (Macrogen, Inc., Seoul, Korea).

4.3. In Silico Preprocessing, Mapping and Expression Analyses

The raw reads, trimmed reads and mapping quality were checked using FastQC software, v. 0.11.7 (Bioinformatics Group at the Babraham Institute, Cambridge, UK; www.bioinformatics.babraham.ac.uk (accessed on 20 June 2022)). The trimming process was performed using the Trimmomatic software, v. 0.38 [82] and included adapter removing, low quality read filtering (Phred cutoff score ≤ 20; with 10-bp frameshifts) and reads trimmed to equal lengths = 90 nucleotides. Preprocessed reads were mapped to the turkey reference genome (version 2.01), using annotation version 2.01.91 and the transcriptome was then additionally annotated using genome version 5.01 with annotation version 5.01.104. The first reference was used because this research is a continuation of the analyses from the previous work [5], and it was crucial that the data could be compared directly and the 5.01.104 genome was used to analyze the release version with new lncRNA annotations.
The mapping was performed with the STAR software, v. 2.4 (Cold Spring Harbor Laboratory, Cold Spring Harbor, NY; https://github.com/alexdobin/STAR (accessed on 30 October 2021) [83] and StringTie, v. 1.3.3 (https://ccb.jhu.edu/software/stringtie (accessed on 8 October 2021)) [84] was used to enrich additional annotations by merging the reference GTF file with BAM files.

4.4. Expression Profiling and Differentially Expressed Genes

Expression profiles as the count matrices were enriched at both levels: transcripts and genes. The analysis of differential expression was performed on the gene level using the Ballgown [85] and Cuffdiff [86] methods. Genes were classified as the differentially expressed if the adjusted p-value was <0.01, the logarithmic value of expression fold-change was >1.0) and were confirmed by both methods. All statistical analyses were done in R Bioconductor packages [87].

4.5. The lncRNA Identification

The transcript sequences were extracted from the reference genome using the annotation file obtained in this research (merged gff files from all samples) with gffread script (https://github.com/gpertea/gffread (accessed on 30 October 2021)). The transcript count matrix was obtained simultaneously with the gene count matrix using merged GTF. It allowed control of the annotations on both gene and transcript levels. The reference turkey genome version 5.01.104 contains 1703 annotated lncRNA transcripts, and this genome was used (in the first step of the novel lncRNA identification process) to annotate (using BLAST similarity > 90%, E-value < 1 × 10−10) and extract known lncRNA transcripts from the transcriptomic dataset (Figure 1). The next steps of the filtering stage involved: removing all annotated protein coding transcripts, removing all transcripts of summarized expression lower than 10 counts, and filtering out all transcripts shorter than 200 nucleotides and consisted of only one exon. The filtered dataset was then analyzed for coding probability using seven methods: CPC2 [88], FEELnc [89], CPAT [90], PLEK [91], CNCI [92], LncFinder [93] and LncRNA Mdeep [94]. Only those sequences that were estimated as non-coding by at least five methods were finally filtered using the pFam [80] database and were classified as novel lncRNAs. Both datasets of novel and known lncRNAs were merged and used for the identification of differentially expressed lncRNAs (DELs) in each comparison. Consequently, lncRNAs were classified as DELs if they met the same conditions as DEGs (adjusted p-value < 0.01, log2(FC) > 1.0).

4.6. The lncRNA Interactions with Sperm Motility-Related DEGs

The merged three sets of DELs and DEGs were used to detect potential correlations and interactions. This work focused on the analysis of lncRNA relationships only with genes associated with sperm motility processes. The set of these genes is presented in Table 3 and has been described in detail in the authors’ previous manuscript [8]. These genes were classified into three groups based on their expression profiles. The first group DD/(T + E) consisted of processes enriched for the genes in which a significant difference of expression was detected between ductus deferens and other two tissues testis and epididymis. The group E/(T + DD) contains genes of differential expression in the epididymis in contrast to testis and ductus deferens, and the last group, T/(E + DD), contained the processes and genes significant in comparisons “testis/ductus deferens” and “testis/epididymis” (Figure 5). DEGs located on the genome closer than 100 kbp from DEL were classified as CIS-acting genes. TRANS-acting genes were those of similar expression profiles calculated as Pearson correlation and selected only those of r coefficient >0.9 and p-value < 0.05.
Direct interactions of DELs with potential target macromolecules were detected using two programs: LncTar v. 1.0 [95] and lncPro v. 1.0 [96] tools. LncTar calculated the binding strength of lncRNA to mRNA, and the results were considered significant when the ndG value was lower than –0.1 for each interaction. The lncPro program calculated the probability of hydrogen-bonding propensities and Van der Waal’s interaction between the secondary structures of the lncRNAs and proteins encoded by DEGs. If the probability was higher than 0.90, such interaction was considered significant.
The STRING [97] interaction network was constructed for proteins detected using the lncPro tool (lncRNA–protein interactions). The minimum required interaction score was set as 0.4, and all disconnected nodes were hidden. Most of the plots (such as Volcano and MA plots, Sankey plots of TRANS-acting relations, DELs heatmap) were prepared using the ggplot2 library in the R Bioconductor environment. The RNA–RNA interactions were visualized using the Cytoscape tool [98].
The whole process of novel lncRNA identification was performed with the homemade R library which is under preparation for publishing. All high-performance computations were performed using a 60-core CPU and a 136 GB RAM server at the Regional IT Center of the University of Warmia and Mazury in Olsztyn.

4.7. Real-Time PCR

As the validation of RNA-seq results, Real-time PCR was performed to detect the expression levels of MSTRG.17922, MSTRG.6075, MSTRG.8989 and MSTRG.25698. Total RNA samples used in the NGS experiment were transcribed to cDNA using the QuantiTect® Reverse Transcription Kit (Qiagen, Valencia, CA, USA), according to the manufacturer’s protocol. The expression of selected lncRNAs in testis, epididymis and ductus deferens were analyzed to validate obtained NGS data. Briefly, Real-time PCR reaction included the following: cDNA template, Power SybrGreen® Master Mix (Life Technologies, Grand Island, NY, USA), RNase-free water to a final volume of 20 μL and primers selected for validation of NGS data. Each Real-time PCR reaction was performed in duplicate. Non-template controls were used for each set of primers to confirm the specificity of the reaction. The relative amplification of selected lncRNAs was calculated using the ΔΔCt method according to the Guide to Performing Relative Quantitation of Gene Expression Using Real-Time Quantitative PCR protocol and normalized using the geometric mean of the expression levels of a high stability reference gene in turkey testis, epididymis and ductus deferens, i.e., glyceraldehyde 3-phosphate dehydrogenase (GAPDH). Statistical analyses and graphs of the Real-time PCR results were made using the ‘pcr’ package [99] in the R environment. Statistical testing of PCR data was performed using the linear regression model [100].

5. Conclusions

We identified 863 new lncRNAs in the ductus deferens, testes and epididymis of turkeys and for the first time we characterize their potential role in the regulation of sperm maturation processes, in particular sperm motility. This offers the opportunity to further investigation of the function and role of lncRNAs in processes of molecular regulation.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ijms23147642/s1.

Author Contributions

Conceptualization, J.P.J.; methodology, J.P.J., L.P., K.G.M. and M.S.; software, J.P.J.; validation, J.B. and S.D.; formal analysis, J.P.J., K.G.M. and L.P.; investigation, J.P.J., M.S. and J.B.; resources, J.P.J., K.K. and M.S.; data curation, J.P.J., K.G.M., L.P. and S.D.; writing—original draft preparation, J.P.J. and A.L.; writing—review and editing, J.P.J., A.L. and M.M.; visualization, J.P.J. and K.G.M.; supervision, J.P.J.; project administration, J.P.J.; funding acquisition, J.P.J. and M.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Science Centre, Poland, Grant No.: 2016/21/B/NZ9/03583.

Institutional Review Board Statement

Not applicable; all experiments were performed according to the Polish Animal Welfare Act. No ethics approval was required for these experiments.

Informed Consent Statement

Not applicable.

Data Availability Statement

The whole project is available in the BioProject database in the NCBI under accession number: PRJNA597008. The raw data are deposited in the SRA database under SRP238360 ID. The expression tables and all transcriptomic sequences in fasta format are available in GEO database (GSE142428). All sequences of predicted lncRNAs were extracted from transcriptome fasta file and submitted to NCBI GenBank under accession numbers ranging from OM983514 to OM984417.

Acknowledgments

The authors thank Stefano Pascarella (Sapienza Università di Roma) for useful discussions, technical and organizational support. Jan Paweł Jastrzębski was the recipient of an internship from the “Development Program of the University of Warmia and Mazury in Olsztyn”, Task 16. “Implementation of internship programs for academic and didactic staff of the UWM in Olsztyn” (POWR.03.05.00-00-Z310/17) funded by the European Social Fund, during which a part of the work described in this paper was conducted.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

ACVR1Activin receptor type-1
ANXA1Annexin A1
APPAmyloid-beta precursor protein
ARMC4Armadillo repeat containing 4
ARPC1-5Actin related protein 2/3 complex subunit from 1 to 5
ATP1A1ATPase Na+/K+ transporting subunit alpha 1
BAG3BCL2 associated athanogene 3
CCDC40Coiled-coil domain-containing protein 40
CDH20Cadherin 20
CDHR2Cadherin related family member 2
CFAP221Cilia and flagella-associated protein 221
CFAP69Cilia and flagella associated protein 69
CGNL1Cingulin Like 1
CLGNCalmegin
CNN3Calponin 3
CTBP1C-terminal binding protein 1
DEUP1Deuterosome assembly protein 1
DNAH5Dynein axonemal heavy chain 5
DNAH7Dynein axonemal heavy chain 7
EVLEna/vasodilator-stimulated phosphoprotein-like protein
EZREzrin
FAT1FAT atypical cadherin 1
GATA3GATA binding protein 3
GSNGelsolin
HEXBHexosaminidase subunit beta
HNF4AHepatocyte nuclear factor 4 alpha
HPGDSHematopoietic prostaglandin D synthase
ITGB1Integrin subunit beta 1
KRT19Keratin 19
LEF1Lymphoid enhancer binding factor 1
LRP2LDL receptor related protein 2
LRRC6Leucine rich repeat containing 6
METMET proto-oncogene. receptor tyrosine kinase
MICAL2Microtubule associated monooxygenase. calponin and LIM domain containing 2
MYH11Myosin Heavy Chain 11
NDC80NDC80. kinetochore complex component
NEDD9Neural precursor cell expressed, developmentally down-regulated 9
NET1Neuroepithelial transforming gene 1
NRP1Neuropilin 1
PACSIN2Protein kinase C and casein kinase substrate in neurons 2
PDLIM3PDZ and LIM domain 3
PIH1D3PIH1 domain containing 3
PLA2G12BPhospholipase A2 group XIIB
PLK4Polo like kinase 4
POF1BPOF1B. actin binding protein
RDH10Retinol dehydrogenase 10
ROPN1LRhophilin associated tail protein 1 like
RSPH9Radial spoke head 9 homolog
SDCBPSyndecan binding protein
SLC9A3R1SLC9A3 regulator 1
SLIT2Slit guidance ligand 2
SORBS1Sorbin and SH3 domain containing 1
SPRY2Sprouty RTK signaling antagonist 2
SPTBN1Spectrin beta, non-erythrocytic 1
TTC26Tetratricopeptide repeat domain 26
TWF1Twinfilin actin binding protein 1
VEGFAVascular endothelial growth factor A
ZMYND10Zinc finger MYND-type containing 10

References

  1. Pardyak, L.; Kaminska, A.; Brzoskwinia, M.; Hejmej, A.; Kotula-Balak, M.; Jankowski, J.; Ciereszko, A.; Bilinska, B. Differences in aromatase expression and steroid hormone concentrations in the reproductive tissues of male domestic turkeys (Meleagris gallopavo) with white and yellow semen. Br. Poult. Sci. 2018, 59, 591–603. [Google Scholar] [CrossRef] [PubMed]
  2. Pardyak, L.; Kaminska, A.; Brzoskwinia, M.; Hejmej, A.; Kotula-Balak, M.; Jankowski, J.; Ciereszko, A.; Bilinska, B. Differential expression of cell-cell junction proteins in the testis, epididymis, and ductus deferens of domestic turkeys (Meleagris gallopavo) with white and yellow semen. Poult. Sci. 2020, 99, 555–566. [Google Scholar] [CrossRef] [PubMed]
  3. Hess, R.A.; Sharpe, R.M.; Hinton, B.T. Estrogens and development of the rete testis, efferent ductules, epididymis and vas deferens. Differentiation 2021, 118, 41–71. [Google Scholar] [CrossRef] [PubMed]
  4. Fraser, L.; Brym, P.; Pareek, C.S.; Mogielnicka-Brzozowska, M.; Paukszto, Ł.; Jastrzębski, J.P.; Wasilewska-Sakowska, K.; Mańkowska, A.; Sobiech, P.; Żukowski, K. Transcriptome analysis of boar spermatozoa with different freezability using RNA-Seq. Theriogenology 2020, 142, 400–413. [Google Scholar] [CrossRef] [PubMed]
  5. Słowińska, M.; Paukszto, Ł.; Jastrzębski, J.P.; Bukowska, J.; Kozłowski, K.; Jankowski, J.; Ciereszko, A. Transcriptome analysis of turkey (Meleagris gallopavo) reproductive tract revealed key pathways regulating spermatogenesis and post-testicular sperm maturation. Poult. Sci. 2020, 99, 6094–6118. [Google Scholar] [CrossRef]
  6. Fraser, L.; Paukszto, Ł.; Mańkowska, A.; Brym, P.; Gilun, P.; Jastrzębski, J.P.; Pareek, C.S.; Kumar, D.; Pierzchała, M. Regulatory potential of long non-coding rnas (LncRNAs) in boar spermatozoa with good and poor freezability. Life 2020, 10, 300. [Google Scholar] [CrossRef]
  7. Joshi, M.; Rajender, S. Long non-coding RNAs (lncRNAs) in spermatogenesis and male infertility. Reprod. Biol. Endocrinol. 2020, 18, 103. [Google Scholar] [CrossRef]
  8. Harvey, S.; Scanes, C.G.; Phillips, J.G. Advances in Experimental Medicine and Biology. In Avian Reproduction; Sasanami, T., Ed.; Springer: Singapore, 2017; Volume 1001, ISBN 978-981-10-3974-4. [Google Scholar]
  9. Johnson, L.; Varner, D.D.; Roberts, M.E.; Smith, T.L.; Keillor, G.E.; Scrutchfield, W.L. Efficiency of spermatogenesis: A comparative approach. Anim. Reprod. Sci. 2000, 60–61, 471–480. [Google Scholar] [CrossRef]
  10. Oliva, R.; Mezquita, C. Marked differences in the ability of distinct protamines to disassemble nucleosomal core particles in vitro. Biochemistry 1986, 25, 6508–6511. [Google Scholar] [CrossRef]
  11. Sprando, R.L.; Russell, L.D. Spermiogenesis in the bluegill (Lepomis macrochirus): A study of cytoplasmic events including cell volume changes and cytoplasmic elimination. J. Morphol. 1988, 198, 165–177. [Google Scholar] [CrossRef]
  12. Aire, T.A.; Olowo-okorun, M.O.; Ayeni, J.S. The seminiferous epithelium in the guinea fowl (Numida meleagris). Cell Tissue Res. 1980, 205, 319–325. [Google Scholar] [CrossRef] [PubMed]
  13. De Reviers, M.; Richetin, C.; Brillard, J.P. LE DÉVELOPPEMENT TESTICULAIRE CHEZ LE COQ. II.—MORPHOLOGIE DE L’ÉPITHÉLIUM SÉMINIFÈRE ET ÉTABLISSEMENT DE LA SPERMATOGENÈSE. Ann. Biol. Anim. Biochim. Biophys. 1971, 11, 531–546. [Google Scholar] [CrossRef]
  14. Gunawardana, V.E. Stages of spermatids in the domestic fowl: A light microscope study using Araldite sections. J. Anat. 1977, 123, 351–360. [Google Scholar] [PubMed]
  15. Tiba, T.; Yoshida, K.; Miyake, M.; Tsuchiya, K.; Kita, I.; Tsubota, T. Regularities and Irregularities in the Structure of the Seminiferous Epithelium in the Domestic Fowl (Gallus domesticus). Anat. Histol. Embryol. 1993, 22, 241–253. [Google Scholar] [CrossRef]
  16. Lin, M.; Jones, R.C.; Blackshaw, A.W. The cycle of the seminiferous epithelium in the Japanese quail (Coturnix coturnix japonica) and estimation of its duration. Reproduction 1990, 88, 481–490. [Google Scholar] [CrossRef]
  17. Lin, M.; Jones, R.C. Spermiogenesis and spermiation in the Japanese quail (Coturnix coturnix japonica). J. Anat. 1993, 183 Pt 3, 525–535. [Google Scholar]
  18. Freitas, M.J.; Vijayaraghavan, S.; Fardilha, M. Signaling mechanisms in mammalian sperm motility. Biol. Reprod. 2017, 96, 2–12. [Google Scholar] [CrossRef] [Green Version]
  19. Gunes, S.; Sengupta, P.; Henkel, R.; Alguraigari, A.; Sinigaglia, M.M.; Kayal, M.; Joumah, A.; Agarwal, A. Microtubular dysfunction and male infertility. World J. Mens. Health 2020, 38, 9–23. [Google Scholar] [CrossRef]
  20. Martins, M.C.; Gonçalves, L.M.; Nonato, A.; Travençolo, B.A.N.; Alves, B.G.; Beletti, M.E. Sperm head morphometry and chromatin condensation are in constant change at seminiferous tubules, epididymis, and ductus deferens in bulls. Theriogenology 2021, 161, 200–209. [Google Scholar] [CrossRef]
  21. Food and Agriculture Organization Statistical Division (FAOSTAT) of the United Nations. Available online: https://www.fao.org/statistics/en/ (accessed on 1 July 2022).
  22. van Wambeke, F.; Huyghebaert, G. Current role of semen storage and artificial insemination in the turkey industry. Br. Poult. Sci. 1989, 30, 461–469. [Google Scholar] [CrossRef]
  23. Sexton, T.J. Oxidative and glycolytic activity of chicken and turkey spermatozoa. Comp. Biochem. Physiol. Part B Comp. Biochem. 1974, 48, 59–65. [Google Scholar] [CrossRef]
  24. Wishart, G.J. The effect of continuous aeration on the fertility of fowl and Turkey semen stored above 0 °C. Br. Poult. Sci. 1981, 22, 445–450. [Google Scholar] [CrossRef] [PubMed]
  25. Sexton, T.J.; Giesen, A.F. Beltsville Poultry Semen Extender. Poult. Sci. 1982, 61, 1202–1208. [Google Scholar] [CrossRef]
  26. Marzoni, M.; Castillo, A.; Sagona, S.; Citti, L.; Rocchiccioli, S.; Romboli, I.; Felicioli, A. A proteomic approach to identify seminal plasma proteins in roosters (Gallus gallus domesticus). Anim. Reprod. Sci. 2013, 140, 216–223. [Google Scholar] [CrossRef]
  27. Słowińska, M.; Nynca, J.; Arnold, G.J.; Fröhlich, T.; Jankowski, J.; Kozłowski, K.; Mostek, A.; Ciereszko, A. Proteomic identification of turkey (Meleagris gallopavo) seminal plasma proteins. Poult. Sci. 2017, 96, 3422–3435. [Google Scholar] [CrossRef]
  28. Long, J.A. Avian Semen Cryopreservation: What Are the Biological Challenges? Poult. Sci. 2006, 85, 232–236. [Google Scholar] [CrossRef]
  29. Bateman, A.; Martin, M.J.; O’Donovan, C.; Magrane, M.; Alpi, E.; Antunes, R.; Bely, B.; Bingley, M.; Bonilla, C.; Britto, R.; et al. UniProt: The universal protein knowledgebase. Nucleic Acids Res. 2017, 45, D158–D169. [Google Scholar] [CrossRef] [Green Version]
  30. Mańkowska, A.; Brym, P.; Paukszto, Ł.; Jastrzębski, J.P.; Fraser, L. Gene polymorphisms in boar spermatozoa and their associations with post-thaw semen quality. Int. J. Mol. Sci. 2020, 21, 1902. [Google Scholar] [CrossRef] [Green Version]
  31. Froman, D.P.; Feltmann, A.J.; Mclean, D.J. Increased Fecundity Resulting from Semen Donor Selection Based Upon in Vitro Sperm Motility. Poult. Sci. 1997, 76, 73–77. [Google Scholar] [CrossRef]
  32. Froman, D.P.; Feltmann, A.J.; Rhoads, M.L.; Kirby, J.D. Sperm Mobility: A Primary Determinant of Fertility in the Domestic Fowl (Gallus domesticus). Biol. Reprod. 1999, 61, 400–405. [Google Scholar] [CrossRef] [Green Version]
  33. FROMAN, D.P.; McLEAN, D.J. Objective Measurement of Sperm Motility Based Upon Sperm Penetration of Accudenz®. Poult. Sci. 1996, 75, 776–784. [Google Scholar] [CrossRef] [PubMed]
  34. Aire, T.A. Spermiogenesis in birds. Spermatogenesis 2014, 4, e959392. [Google Scholar] [CrossRef] [PubMed]
  35. Manier, M.K.; Welch, G.; Van Nispen, C.; Bakst, M.R.; Long, J. Low-mobility sperm phenotype in the domestic turkey: Impact on sperm morphometry and early embryonic death. Reprod. Domest. Anim. 2019, 54, 613–621. [Google Scholar] [CrossRef]
  36. Kornienko, A.E.; Guenzl, P.M.; Barlow, D.P.; Pauler, F.M. Gene regulation by the act of long non-coding RNA transcription. BMC Biol. 2013, 11, 59. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Asano, A.; Tajima, A. Development and Preservation of Avian Sperm. Adv. Exp. Med. Biol. 2017, 1001, 59–73. [Google Scholar] [PubMed]
  38. Goley, E.D.; Welch, M.D. The ARP2/3 complex: An actin nucleator comes of age. Nat. Rev. Mol. Cell Biol. 2006, 7, 713–726. [Google Scholar] [CrossRef] [PubMed]
  39. Campellone, K.G.; Welch, M.D. A nucleator arms race: Cellular control of actin assembly. Nat. Rev. Mol. Cell Biol. 2010, 11, 237–251. [Google Scholar] [CrossRef] [Green Version]
  40. Dayel, M.J.; Holleran, E.A.; Mullins, R.D. Arp2/3 complex requires hydrolyzable ATP for nucleation of new actin filaments. Proc. Natl. Acad. Sci. USA 2001, 98, 14871–14876. [Google Scholar] [CrossRef] [Green Version]
  41. Dayel, M.J.; Mullins, R.D. Activation of Arp2/3 Complex: Addition of the First Subunit of the New Filament by a WASP Protein Triggers Rapid ATP Hydrolysis on Arp2. PLoS Biol. 2004, 2, e91. [Google Scholar] [CrossRef] [Green Version]
  42. Martin, A.C.; Xu, X.-P.; Rouiller, I.; Kaksonen, M.; Sun, Y.; Belmont, L.; Volkmann, N.; Hanein, D.; Welch, M.; Drubin, D.G. Effects of Arp2 and Arp3 nucleotide-binding pocket mutations on Arp2/3 complex function. J. Cell Biol. 2005, 168, 315–328. [Google Scholar] [CrossRef]
  43. Robinson, R.C.; Turbedsky, K.; Kaiser, D.A.; Marchand, J.-B.; Higgs, H.N.; Choe, S.; Pollard, T.D. Crystal Structure of Arp2/3 Complex. Science 2001, 294, 1679–1684. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Lee, J.S.; Kwon, W.S.; Rahman, M.S.; Yoon, S.J.; Park, Y.J.; Pang, M.G. Actin-related protein 2/3 complex-based actin polymerization is critical for male fertility. Andrology 2015, 3, 937–946. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Cannarella, R.; Condorelli, R.A.; Mongioì, L.M.; La Vignera, S.; Calogero, A.E. Molecular biology of spermatogenesis: Novel targets of apparently idiopathic male infertility. Int. J. Mol. Sci. 2020, 21, 1728. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  46. Xiong, W.; Sun, Y.; Zeng, Z. Antimicrobial use and antimicrobial resistance in food animals. Environ. Sci. Pollut. Res. 2018, 25, 18377–18384. [Google Scholar] [CrossRef]
  47. Kozyreva, V.K.; McLaughlin, S.L.; Livengood, R.H.; Calkins, R.A.; Kelley, L.C.; Rajulapati, A.; Ice, R.J.; Smolkin, M.B.; Weed, S.A.; Pugacheva, E.N. NEDD9 regulates actin dynamics through cortactin deacetylation in an AURKA/HDAC6-dependent manner. Mol. Cancer Res. 2014, 12, 681–693. [Google Scholar] [CrossRef] [Green Version]
  48. Cai, X.; Yu, S.; Mipam, T.D.; Yang, F.; Zhao, W.; Liu, W.; Cao, S.Z.; Shen, L.; Zhao, F.; Sun, L.; et al. Comparative analysis of testis transcriptomes associated with male infertility in cattleyak. Theriogenology 2017, 88, 28–42. [Google Scholar] [CrossRef]
  49. Zhao, B.; Wang, J.; Liu, L.; Li, X.; Liu, S.; Xia, Q.; Shi, J. Annexin A1 translocates to nucleus and promotes the expression of pro-inflammatory cytokines in a PKC-dependent manner after OGD/R. Sci. Rep. 2016, 6, 27028. [Google Scholar] [CrossRef] [Green Version]
  50. Munuce, M.J.; Marini, P.E.; Teijeiro, J.M. Expression profile and distribution of Annexin A1, A2 and A5 in human semen. Andrologia 2019, 51, e13224. [Google Scholar] [CrossRef]
  51. Sargent, K.M.; Lu, N.; Clopton, D.T.; Pohlmeier, W.E.; Brauer, V.M.; Ferrara, N.; Silversides, D.W.; Cupp, A.S. Loss of Vascular Endothelial Growth Factor A (VEGFA) Isoforms in Granulosa Cells Using pDmrt-1-Cre or Amhr2-Cre Reduces Fertility by Arresting Follicular Development and by Reducing Litter Size in Female Mice. PLoS ONE 2015, 10, e0116332. [Google Scholar] [CrossRef] [Green Version]
  52. Sá, R.; Miranda, C.; Carvalho, F.; Barros, A.; Sousa, M. Expression of stem cell markers: OCT4, KIT, ITGA6, and ITGB1 in the male germinal epithelium. Syst. Biol. Reprod. Med. 2013, 59, 233–243. [Google Scholar] [CrossRef]
  53. Wang, J.; Li, J.; Xu, W.; Xia, Q.; Gu, Y.; Song, W.; Zhang, X.; Yang, Y.; Wang, W.; Li, H.; et al. Androgen promotes differentiation of PLZF+ spermatogonia pool via indirect regulatory pattern. Cell Commun. Signal. 2019, 17, 57. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  54. Iwamori, N.; Tominaga, K.; Sato, T.; Riehle, K.; Iwamori, T.; Ohkawa, Y.; Coarfa, C.; Ono, E.; Matzuk, M.M. MRG15 is required for pre-mRNA splicing and spermatogenesis. Proc. Natl. Acad. Sci. USA 2016, 113, E5408–E5415. [Google Scholar] [CrossRef] [Green Version]
  55. Stöckl, J.B.; Schmid, N.; Flenkenthaler, F.; Drummer, C.; Behr, R.; Mayerhofer, A.; Arnold, G.J.; Fröhlich, T. Age-related alterations in the testicular proteome of a non-human primate. Cells 2021, 10, 1306. [Google Scholar] [CrossRef] [PubMed]
  56. Dunleavy, J.E.M.; O’Bryan, M.K.; Stanton, P.G.; O’Donnell, L. The cytoskeleton in spermatogenesis. Reproduction 2019, 157, R53–R72. [Google Scholar] [CrossRef] [Green Version]
  57. Peruquetti, R.L.; Taboga, S.R.; Azeredo-Oliveira, M.T.V.D. Morphological Changes of Mammalian Nucleoli during Spermatogenesis and Their Possible Role in the Chromatoid Body Assembling. ISRN Cell Biol. 2012, 2012, 829854. [Google Scholar] [CrossRef] [Green Version]
  58. Zuo, Y.; Ulu, A.; Chang, J.T.; Frost, J.A. Contributions of the RhoA guanine nucleotide exchange factor Net1 to polyoma middle T antigen-mediated mammary gland tumorigenesis and metastasis. Breast Cancer Res. 2018, 20, 1–16. [Google Scholar] [CrossRef]
  59. Gross, N.; Strillacci, M.G.; Peñagaricano, F.; Khatib, H. Characterization and functional roles of paternal RNAs in 2–4 cell bovine embryos. Sci. Rep. 2019, 9, 20347. [Google Scholar] [CrossRef] [Green Version]
  60. Yoon, J.; Wu, H.; Hung, R.-J.; Terman, J.R. Enhanced Production of the Mical Redox Domain for Enzymology and F-actin Disassembly Assays. Int. J. Mol. Sci. 2021, 22, 1991. [Google Scholar] [CrossRef]
  61. Cabodi, S.; del Pilar Camacho-Leal, M.; Di Stefano, P.; Defilippi, P. Integrin signalling adaptors: Not only figurants in the cancer story. Nat. Rev. Cancer 2010, 10, 858–870. [Google Scholar] [CrossRef]
  62. Chandra, V.; Huang, P.; Potluri, N.; Wu, D.; Kim, Y.; Rastinejad, F. Multidomain integration in the structure of the HNF-4α nuclear receptor complex. Nature 2013, 495, 394–398. [Google Scholar] [CrossRef] [Green Version]
  63. He, J.; Li, L.; Yu, Y.; Hu, X.; Zhang, H.; Liu, R.; Wang, R. Two mutations in the axonemal dynein heavy chain gene 5 in a Chinese asthenozoospermia patient: A case report. Medicine 2020, 99, e20813. [Google Scholar] [CrossRef] [PubMed]
  64. Precone, V.; Cannarella, R.; Paolacci, S.; Busetto, G.M.; Beccari, T.; Stuppia, L.; Tonini, G.; Zulian, A.; Marceddu, G.; Calogero, A.E.; et al. Male Infertility Diagnosis: Improvement of Genetic Analysis Performance by the Introduction of Pre-Diagnostic Genes in a Next-Generation Sequencing Custom-Made Panel. Front. Endocrinol. 2021, 11, 605237. [Google Scholar] [CrossRef] [PubMed]
  65. Fang, D.A.; Yang, Q.Z.; Duan, J.R.; Wang, Q.; Zhang, M.Y.; Zhou, Y.F.; Liu, K.; Shi, W.G. Characteristic of PGDS potential regulation role on spermatogenesis in the Chinese mitten crab Eriocheir sinensis. Gene 2014, 543, 244–252. [Google Scholar] [CrossRef]
  66. Pasini, M.E.; Intra, J.; Gomulski, L.M.; Calvenzani, V.; Petroni, K.; Briani, F.; Perotti, M.E. Identification and expression profiling of Ceratitis capitata genes coding for β-hexosaminidases. Gene 2011, 473, 44–56. [Google Scholar] [CrossRef] [PubMed]
  67. Nixon, B.; Ewen, K.A.; Krivanek, K.M.; Clulow, J.; Kidd, G.; Ecroyd, H.; Jones, R.C. Post-testicular sperm maturation and identification of an epididymal protein in the Japanese quail (Coturnix coturnix japonica). Reproduction 2014, 147, 265–277. [Google Scholar] [CrossRef] [Green Version]
  68. Cheng, W.; Ip, Y.T.; Xu, Z. Gudu, an Armadillo repeat-containing protein, is required for spermatogenesis in Drosophila. Gene 2013, 531, 294–300. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  69. Wu, Y.; Hu, X.; Li, Z.; Wang, M.; Li, S.; Wang, X.; Lin, X.; Liao, S.; Zhang, Z.; Feng, X.; et al. Transcription Factor RFX2 Is a Key Regulator of Mouse Spermiogenesis. Sci. Rep. 2016, 6, 20435. [Google Scholar] [CrossRef]
  70. Wang, Y.Y.; Ke, C.C.; Chen, Y.L.; Lin, Y.H.; Yu, I.S.; Ku, W.C.; O’Bryan, M.K.; Lin, Y.H. Deficiency of the Tbc1d21 gene causes male infertility with morphological abnormalities of the sperm mitochondria and flagellum in mice. PLoS Genet. 2020, 16, e1009020. [Google Scholar] [CrossRef]
  71. Becker-Heck, A.; Zohn, I.E.; Okabe, N.; Pollock, A.; Lenhart, K.B.; Sullivan-Brown, J.; McSheene, J.; Loges, N.T.; Olbrich, H.; Haeffner, K.; et al. The coiled-coil domain containing protein CCDC40 is essential for motile cilia function and left-right axis formation. Nat. Genet. 2011, 43, 79–84. [Google Scholar] [CrossRef] [Green Version]
  72. Antony, D.; Becker-Heck, A.; Zariwala, M.A.; Schmidts, M.; Onoufriadis, A.; Forouhan, M.; Wilson, R.; Taylor-Cox, T.; Dewar, A.; Jackson, C.; et al. Mutations in CCDC39 and CCDC40 are the Major Cause of Primary Ciliary Dyskinesia with Axonemal Disorganization and Absent Inner Dynein Arms. Hum. Mutat. 2013, 34, 462–472. [Google Scholar] [CrossRef] [Green Version]
  73. Sheybani, N.; Bakhtiarizadeh, M.R.; Salehi, A. An integrated analysis of mRNAs, lncRNAs, and miRNAs based on weighted gene co-expression network analysis involved in bovine endometritis. Sci. Rep. 2021, 11, 18050. [Google Scholar] [CrossRef] [PubMed]
  74. Wang, K.C.; Chang, H.Y. Molecular Mechanisms of Long Noncoding RNAs. Mol. Cell 2011, 43, 904–914. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  75. Zhang, X.; Wang, W.; Zhu, W.; Dong, J.; Cheng, Y.; Yin, Z.; Shen, F. Mechanisms and Functions of Long Non-Coding RNAs at Multiple Regulatory Levels. Int. J. Mol. Sci. 2019, 20, 5573. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  76. Gil, N.; Ulitsky, I. Regulation of gene expression by cis-acting long non-coding RNAs. Nat. Rev. Genet. 2020, 21, 102–117. [Google Scholar] [CrossRef] [PubMed]
  77. Gong, J.; Shao, D.; Xu, K.; Lu, Z.; Lu, Z.J.; Yang, Y.T.; Zhang, Q.C. RISE: A database of RNA interactome from sequencing experiments. Nucleic Acids Res. 2018, 46, D194–D201. [Google Scholar] [CrossRef] [PubMed]
  78. Graindorge, A.; Pinheiro, I.; Nawrocka, A.; Mallory, A.C.; Tsvetkov, P.; Gil, N.; Carolis, C.; Buchholz, F.; Ulitsky, I.; Heard, E.; et al. In-cell identification and measurement of RNA-protein interactions. Nat. Commun. 2019, 10, 5317. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  79. Statello, L.; Guo, C.J.; Chen, L.L.; Huarte, M. Gene regulation by long non-coding RNAs and its biological functions. Nat. Rev. Mol. Cell Biol. 2021, 22, 96–118. [Google Scholar] [CrossRef]
  80. Wang, M.; Liu, X.; Chang, G.; Chen, Y.; An, G.; Yan, L.; Gao, S.; Xu, Y.; Cui, Y.; Dong, J.; et al. Single-Cell RNA Sequencing Analysis Reveals Sequential Cell Fate Transition during Human Spermatogenesis. Cell Stem Cell 2018, 23, 599–614. [Google Scholar] [CrossRef] [Green Version]
  81. Słowińska, M.; Liszewska, E.; Nynca, J.; Bukowska, J.; Hejmej, A.; Bilińska, B.; Szubstarski, J.; Kozłowski, K.; Jankowski, J.; Ciereszko, A. Isolation and Characterization of an Ovoinhibitor, a Multidomain Kazal-Like Inhibitor from Turkey (Meleagris gallopavo) Seminal Plasma1. Biol. Reprod. 2014, 91, 108. [Google Scholar] [CrossRef]
  82. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 2014, 30, 2114–2120. [Google Scholar] [CrossRef] [Green Version]
  83. Dobin, A.; Davis, C.A.; Schlesinger, F.; Drenkow, J.; Zaleski, C.; Jha, S.; Batut, P.; Chaisson, M.; Gingeras, T.R. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 2013, 29, 15–21. [Google Scholar] [CrossRef]
  84. Pertea, M.; Pertea, G.M.; Antonescu, C.M.; Chang, T.-C.; Mendell, J.T.; Salzberg, S.L. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 2015, 33, 290–295. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  85. Frazee, A.C.; Pertea, G.; Jaffe, A.E.; Langmead, B.; Salzberg, S.L.; Leek, J.T. Ballgown bridges the gap between transcriptome assembly and expression analysis. Nat Biotechnol. 2015, 33, 243–246. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  86. Trapnell, C.; Roberts, A.; Goff, L.; Pertea, G.; Kim, D.; Kelley, D.R.; Pimentel, H.; Salzberg, S.L.; Rinn, J.L.; Pachter, L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat. Protoc. 2012, 7, 562–578. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  87. R Core Team. R: A Language and Environment for Statistical Computing, version 4.2.1; R Foundation for Statistical Computing: Vienna, Austria, 2013; Available online: http://www.R-project.org/ (accessed on 20 June 2022).
  88. Kang, Y.J.; Yang, D.C.; Kong, L.; Hou, M.; Meng, Y.Q.; Wei, L.; Gao, G. CPC2: A fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res. 2017, 45, W12–W16. [Google Scholar] [CrossRef] [Green Version]
  89. Wucher, V.; Legeai, F.; Hédan, B.; Rizk, G.; Lagoutte, L.; Leeb, T.; Jagannathan, V.; Cadieu, E.; David, A.; Lohi, H.; et al. FEELnc: A tool for long non-coding RNA annotation and its application to the dog transcriptome. Nucleic Acids Res. 2017, 45, e57. [Google Scholar] [CrossRef] [Green Version]
  90. Wang, L.; Park, H.J.; Dasari, S.; Wang, S.; Kocher, J.-P.; Li, W. CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. Nucleic Acids Res. 2013, 41, e74. [Google Scholar] [CrossRef]
  91. Li, A.; Zhang, J.; Zhou, Z. PLEK: A tool for predicting long non-coding RNAs and messenger RNAs based on an improved k-mer scheme. BMC Bioinform. 2014, 15, 311. [Google Scholar] [CrossRef] [Green Version]
  92. Sun, L.; Luo, H.; Bu, D.; Zhao, G.; Yu, K.; Zhang, C.; Liu, Y.; Chen, R.; Zhao, Y. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013, 41, e166. [Google Scholar] [CrossRef]
  93. Han, S.; Liang, Y.; Ma, Q.; Xu, Y.; Zhang, Y.; Du, W.; Wang, C.; Li, Y. LncFinder: An integrated platform for long non-coding RNA identification utilizing sequence intrinsic composition, structural information and physicochemical property. Brief. Bioinform. 2019, 20, 2009–2027. [Google Scholar] [CrossRef]
  94. Fan, X.N.; Zhang, S.W.; Zhang, S.Y.; Ni, J.J. Lncrna_mdeep: An alignment-free predictor for distinguishing long non-coding rnas from protein-coding transcripts by multimodal deep learning. Int. J. Mol. Sci. 2020, 21, 5222. [Google Scholar] [CrossRef] [PubMed]
  95. Li, J.; Ma, W.; Zeng, P.; Wang, J.; Geng, B.; Yang, J.; Cui, Q. LncTar: A tool for predicting the RNA targets of long noncoding RNAs. Brief. Bioinform. 2014, 16, 806–812. [Google Scholar] [CrossRef] [PubMed]
  96. Lu, Q.; Ren, S.; Lu, M.; Zhang, Y.; Zhu, D.; Zhang, X.; Li, T. Computational prediction of associations between long non-coding RNAs and proteins. BMC Genom. 2013, 14, 651. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  97. Szklarczyk, D.; Gable, A.L.; Lyon, D.; Junge, A.; Wyder, S.; Huerta-Cepas, J.; Simonovic, M.; Doncheva, N.T.; Morris, J.H.; Bork, P.; et al. STRING v11: Protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019, 47, D607–D613. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  98. Shannon, P.; Markiel, A.; Ozier, O.; Baliga, N.S.; Wang, J.T.; Ramage, D.; Amin, N.; Schwikowski, B.; Ideker, T. Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res. 2003, 13, 2498–2504. [Google Scholar] [CrossRef]
  99. Ahmed, M.; Kim, D.R. pcr: An R package for quality assessment, analysis and testing of qPCR data. PeerJ 2018, 6, e4473. [Google Scholar] [CrossRef]
  100. Yuan, J.S.; Reed, A.; Chen, F.; Stewart, C.N. Statistical analysis of real-time PCR data. BMC Bioinform. 2006, 7, 85. [Google Scholar] [CrossRef] [Green Version]
Ijms 23 07642 g001 550
Figure 1. The flowchart of the lncRNA identification process. The process consists of several steps (described in the main text) grouped in the three stages: “input data”—preparing and reading the data from the previous steps (colored orange) and reference data (colored blue), “filtering steps”—filtering the transcripts basing on the annotation (i.e., extracting known lncRNAs and removing protein-coding transcripts) and structural features such as sequence length and the number of exons and “coding probability”—screening the sequence dataset based on the coding probability.
Figure 1. The flowchart of the lncRNA identification process. The process consists of several steps (described in the main text) grouped in the three stages: “input data”—preparing and reading the data from the previous steps (colored orange) and reference data (colored blue), “filtering steps”—filtering the transcripts basing on the annotation (i.e., extracting known lncRNAs and removing protein-coding transcripts) and structural features such as sequence length and the number of exons and “coding probability”—screening the sequence dataset based on the coding probability.
Ijms 23 07642 g001
Ijms 23 07642 g002 550
Figure 2. The Venn plots show the number of transcripts in the lncRNA identification process in the two stages: (A) filtering by basic parameters: expression level (expr), number of exons, sequence length, annotation to protein-coding genes (nonPC), (B) numbers of sequences predicted as non-coding using seven tools.
Figure 2. The Venn plots show the number of transcripts in the lncRNA identification process in the two stages: (A) filtering by basic parameters: expression level (expr), number of exons, sequence length, annotation to protein-coding genes (nonPC), (B) numbers of sequences predicted as non-coding using seven tools.
Ijms 23 07642 g002
Ijms 23 07642 g003 550
Figure 3. MA and volcano plots of expressed transcripts (presented on a gene level) in three comparisons: (A,D) epididymis vs. testis, (B,E) ductus deferens vs. testis and (C,F) epididymis vs. ductus deferens. Each dot presents expressed transcripts, where dark grey is DEGs, red is upregulated lncRNAs and green is downregulated lncRNAs. The triangles and rectangles indicate off the chart values.
Figure 3. MA and volcano plots of expressed transcripts (presented on a gene level) in three comparisons: (A,D) epididymis vs. testis, (B,E) ductus deferens vs. testis and (C,F) epididymis vs. ductus deferens. Each dot presents expressed transcripts, where dark grey is DEGs, red is upregulated lncRNAs and green is downregulated lncRNAs. The triangles and rectangles indicate off the chart values.
Ijms 23 07642 g003
Ijms 23 07642 g006 550
Figure 6. Trans-correlations (A) and direct lncRNA-mRNA interactions (B) identified in the group of processes with significantly expressed genes between DEGs and DELs associated with the E/(T + DD) group of molecular functions (see Table 4 and Figure 4). The first column in Figure A represents the biological processes. The second column shows the genes associated with the processes in the first column (the lines of the Sankey plot colored corresponding to the processes represent the gene–process relationship). The third column presents lncRNAs, and the lines connecting them with genes symbolize the expression correlation (blue—positive correlation, red—negative correlation). The green color in Figure B indicates lncRNAs, and red indicates the protein-coding genes. The lines symbolize the predicted direct interactions of the lncRNA-mRNA.
Figure 6. Trans-correlations (A) and direct lncRNA-mRNA interactions (B) identified in the group of processes with significantly expressed genes between DEGs and DELs associated with the E/(T + DD) group of molecular functions (see Table 4 and Figure 4). The first column in Figure A represents the biological processes. The second column shows the genes associated with the processes in the first column (the lines of the Sankey plot colored corresponding to the processes represent the gene–process relationship). The third column presents lncRNAs, and the lines connecting them with genes symbolize the expression correlation (blue—positive correlation, red—negative correlation). The green color in Figure B indicates lncRNAs, and red indicates the protein-coding genes. The lines symbolize the predicted direct interactions of the lncRNA-mRNA.
Ijms 23 07642 g006
Ijms 23 07642 g007 550
Figure 7. The STRING network of proteins involved in sperm motility. DD/(T + E), E/(T + DD) and T/(E + DD) groups correspond to the three groups of analyzed processes (see Table 4 and Figure 5). The color of the lines symbolizes the type of interaction and the source of information: light blue and pink-known interactions, green, red and blue-predicted interactions, yellow—text mining, black—co-expression, and purple—protein homology.
Figure 7. The STRING network of proteins involved in sperm motility. DD/(T + E), E/(T + DD) and T/(E + DD) groups correspond to the three groups of analyzed processes (see Table 4 and Figure 5). The color of the lines symbolizes the type of interaction and the source of information: light blue and pink-known interactions, green, red and blue-predicted interactions, yellow—text mining, black—co-expression, and purple—protein homology.
Ijms 23 07642 g007
Ijms 23 07642 g008 550
Figure 8. The real-time PCR results. The statistical significance of the expression differences was confirmed at the level <0.05 and is marked with a star. Data are shown as mean values ± SD.
Figure 8. The real-time PCR results. The statistical significance of the expression differences was confirmed at the level <0.05 and is marked with a star. Data are shown as mean values ± SD.
Ijms 23 07642 g008
Ijms 23 07642 g009 550
Figure 9. A visualization of the 3D model of the human HNF4A protein (PDB id 4IQR). The protein part is presented as ribbons, and the conserved DNA-binding domains are colored red.
Figure 9. A visualization of the 3D model of the human HNF4A protein (PDB id 4IQR). The protein part is presented as ribbons, and the conserved DNA-binding domains are colored red.
Ijms 23 07642 g009
Table
Table 1. Overall statistics of sequenced, mapped and processed data. The read values are presented in millions. Min/max are the minimum and maximum values. Mean values are calculated for all six samples for each tissue.
Table 1. Overall statistics of sequenced, mapped and processed data. The read values are presented in millions. Min/max are the minimum and maximum values. Mean values are calculated for all six samples for each tissue.
TEDD
Min/MaxMeanMin/MaxMeanMin/MaxMean
Raw reads127.9/165.8146.9126.4/176.9152.0102.6/167.7145.9
Trimmed reads110.7/140.6123.896/151.9127.586.2/137.6121.4
% of trimmed reads81.1/86.584.376/87.283.582/84.783.2
Mapped reads78/10087.664.8/117.393.062.3/99.984.8
Uniquely mapped reads76.2/97.685.962.6/115.891.661.3/98.483.5
% of uniquely mapped reads97.6/98.698.0096.6/98.998.398.2/98.698.4
Expressed transcripts59,151/73,56967,594.322,109/61,86847,862.833,355/47,86141,412.5
Expressed genes36,879/43,51040,405.317,596/34,71128,438.521,196/27,52624,970.2
Table
Table 2. The number of up- and down-regulated differentially expressed lncRNAs in three comparisons.
Table 2. The number of up- and down-regulated differentially expressed lncRNAs in three comparisons.
E vs. TDD vs. TE vs. DD
All13918225
Up-regulated9313225
Down-regulated46500
Table
Table 3. Differentially expressed genes (DEGs) associated with molecular processes engaged in sperm motility and grouped into three sets (T/(E + DD), E/(T + DD) and DD/(T + E) based on expression profiles presented in Figure 5.
Table 3. Differentially expressed genes (DEGs) associated with molecular processes engaged in sperm motility and grouped into three sets (T/(E + DD), E/(T + DD) and DD/(T + E) based on expression profiles presented in Figure 5.
GroupProcessGenes
T/(E + DD)actin binding
actin filament binding
actin cytoskeleton organization
actin filament-based process
actin filament organization
regulation of actin filament length
regulation of protein phosphorylation
apoptotic process
response to oxygen-containing compound
response to stress
solute: cation symporter activity
response to estrogen		ARPC5, CALD1, CNN3, EVL, EZR, GSN, ITGB1, MICAL2, MYH11, PARVA, POF1B, SPTBN1, TWF1, USH1C, CARMIL1, CGNL1, EPHA1, FAT1, MET, NRP1, PACSIN2, PDLIM3, PHLDB2, SDCBP, SLC9A3R1, SLIT2, SORBS1, WDR1, ANXA1, IQGAP2, ACVR1, APP, BMP7, CTNNB1, HBEGF, NEDD9, PKD2, SPRY2, THBS1, ACKR3, BAG3, GATA3, LEF1, NET1, PTPRK, ATP1A1, CAT, DCN, MMP2, GCNT1, SOD, PRDX1, ATP6V0A4, PTPPK, FIGF, MYLK, SLC15A1, SLC34A2, KRT19
E/(T + DD)calcium ion binding
reproductive process
lipid metabolism		CDHR2, CLGN, HPGDS, LRP2, CDHR4, PLA2G12B, CDH20, DNAH5, HEXB, HNF4A, RDH10
DD/(T + E)axonemal dynein complex assembly
microtubule cytoskeleton organization
microtubule-based movement
inner dynein arm assembly
outer dynein arm assembly
flagellated sperm motility		ARMC4, CCDC40, DNAH7, LRRC6, PIH1D3, ZMYND10, DEUP1, NDC80, PLK4, RSPH9, TTC26, CFAP221, ROPN1L, CFAP69
Table
Table 4. Cumulative correlations between DELs and DEGs. Significant correlations between lncRNA and mRNA or proteins in each comparison (number of interactions associated with the sperm motility/all identified interactions).
Table 4. Cumulative correlations between DELs and DEGs. Significant correlations between lncRNA and mRNA or proteins in each comparison (number of interactions associated with the sperm motility/all identified interactions).
Relations/InteractionsCisTranslncRNA-mRNAlncRNA-Protein
E vs. T0/23519/15,1903302/10,1715845/9712
DD vs. T0/30749/22,9834875/13,9387335/14,172
E vs. DD0/199/759221/317405/512
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Jastrzebski, J.P.; Lipka, A.; Majewska, M.; Makowczenko, K.G.; Paukszto, L.; Bukowska, J.; Dorocki, S.; Kozlowski, K.; Slowinska, M. In Silico Identification of lncRNAs Regulating Sperm Motility in the Turkey (Meleagris gallopavo L.). Int. J. Mol. Sci. 2022, 23, 7642. https://doi.org/10.3390/ijms23147642

AMA Style

Jastrzebski JP, Lipka A, Majewska M, Makowczenko KG, Paukszto L, Bukowska J, Dorocki S, Kozlowski K, Slowinska M. In Silico Identification of lncRNAs Regulating Sperm Motility in the Turkey (Meleagris gallopavo L.). International Journal of Molecular Sciences. 2022; 23(14):7642. https://doi.org/10.3390/ijms23147642

Chicago/Turabian Style

Jastrzebski, Jan Pawel, Aleksandra Lipka, Marta Majewska, Karol G. Makowczenko, Lukasz Paukszto, Joanna Bukowska, Slawomir Dorocki, Krzysztof Kozlowski, and Mariola Slowinska. 2022. "In Silico Identification of lncRNAs Regulating Sperm Motility in the Turkey (Meleagris gallopavo L.)" International Journal of Molecular Sciences 23, no. 14: 7642. https://doi.org/10.3390/ijms23147642

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop