Characterization of Alternative Splicing (AS) Events during Chicken (Gallus gallus) Male Germ-Line Stem Cell Differentiation with Single-Cell RNA-seq

Simple Summary Studies have shown that alternative splicing (AS) has been utilized in a wide variety of life processes. However, there are very few studies on AS during germ cell development. In this study, we preliminarily investigated the variation of variable shear events during the formation of chicken germ cells through the RNA-seq data analysis of embryonic stem cells (ESCs), gonad PGCs (gPGCs), and spermatogonia stem cells (SSCs), and the critical AS mode for several crucial stage-specific genes, which were identified during germ cell development. The results of this study lay a theoretical foundation for further analysis of the regulation mechanism of key genes involved in germ cell formation. Abstract Alternative splicing (AS) is a ubiquitous, co-transcriptional, and post-transcriptional regulation mechanism during certain developmental processes, such as germ cell differentiation. A thorough understanding of germ cell differentiation will help us to open new avenues for avian reproduction, stem cell biology, and advances in medicines for human consumption. Here, based on single-cell RNA-seq, we characterized genome-wide AS events in manifold chicken male germ cells: embryonic stem cells (ESCs), gonad primordial germ cells (gPGCs), and spermatogonia stem cells (SSCs). A total of 38,494 AS events from 15,338 genes were detected in ESCs, with a total of 48,955 events from 14,783 genes and 49,900 events from 15,089 genes observed in gPGCs and SSCs, respectively. Moreover, this distribution of AS events suggests the diverse splicing feature of ESCs, gPGCs, and SSCs. Finally, several crucial stage-specific genes, such as NANOG, POU5F3, LIN28B, BMP4, STRA8, and LHX9, were identified in AS events that were transmitted in ESCs, gPGCs, and SSCs. The gene expression results of the RNA-seq data were validated by qRT-PCR. In summary, we provided a comprehensive atlas of the genome-wide scale of the AS event landscape in male chicken germ-line cells and presented its distribution for the first time. This research may someday improve treatment options for men suffering from male infertility.


Introduction
Alternative splicing (AS) is an important, co-transcriptional, and post-transcriptional regulatory mechanism that results in multiple mRNA and protein isoforms from a single

RNA-Seq Reads Mapping
Filtering and quality control checks of raw reads from RNA-seq were performed utilizing the Fast QC program (http://www.bioinformatics.babraham.ac.uk/projects/ fastqc/, accessed on 1 August 2019). The clean reads were mapped to the chicken reference genome (Gallus gallus, galGal4.0) using Tophat (version 2.1.1) [32]. The mapped reads were assembled into contigs with Cufflinks (version 2.2.1). The expression value for each gene was calculated and normalized to FPKM (fragments per kilobase of exon per million fragments mapped) values using Cuffnorm (version 2.2.1) in the Cufflinks package [33].

Identification, Classification, and Quantification of AS Events
After utilizing Cufflinks, the identification and quantification of AS events was accomplished using ASprofile [34]. Initially, the extract AS program compares all pairs of transcripts within a gene to determine exon-intron structural differences that indicate an AS event. The following types of events have been implemented: exon skipping (SKIP), cassette exons (MSKIP), alternative transcript start and termination (TSS, TTS), retention of single or multiple introns (IR, MIR), and alternative exons (AEs). Additionally, variations on these classes have allowed for some wiggle room at the boundaries of the surrounding features, and they have been reported and prefixed with an X. Therefore, we classified these AS events into 5 types: alternative 5 splice Site (A5 , including TSS and XTSS), alternative 3 splice Site (A3 , including TTS and XTTS), exon skip (ES, including SKIP, MSKIP, XSKIP, and XMSKIP), retained intron (RI, including IR, XMIR, XIR, and XMIR), and alternative exon (AE, including AE and XAE). Finally, the extract-as-fpkm program calculates the fragments per kilobase million (FPKM) of each event from those of transcripts harboring an event in each sample. The results can be used to observe trends in dynamics of AS profiles.

Intersections and Interaction Network Analysis of the AS Event Gene
The intersections (Upset plot and Veen plot) between different types of AS events were also investigated and illustrated by employing the UpSetR package (version 1.4.0, https: //www.r-project.org/web/packages/UpSetR, accessed on 22 May 2019) of R software (version 3.6.2, https://www.r-project.org/, accessed on 12 December 2019). After selection by Veen plot, the genes of AS events were submitted to String (version 11.0, https://stringdb.org/, accessed on 19 January 2019) for the protein-protein interaction (PPI) analysis. Cystoscope (version 3.7.2, https://cytoscape.org/, accessed on 5 September 2019) was then applied to illustrate and identify neighborhoods where proteins were densely connected in the PPI network.

GO Enrichment and KEGG Pathway Analysis of AS Events Gene
Gene ontology (GO) enrichment analysis was performed using DAVID (v6.8, https: //david.ncifcrf.gov/, accessed on 1 September 2020) [35,36], and GO terms were considered as statistically significant with p-values < 0.05. GO terms have been classified by the GO database (ftp://ftp.ncbi.nih.gov/gene/DATA/gene2-go.gz, accessed on 10 May 2017), which is an international standard gene functional classification system. In addition, BinGO, which illustrates relationships among GO terms, was explained using Cytoscape 3.7.2 version [37]. The Kyoto Encyclopedia of Genes and Genomes (KEGG) helps researchers to better understand the biological functions of genes based on large-scale molecular datasets (http://www.genome.jp/kegg/, accessed on 10 May 2017). In this study, the enrichment of KEGG was analyzed by KOBAS 3.0 (http://kobas.cbi.pku.edu.cn/kobas3 /annotate/, accessed on 20 September 2018) and a p-value < 0.05 was set as the cutoff for significantly enriched KEGG terms. Illustrations were performed through the ggplot2 package in R [38].

RNA Preparation and RT-qPCR
Isolated cells were homogenized using TRIzol, and the total RNA was isolated according to the manufacturer's instructions (QIAGEN, Beijing, China, DP424). qRT-PCR was performed using the FastKing One-Step RT-PCR Kit with SYBR green (QIAGEN, Beijing, China, KR123), and the mRNA levels were determined by a CFX-Connect Real-time PCR detection system (BIO-RAD, Hercules, CA, USA, 7500fast). Quantification results were normalized to the relative housing keeping gene (GAPDH) using the 2-∆∆Ct method. The sequences of qRT-PCR primers are listed in Table 1.

The AS Events during Male Chicken Germ Cell Differentiation
To study AS events in male chicken germ cell differentiation, we performed genomewide RNA-seq for ESCs, gPGCs, and SSCs isolated from chicken embryonic blastoderms (E0, HH stage X), gonads (E4.5, HH stage [28][29][30], and testes (E18.5, HH stage 44), respectively ( Figure 1a). Our results showed that the sequencing of the three cells transcriptome provided 195.7 million reads of ESCs, 170.1 million reads of gPGCs, and 143.2 million reads of SSCs. After quality control and mapping to the entire annotated chicken genome (galGal 4.0), 28.7%, 31.47% and 28.13% mapped splice reads were observed in ESCs, gPGCs, and SSCs, respectively (see Figure 1b). To refine and explore AS event changes during germ cell differentiation, we divided AS events into alternative 5 splice site (A5 ), alternative 3 splice site (A3 ), exon skip (ES), retained intron (RI), and alternative exon (AE) (see Figure 1c). In ESCs, a total of 38,494 AS events from 15,338 genes were detected, indicating that each gene might have had 2-3 AS events on average. However, 48,955 AS events from 14,783 genes and 49,900 AS events from 15,089 genes were found in gPGCs and SSCs, respectively. This indicated an average of 3-4 AS events from each gene, which was a notably higher frequency than that in ESCs (see Table 2). In detail, most AS events detected during germ cell differentiation were A5 and A3 followed by ES, RI, and AE, which further implies that gPGCs and SSCs had a higher frequency of AS events than ESCs across all five kinds of AS events, especially ES, RI, and AE (see Figure 1d). Moreover, ES, RI, and AE were much less prevalent than A5 /A3 in three germ cell types, which suggests that A5 /A3 are major AS events during chicken male germ cell differentiation.
provided 195.7 million reads of ESCs, 170.1 million reads of gPGCs, and 143.2 million reads of SSCs. After quality control and mapping to the entire annotated chicken genome (galGal 4.0), 28.7%, 31.47% and 28.13% mapped splice reads were observed in ESCs, gPGCs, and SSCs, respectively (see Figure 1b). To refine and explore AS event changes during germ cell differentiation, we divided AS events into alternative 5′ splice site (A5′), alternative 3′ splice site (A3′), exon skip (ES), retained intron (RI), and alternative exon (AE) (see Figure 1c). In ESCs, a total of 38,494 AS events from 15,338 genes were detected, indicating that each gene might have had 2-3 AS events on average. However, 48,955 AS events from 14,783 genes and 49,900 AS events from 15,089 genes were found in gPGCs and SSCs, respectively. This indicated an average of 3-4 AS events from each gene, which was a notably higher frequency than that in ESCs (see Table 2). In detail, most AS events detected during germ cell differentiation were A5′ and A3′ followed by ES, RI, and AE, which further implies that gPGCs and SSCs had a higher frequency of AS events than ESCs across all five kinds of AS events, especially ES, RI, and AE (see Figure 1d). Moreover, ES, RI, and AE were much less prevalent than A5′/A3′ in three germ cell types, which suggests that A5′/A3′ are major AS events during chicken male germ cell differentiation.

AS Event Distribution Analysis during Male Chicken Germ Cell Differentiation
To investigate the splicing feature of each germ cell, we generalized the common gene distribution of AS events during male germ cell differentiation. The results showed that the A5 /A3 gene group contained 12,453, 8081, and 9031 genes in ESCs, gPGCs, and SSCs, respectively. This was the intersection group with the largest number of genes involved in two types of events. However, the A5 /A3 /ES and AT groups contained 749 and 390 genes which followed the A5 /A3 group in ESCs. In gPGCs, the following groups were the A5 /A3 groups, which contained 1593 and 1338 genes, respectively. Additionally, SSCs, accounting for 1413 and 1175 genes, were contained in the A5 /A3 /ES group, and the A3 group followed the A5 /A3 group. In each type of male germ cell, different splicing features were present; more details can be found in Figure 2 and Tables S1-S3. Overall, the male germ cell exhibited different distribution characteristics of AS events, suggesting that AS is involved in the male germ cell differentiation.
gene distribution of AS events during male germ cell differentiation. The results showed that the A5′/A3′ gene group contained 12,453, 8081, and 9031 genes in ESCs, gPGCs, and SSCs, respectively. This was the intersection group with the largest number of genes involved in two types of events. However, the A5′/A3′/ES and AT groups contained 749 and 390 genes which followed the A5′/A3′ group in ESCs. In gPGCs, the following groups were the A5′/A3′ groups, which contained 1593 and 1338 genes, respectively. Additionally, SSCs, accounting for 1413 and 1175 genes, were contained in the A5′/A3′/ES group, and the A3′ group followed the A5′/A3′ group. In each type of male germ cell, different splicing features were present; more details can be found in Figure 2 and Tables S1-S3. Overall, the male germ cell exhibited different distribution characteristics of AS events, suggesting that AS is involved in the male germ cell differentiation.

A Cell-Type-Specific Spliced Gene in ESCs, GPGCs, and SSCs
To find out the potential influence of AS in male germ cell gene expressions, we compared the total gene and spliced gene expression patterns in ESCs, gPGCs, and SSCs. They shared a similar expression pattern (see Figure 3a; Tables S4 and S5). In addition to the expression pattern comparison at the gene level, we studied differentially spliced genes during the male germ cell differentiation and compared AS genes of three germ cells. Their overlap had 14,168 genes, and 222, 83, and 79 of these genes were singlehandedly present in ESCs, gPGCs, and SSCs, respectively (see Figure 3b). To determine whether the genes were affected by male germ cell differentiation, we analyzed these cell-specific gene groups. The network results revealed that ESCs might contain four interaction networks, which include hub genes at MYC, HOXA, RARS, etc. However, the GPGCs and SSCs may exhibit fewer networks and hub genes (see Figure 3c). Taken together, these results imply that differentially spliced genes appear in male germ cell differentiations. pared the total gene and spliced gene expression patterns in ESCs, gPGCs, and SSCs. They shared a similar expression pattern (see Figure 3a; Tables S4 and S5). In addition to the expression pattern comparison at the gene level, we studied differentially spliced genes during the male germ cell differentiation and compared AS genes of three germ cells. Their overlap had 14,168 genes, and 222, 83, and 79 of these genes were singlehandedly present in ESCs, gPGCs, and SSCs, respectively (see Figure 3b). To determine whether the genes were affected by male germ cell differentiation, we analyzed these cell-specific gene groups. The network results revealed that ESCs might contain four interaction networks, which include hub genes at MYC, HOXA, RARS, etc. However, the GPGCs and SSCs may exhibit fewer networks and hub genes (see Figure 3c). Taken together, these results imply that differentially spliced genes appear in male germ cell differentiations.

Function Enrichment Results of Spliced Genes during Male Germ Cell Differentiation
To verify the function of these spliced genes, we performed GO and KEGG pathway enrichment analyses universally for the process of spliced genes during male germ cell differentiation. In total, 22 first-class GO terms were significantly identified: 7 GO terms were enriched in biological processes, including cellular process (GO:0009987), developmental process (GO:0032502), and metabolic process (GO:0008152). Additionally, we found 10 and 5 GO terms in the cellular component and molecular function parts, respectively (see Figure 4a, Table S6). Moreover, GO terms such as cellular process (GO:0009987), developmental process (GO:0032502), and metabolic process (GO:0008152) stand out, which play important roles in the biological process (Figure 4b). A total of 166 KEGG pathways were significantly enriched by the overall spliced genes (Table S7). The top 30 pathways are shown in Figure 4c, containing metabolic pathways (gga01100), endocytosis (gga04144), MAPK signaling pathway (gga04010), calcium signaling pathway (gga04218), Wnt signaling pathway (gga04310), and spliceosome (gga03040) (Figure 4c). These results indicate that spliced genes widely exist in the cell differentiation process.

Function Enrichment Results of Spliced Genes during Male Germ Cell Differentiation
To verify the function of these spliced genes, we performed GO and KEGG pathway enrichment analyses universally for the process of spliced genes during male germ cell differentiation. In total, 22 first-class GO terms were significantly identified: 7 GO terms were enriched in biological processes, including cellular process (GO:0009987), developmental process (GO:0032502), and metabolic process (GO:0008152). Additionally, we found 10 and 5 GO terms in the cellular component and molecular function parts, respectively (see Figure 4a, Table S6). Moreover, GO terms such as cellular process (GO:0009987), developmental process (GO:0032502), and metabolic process (GO:0008152) stand out, which play important roles in the biological process (Figure 4b). A total of 166 KEGG pathways were significantly enriched by the overall spliced genes (Table S7). The top 30 pathways are shown in Figure 4c, containing metabolic pathways (gga01100), endocytosis (gga04144), MAPK signaling pathway (gga04010), calcium signaling pathway (gga04218), Wnt signaling pathway (gga04310), and spliceosome (gga03040) (Figure 4c). These results indicate that spliced genes widely exist in the cell differentiation process.

Validation of AS Events and Crucial Gene Expression
Even previous results showed that spliced genes have similar expression patterns to total genes, but details of the relationship between AS events and spliced genes are still unknown. We displayed AS events from six genes which play a crucial role in male chicken germ cell differentiation. NANOG and POU5F3, highly expressed in ESCs, exhibited A5 and A3 events [39,40]. The LIN28B had more types of AS events in gPGCs, and both of the BMP4 had specific expression in gPGCs [41,42] (see Figure 5a). STRA8 and LHX9 exhibited high expressions in SSCs and showed different patterns of AS events during male germ cell differentiation [43,44] (see Figure 5b). Transcript levels of these genes in ESCs, gPGCs, and

Validation of AS Events and Crucial Gene Expression
Even previous results showed that spliced genes have similar expression patterns to total genes, but details of the relationship between AS events and spliced genes are still unknown. We displayed AS events from six genes which play a crucial role in male chicken germ cell differentiation. NANOG and POU5F3, highly expressed in ESCs, exhibited A5′ and A3′ events [39,40]. The LIN28B had more types of AS events in gPGCs, and both of the BMP4 had specific expression in gPGCs [41,42] (see Figure 5a). STRA8 and LHX9 exhibited high expressions in SSCs and showed different patterns of AS events during male germ cell differentiation [43,44] (see Figure 5b). Transcript levels of these genes in ESCs, gPGCs, and SSCs were measured by qRT-PCR. The consistent results from both qRT-PCR and RNA-Seq suggested the reliability of RNA-Seq data and its accuracy in quantifying gene expressions in male germ cells (see Figure 5c).

Validation of AS Events and Crucial Gene Expression
Even previous results showed that spliced genes have similar expression patterns to total genes, but details of the relationship between AS events and spliced genes are still unknown. We displayed AS events from six genes which play a crucial role in male chicken germ cell differentiation. NANOG and POU5F3, highly expressed in ESCs, exhibited A5′ and A3′ events [39,40]. The LIN28B had more types of AS events in gPGCs, and both of the BMP4 had specific expression in gPGCs [41,42] (see Figure 5a). STRA8 and LHX9 exhibited high expressions in SSCs and showed different patterns of AS events during male germ cell differentiation [43,44] (see Figure 5b). Transcript levels of these genes in ESCs, gPGCs, and SSCs were measured by qRT-PCR. The consistent results from both qRT-PCR and RNA-Seq suggested the reliability of RNA-Seq data and its accuracy in quantifying gene expressions in male germ cells (see Figure 5c).

Discussion
The phenomenon of AS was discovered in the twentieth century, but it has not been systemically analyzed [45]. Advances in high-throughput sequencing provide an efficient technology to acquire AS events at the genome-wide level [46]. In this study, we reported and characterized genome-wide AS events during male chicken germ cell differentiation through a systemic comparison of single-cell RNA-seq. Results showed that AS events, especially A5 /A3 events, are comprehensively involved in the process of male germ cell differentiation, while showing different patterns in ESCs, gPGCs, and SSCs. Functional enrichments in GO and KEGG pathways are closely related to cell differentiation. Moreover, some crucial gene expression patterns are differential AS events. Collectively, all these results provide a glimpse of dynamic involvements and systematical profiling of AS events in chicken ESCs, gPGCs, and SSCs, which might arouse the research use of numerous, potential biomarkers and AS events in the study of male germ cell differentiation.
A diverse splicing pattern in one gene leads to a variety of isoforms and expression patterns, which makes AS and its regulation mechanism more complex in the germ cellular process [47]. In this study, a total of 38,494 AS events out of 15,338 genes in ESCs, 48,955 AS events out of 14,783 genes in gPGCs, and 49,900 out of 15,089 genes in SSCs were detected. This is significant, in that AS is a commonly used process during male chicken germ cell differentiation. In addition, all cells generated the largest number of A5 and A3 , while the smallest number of AEs was generated. It is possible that the A5 /A3 is major factor in determining the transcript splicing, thus each alternative transcript is somewhat defined by the A5 /A3 . After comparing the difference of AS preference through ESCs, gPGCs, and SSCs, five types of AS events were increased, and ES, RI, and AE were all significantly enhanced in gPGCs and SSCs. The difference in AS has also been described in mammalian spermatogenesis, giving support to the agreement that changes in AS may contribute to male germ cell differentiation [48,49]. Moreover, the distribution of AS events demonstrated differential AS event patterns in ESCs, gPGCs, and SSCs. It is obvious that the splicing features of AS events may lead to direct differentiation.
To further understand the gene function affected by AS through influencing male chicken germ cell differentiation, we performed PPI, GO, and KEGG analyses of AS genes, and found that they were enriched in some biological processes, such as the developmental process, cellular process, and metabolic process. These results further suggest that AS events play an important role during male chicken germ cell differentiation. An interesting finding is that the metabolic process was not only enriched in GO terms, but also highly enriched in KEGG pathways. Similarly, many reports conclude that metabolic changes accompany the SSC differentiation [50,51]. We suspect that AS events induce the pathway transmission during male germ cell formation via metabolic transform. In addition, the calcium signing pathway, a key player in cell differentiation and development, was found to be regulated by some spliced factors in our study [52,53]. Being aware of many pathways, such as Wnt, MAPK, and mTOR, involved in the chicken germ cell formation process, we observed that abundant genes within these crucial pathways exhibited AS events. The implication is that AS events may affect germ cell differentiation via these pathways [54][55][56]. Moreover, germ cell-specific genes such as NANOG, POU5F3, LIN28B, BMP4, STRA8, and LHX9, which act as stage markers of germ cells, were developmentally regulated and stage-specific. In our results, the expression of these genes was associated with AS events transmitted through male chicken germ cell differentiation [39][40][41][42][43][44]. AS events may regulate the expression of these genes and subsequently impact the differentiation and developmental cell process. Taken together, AS is a complex and important regulatory mechanism, which may play a crucial role during male germ cell differentiation.

Conclusions
In summary, we have provided a comprehensive atlas on the genome-wide scale of the AS event landscape in male chicken germ-line cells (ESCs, gPGCs, and SSCs) and presented the distribution of these AS events. AS events have reportedly been involved in the germ cell differentiation process since they were first discovered. This helps us understand how AS events are associated with the male germ cell differentiation process, and to possibly improve treatment options for male infertility in animals and humans.