Transcriptome Sequencing and Comparative Analysis of Amphoteric ESCs and PGCs in Chicken (Gallus gallus)

Simple Summary The study of chicken embryonic stem cells (ESCs) and primordial germ cells (PGCs) showed the potential application of developmental biology and translational medicine. However, the difference between amphoteric ESCs and PGCs is still elusive, limiting the accuracy of correlative research. In this paper, chicken amphoteric ESCs and PGCs were isolated, separated, and sequenced to explore their dynamic transcriptomes. Our results provide a knowledge base of transcriptomes in chicken amphoteric ESCs and PGCs, which will help other researchers interested in studying relative biological processes. Abstract Chicken (Gallus gallus) pluripotent embryonic stem cells (ESCs) and primordial germ cells (PGCs) can be broadly applied in the research of developmental and embryonic biology, but the difference between amphoteric ESCs and PGCs is still elusive. This study determined the sex of collected samples by identifying specific sex markers via polymerase chain reaction (PCR) and fluorescence activated cell sorter (FACS). RNA-seq was utilized to investigate the transcriptomic profile of amphoteric ESCs and PGCs in chicken. The results showed no significant differentially expressed genes (DEGs) in amphoteric ESCs and 227 DEGs exhibited in amphoteric PGCs. Moreover, those 227 DEGs were mainly enriched in 17 gene ontology (GO) terms and 27 pathways according to Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. Furthermore, qRT-PCR was performed to verify RNA-seq results, and the results demonstrated that Notch1 was highly expressed in male PGCs. In summary, our results provided a knowledge base of chicken amphoteric ESCs and PGCs, which is helpful for future research in relevant biological processes.


Introduction
Chicken (Gallus gallus) is a classical model in developmental and embryonic biology [1]. For decades, chickens have played a vital role in animal research as alternatives and outbred

RNA Extraction and Sequencing
RNA was extracted with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. To effectively remove the genomic DNA, we added RNase-free DNase I (Takara, Dalian, China) to the reaction mixture for at least 10 min. The extracted RNA was quantified with the Nanodrop system (Thermo, Wilmington, DE, USA), and the fragment size distribution was checked by 1.5% agarose gel electrophoresis. RNA concentration was assessed using a Thermo NanoDrop2000TM spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA), and the RNA integrity number (RIN) was determined using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) ( Table 1). The extracted RNA was stored for reverse-transcription and sequencing. The RNA Libraries pools of four kinds of cells were established following the protocol of Illumina mRNA-seq with 50 ng RNA, and the experiments were performed by the Shanghai OE Biotech Co., Ltd. (Shanghai, China).

Differential Expression Analysis and Functional Enrichment
Filtering and quality control checks of raw reads from RNA-seq were done by FastQC [22]. The clean reads were mapped to reference sequences using SOAP2 [23]. The gene expression levels were quantified by calculating the RPKM (Reads Per kb transcriptome per Million reads) values. The gene expression pattern and PCA analyses were performed on R [24]. The expression of genes with a fold change > 2 and FDR < 0.001 were filtered as differentially expressed genes (DEGs). The functional analysis of DEGs was carried out by GO analysis and KEGG pathway enrichment analysis using DAVID [25].

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)
Tissues and cells were homogenized in TRIzol Reagent, and total RNA was isolated according to the manufacturer's instruction (QIAGEN, Beijing, China, DP424). qRT-PCR was performed using the FastKing One-Step RT-PCR Kit with SYBR green (QIAGEN, Beijing, China, KR123). mRNA levels of related genes were determined by CFX-Connect Real-Time PCR detection system (BIO-RAD, California, USA, 7500fast). Results were represented as quantification normalized to relative housekeeping genes (β-Actin) using the 2-∆∆Ct method. The sequences of qRT-PCR primers are listed in Table 2.

Amphoteric ESCs and PGCs Confirmation
For sample preparation for transcriptome analyses, the blastoderms and gonads were collected from incubated chicken eggs at day 0 and day 4.5 ( Figure 1a). The sex of collected tissues was identified by PCR as described previously. PCR products with two bands (600 and 450 bp) were observed in tissues collected from females (ZW) and only one band (600 bp) in males (ZZ). In total, ESCs were isolated from 9699 eggs, including 4845 females and 4854 males; PGCs were obtained from 3150 eggs, including 1594 males and 1556 females. Samples of same-sex eggs from each stage were collected for isolation and culturing. The result showed that no morphological differences between the cultured amphoteric ESCs and PGCs. To further confirm the cell type, different cell surface markers

Amphoteric ESCs and PGCs Confirmation
For sample preparation for transcriptome analyses, the blastoderms and gonads were collected from incubated chicken eggs at day 0 and day 4.5 ( Figure 1a). The sex of collected tissues was identified by PCR as described previously. PCR products with two bands (600 and 450 bp) were observed in tissues collected from females (ZW) and only one band (600 bp) in males (ZZ). In total, ESCs were isolated from 9699 eggs, including 4845 females and 4854 males; PGCs were obtained from 3150 eggs, including 1594 males and 1556 females. Samples of same-sex eggs from each stage were collected for isolation and culturing. The result showed that no morphological differences between the cultured amphoteric ESCs and PGCs. To further confirm the cell type, different cell surface markers

Transcriptome Characterization of Amphoteric ESCs and PGCs in Chickens
Eight sequencing libraries were constructed for the Illumina platform with two biological repeats to obtain the genome-wide gene expression profile of amphoteric ESCs and PGCs. In total, 36,221 transcripts could be detected in combined transcriptomes, and gene expression was quantified by calculating RPKM values. A heatmap illustrating the hierarchical clustering of RPKM values was generated to visualize the overall gene expression pattern (Figure 2a). The PCA analysis was utilized to explore gene expression patterns among amphoteric ESCs and PGCs, revealing different patterns between ESCs and PGCs. Moreover, it showed that amphoteric ESCs clustered together closely, while amphoteric PGCs were relatively dispersed, indicating amphoteric ESCs exhibit similar patterns, but PGCs exhibit different patterns (Figure 2b). Differential analysis of transcriptome data demonstrated no DEGs in amphoteric ESCs, while 227 DEGs were found in amphoteric PGCs (Figure 2c and Figure S1). Among 227 DEGs, 52 DEGs were highly expressed in female PGCs and 172 in male PGCs (Table S1). Furthermore, no known or putative regulator of gonadal sexual determination and differentiation was identified among 227 DEGs, which implied that major components of sex determination and differentiation pathways are activated after this stage, but the amphoteric PGCs initiated the differentiated expression.
Animals 2020, 10, x 6 of 12 morphology (right medium, the red arrow points to typical ESCs clones, and the black arrow points to emblematic PGCs mass), and fluorescence activated cell sorter (FACS) (right bottom).

Transcriptome Characterization of Amphoteric ESCs and PGCs in Chickens
Eight sequencing libraries were constructed for the Illumina platform with two biological repeats to obtain the genome-wide gene expression profile of amphoteric ESCs and PGCs. In total, 36,221 transcripts could be detected in combined transcriptomes, and gene expression was quantified by calculating RPKM values. A heatmap illustrating the hierarchical clustering of RPKM values was generated to visualize the overall gene expression pattern (Figure 2a). The PCA analysis was utilized to explore gene expression patterns among amphoteric ESCs and PGCs, revealing different patterns between ESCs and PGCs. Moreover, it showed that amphoteric ESCs clustered together closely, while amphoteric PGCs were relatively dispersed, indicating amphoteric ESCs exhibit similar patterns, but PGCs exhibit different patterns (Figure 2b). Differential analysis of transcriptome data demonstrated no DEGs in amphoteric ESCs, while 227 DEGs were found in amphoteric PGCs (Figure 2c and Figure  S1). Among 227 DEGs, 52 DEGs were highly expressed in female PGCs and 172 in male PGCs (Table  S1). Furthermore, no known or putative regulator of gonadal sexual determination and differentiation was identified among 227 DEGs, which implied that major components of sex determination and differentiation pathways are activated after this stage, but the amphoteric PGCs initiated the differentiated expression.

Enrichment Analysis of DEGs in Amphoteric PGCs in Chickens
The gene ontology database [26] was used to perform functional annotation on DEGs transcripts' three components, including biological processes (BP), molecular function (MF), and cellular component (CC). The clusterProfiler R package was used to analyze GO enrichment of DEGs transcripts. GO terms (corrected p-value < 0.05) were considered significantly enriched. The 17 significantly enriched GO terms were annotated as biological processes (9; 52.9%), cellular component (3; 17.64%), and molecular function (5; 29.41%), as shown in Figure 3a and Table S2. Interestingly, almost all GO terms in biological processes (8/9) were sub-terms of the developmental process, which indicated the developmental-difference process is active in amphoteric PGCs ( Figure  3b). Moreover, we concluded that Notch1 is a key gene in developmental process terms according to all genes' interaction prediction analysis (Figure 3c). Pathway significant enrichment analysis was performed to identify host genes involved in major biochemical metabolic pathways and signal

Enrichment Analysis of DEGs in Amphoteric PGCs in Chickens
The gene ontology database [26] was used to perform functional annotation on DEGs transcripts' three components, including biological processes (BP), molecular function (MF), and cellular component (CC). The clusterProfiler R package was used to analyze GO enrichment of DEGs transcripts. GO terms (corrected p-value < 0.05) were considered significantly enriched. The 17 significantly enriched GO terms were annotated as biological processes (9; 52.9%), cellular component (3; 17.64%), and molecular function (5; 29.41%), as shown in Figure 3a and Table S2. Interestingly, almost all GO terms in biological processes (8/9) were sub-terms of the developmental process, which indicated the developmental-difference process is active in amphoteric PGCs (Figure 3b). Moreover, we concluded that Notch1 is a key gene in developmental process terms according to all genes' interaction prediction analysis (Figure 3c). Pathway significant enrichment analysis was performed to identify host genes involved in major biochemical metabolic pathways and signal transduction pathways in KEGG database [27]. The clusterProfiler R package was used to analyze the statistical enrichment of genes differentially expressed in KEGG Animals 2020, 10, 2228 7 of 11 pathways. Specifically, 24 KEGG terms were significantly enriched in chicken amphoteric PGCs (Figure 3d and Table S3). In addition, the Notch pathway was also significantly enriched in amphoteric PGCs, which is consistent with the interaction prediction analysis results.

Validation of Transcriptome Data by Real-Time qRT-PCR and Notch1 Specific Expression in Male PGCs
A total of 17 genes were selected to validate the result obtained from RNA-Seq data (Figure 4a). Consistent qPCR results indicated the reliability of RNA-Seq data and quantified gene expression accuracy in amphoteric PGCs.   High expression of Notch1 was found in male PGCs at both mRNA and protein levels, which demonstrated the Notch1 is a differential expression factor in amphoteric PGCs (Figure 4b).

Discussion
In this study, we conducted the transcriptomic gene expression analysis of amphoteric ESCs and PGCs in chickens. No genes were differentially expressed in amphoteric ESCs; however, differentially expressed genes were identified in amphoteric PGCs and also found to be highly enriched in significant GO terms and KEGG pathways. Moreover, we confirmed that Notch1 is highly expressed in male PGCs.
Genetic sex determination in mammals and birds occurs at fertilization with the differential inheritance of sex chromosomes. However, sex determination and differentiation are complicated processes that occur during embryonic development [28,29]. As omnipotent or pluripotent cells, ESCs, especially chicken ESCs, have been reported to resemble mice ESCs, which retain the stage of naive and primed [30][31][32]. Regarding mice, expression patterns in males and females vary in germ cells, which indicated that sexual dimorphisms occur in mouse ESCs [33]. However, our finding showed that chicken amphoteric ESCs share a great number of genes and have similar expression patterns, which indicates the homogeneity and highly undifferentiated stages of chicken amphoteric ESCs. Although the sex of both chickens and mice is controlled by sex chromosomes (W/Z or X/Y system), the sex determination in chickens is still elusive because of the absence of X chromosome inactivation [34][35][36]. We propose that might exhibit a potential specific mechanism cause the different expression patterns in amphoteric ESCs between chickens and mice.
Genital ridges are precursor organs of gonads (ovary and testis) containing PGCs and somatic cells [37][38][39]. In gonads, the transition of germ cells from PGCs to sperm cells or eggs is not determined by the genotype but by the somatic cell type of gonads [40,41]. The RNA-seq data revealed sexually dimorphic gene expression in gonad (genital ridge) tissues (day 4.5) before gonadal differentiation (day 6), including some sexual differentiation genes such as FOXL2 [19]. In this paper, the amphoteric PGCs initiate dimorphic gene expression, but these DEGs not including known or putative sex genes. Compared with gonad tissues at day 4.5, amphoteric PGCs represent lower heterogeneity, suggesting the uninitiated germ cell sex determination. A recent study has demonstrated that the amphoteric chicken PGCs could be induced by spermatogonial stem cells (SSCs) or oocytes via co-culture with sertoli cells or granulosa [42][43][44], which showed that PGC sex differentiation could be activated by somatic cells; this may explain different expression patterns of amphoteric gonad tissues with PGCs.
Furthermore, according to our gene ontology analysis, many DEGs are involved in the developmental process and KEGG pathways showed much cell-to-cell interaction and restructuring of the gonad into the morphologically distinct cytokine-cytokine receptor interaction, Toll-like receptor signaling pathway, focal adhesion, and ECM-receptor interaction. Chicken gonadal sexually dimorphic gene expression has been reported to be involved in ECM and cytoskeleton formation and regulation to facilitate testis or ovary formation [20]. Therefore, we conclude that amphoteric PGCs could interact with gonad somatic cells to drive the amphoteric morphological change in E4.5. Moreover, we identified the notch signal key factor, Notch1, is highly expressed in male PGCs and plays a critical role in the interaction network of DGEs. It is well known that Notch1 and notch signal are essential in cell-cell communication and further regulate embryonic development. In mice, Notch1 was reported to induce Sox9 expression during the early stages of embryo differentiation, and Sox9 is a crucial regulator through the process of masculinity [45]. Similarly, constitutive activation of Notch1 signaling in sertoli cells causes the exit of gonocytes from quiescence, indicating Notch1 promotes masculinity during sex differentiation [46]. In chickens, our previous study showed that the Notch1 and notch signal plays a crucial role in chicken PGC formation and production of spermatogonium [47]. Thus, our future studies will be focusing on investigating the effect of amphoteric PGCs and gonad somatic cells during sex differentiation and the functional role of Notch1 in amphoteric PGCs.

Conclusions
Chicken amphoteric ESCs and PGCs were isolated, separated, and sequenced to explore their dynamic transcriptomes. According to our results, amphoteric ESCs share many genes and exhibit similar expression patterns, implying that chicken amphoteric ESCs are homogenous and still in highly undifferentiated stages. In contrast, chicken amphoteric PGCs are heterogeneous and exhibit differentially expressed genes enriched in GO terms and the KEGG pathway involved in developmental processes. High expression of Notch1 in male PGCs might play a key role in amphoteric PGC differentiation. Overall, our results provide a knowledge base of transcriptomes in chicken amphoteric ESCs and PGCs, which would help other researchers interested in relevant biological processes.