Small RNA deep sequencing revealed microRNAs’ involvement in modulating cellular senescence and immortalization state

Unlike rodent cells, spontaneous immortalization of avian cells and human cells is a very rare event. According to patent publications and current literature, there are no more than 4 spontaneously immortalized chicken embryo fibroblast (CEF) cell lines established up to date. One of those cell lines is ADOL (Avian Disease and Oncology Laboratory) ZS-1 cell line, which was established by continuous passaging of the CEFs derived from the specific pathogen free (SPF) 0.TVB*S1 (commonly known as rapid feathering susceptible or RFS) genetic line of chickens. The RFS genetic line of chickens was developed and has been maintained on the SPF chicken farm of USDA-ARS facility, ADOL, in East Lansing, Michigan, which is known as one of a few lines of chickens that are free of any known avian endogenous virus genes. To explore potential roles that epigenetic factors may play in modulating cellular senescence processes and spontaneous immortalization state, total RNAs extracted from samples of the RFS primary CEFs, RFS CEFs reached the 21st passage, and the ZS-1 cells were subjected to small RNA sequencing. Collectively, a total of 531 miRNAs was identified in the 3 types of samples. In contrast to the primary CEF samples, 50 miRNAs were identified with significantly differential expression only in the 21st passage samples; a different subset of 63 differentially expressed miRNAs was identified only in the ZS-1 samples; the majority of differentially expressed miRNAs identified in both the 21st passage CEF and the ZS-1 samples were more or less directionally consistent. Gene Ontology analysis results suggested that the epigenetic factor, miRNAs, plays a role in modulating the cellular senescence and spontaneous immortalization processes through various bioprocesses and key pathways including ErbB and MAPK signaling pathways. These findings provided the experimental and bioinformatic evidence for a better understanding on the epigenetic factor of miRNAs in association with cellular senescence and spontaneous immortalization process in avian cells.


INTRODUCTION
Cell lines play critical roles in biomedical research and in prevention, diagnosis, and even therapeutic treatment of diseases, including infectious diseases, of human and 1 immortalization and are free of endogenous retroviruses. Cells of such types are also critical and ideal for biomanufacturing vital bioreagents, such as antibodies and vaccines (Garber et al., 2009;Giotis et al., 2019;Lee et al., 2020Lee et al., , 2021Lin et al., 2020;Chungu et al., 2021;He et al., 2021).
There are only a few spontaneously immortalized chicken cell lines reported (Kool, 2019). The very first non-vial, non-viral protein, and non-chemically transformed avian immortalized cell line is known as the DF-1 cell line (Foster and Foster, 1997), which was derived by continuously passaging of chicken embryo fibroblastic (CEF) cells of the line-0 chickens, a EV free line of specific pathogen free (SPF) chickens that has been developed and maintained on the experimental farm of the USDA-ARS, Avian Disease and Oncology Laboratory (ADOL) at East Lansing, Michigan, since 1985 (Bacon et al., 2000).
The DF-1 line of cells has been widely used as substrates in virus propagation, recombinant protein expression, recombinant virus production, verification of gene expression and gene overexpression, modeling of gene functions and functions of foreign genes, and exploring host-pathogen interactions. DF-1 cells have also been used in CRISPR/Cas9 gene editing studies (Zhao et al., 2017;Sid and Schusser, 2018). In particular to study of host-pathogen interactions between avian cells and avian leukosis viruses (ALV), it was confirmed that the spontaneously immortalized DF-1 cells are more suitable than regular primary CEF (Maas et al., 2006).
Another EV free SPF line of chickens has been developed and maintained on the USDA-ARS, ADOL farm (Zhang et al., 2008). This genetic line of chickens was derived from a F 2 flock that was initiated from a cross between the line-0 and another ADOL line of chickens known as 0.44-VB*S1-EV21 (Bacon, Hunt and Cheng, 2000). This new addition of the SPF line of chickens was named as 0.TVB*S1 (commonly known as rapid feathering susceptible or RFS). The TVB refers to the tumor virus B locus of chicken genome, a gene that encodes the cellular receptors necessary to mediate successful infection of the subgroup B, D, and E ALV. There are 3 major alleles at this locus, TVB*S1, TVB*S3, and TVB*R, that have been commonly detected and described in chicken (Klucking et al., 2002;Zhang et al., 2007). Briefly, TVB*S1 encodes receptors capable of mediating subgroup B, D, and E ALV infection. This is the most susceptible TVB allele. TVB*S3 encodes receptors that permit subgroup B D ALV infection but block subgroup E ALV from entering an infection cycle (Adkins et al., 2001). TVB*R encodes a dysfunctional receptor that is incapable of mediating any B, D, or E subgroup ALV infection (Klucking, Adkins and Young, 2002). The line-0 is TVB*S3/S*3 homozygous, susceptible to subgroup B and D ALV but resistant to subgroup E ALV infection (C/E). The new EV free SPF line, 0. TVB*S1, is TVB*S1/S*1 homozygous, susceptible to all subgroup B, D, and E ALV infection (C/0) (Zhang, Bacon and Fadly, 2008).
Use of the primary CEF of the 0.TVB*S1 line of chickens, a spontaneously immortalized chicken cell line, ADOL ZS-1, was developed (Zhang, 2016). The primary CEFs were dissociated and grown in culture with continuous passage until senescence began obviously observable. The cells were then concentrated during the cell senescence phase to maintain a 30% to 69% culture confluence. Eventually, the remained subpopulation of cells continued proliferation over 48 passages in culture. The ADOL ZS-1 cell line is free of EV and susceptible to all subgroups of ALV, including the subgroup E. The cells of the ADOL ZS-1 cell line and its subclones thereof may be used for inter alia the production of viral agents, including recombinant viral agents, expression or overexpression of recombinant proteins, exploring or validating gene or foreign gene functions, host-virus interactions, and diagnostic assays. The ADOL ZS-1 cell line is available at ATCC (https://www.atcc.org/products/pta-121623).
The epigenetic factor, miRNAs, is reportedly a class of key modulators of cellular senescence, and cellular senescence is a key player in tumor suppression processes (Liu et al., 2012;Abdelmohsen and Gorospe, 2015). This study aimed to explore miRNA expression of the ADOL ZS-1 cell line, and to provide bioinformatic evidence to elucidate the epigenetic factor's roles that may be involved with cellular senescence and/or the immortalization process free of all known transformation agents.

Cell Samples and Total RNA Preparation
Three types of cell samples were used in this study, which were primary CEFs, 21st passage CEFs, and the ADOL ZS-1 line of cells. Primary CEFs were recovered from fertilized eggs of the 0.TVB*S1 line by aseptically removing the embryonic torso of the 10-day-old embryos, mincing the tissue, dissociating the cells in a trypsin-containing Gibco Dulbecco's Modified Eagle Medium (DMEM), and collecting and growing the dissociated cells in enriched DMEM culture medium. After incubation, cells of this primary culture were recovered by trypsinization. Three independent preparations of the primary culture were sampled for total RNA extractions. The cells of the primary culture were also plated in fresh culture medium containing 6% fetal calf serum for additional passages. Between the 10th and 25th passages, a working solution of Trypsin (0.25% Trypsin in PBS) was used to briefly loosen and to discard some of the senescent (aging) cells from the surface of the petri dish before completely dissociating the cells in the petri dish by trypsinization to passage the cells. Three random plates of cells of the 21st passage were sampled for total RNA extractions. After repeated passage, spontaneously immortalized cells were obtained from the mixed cell population in a pure form. These cells were recovered and designated ZS-1. Cells of the ZS-1 cell line are immortal and have been passaged over 40 passages. Three plates of cells of the ADOL ZS-1 line were sampled for total RNA extractions. The total RNA samples of the cell samples from the primary culture, the 21st passage, and the 48th passage ADOL ZS-1 cell line, which all shared the same genetic lineage of the 0. TVB*S1 line of chickens, were extracted with TRIzol reagent (Invitrogen, Carlsbad, CA) following the manufacturer's instructions.

Small RNA Sequencing
Small RNA libraries were prepared using the Illumina Small RNA Library Preparation Kit (Illumina, San Diego, CA, Cat. 63600). The completed libraries were quantified, and quality controlled (QC) using a combination of Qubit dsDNA High Sensitivity (HS), Caliper LabChipGX HS DNA and Kapa Illumina Library qPCR Quantification assays (Science, 2012;Mehta et al., 2018;Bronkhorst et al., 2020). The small RNA libraries were then sequenced on an Illumina HiSeq 2500 platform in a single-end 50 base format. All these procedures were carried out at and by the Genomics group of the Research Technology Support Facility on Michigan State University campus (Genomics -Research Technology Support Facility (msu.edu)).

Small RNA Sequence Reads Analysis for miRNAs and Target Gene Prediction
The small RNA_Seq reads that passed QC were analyzed following the similar procedures described earlier (Zhang et al., 2020a). Briefly, the pass-filter (PF) reads (the number of clusters that passed Illumina's "Chastity filter") were analyzed with miRDeep* software package (v3.8) using the default parameters (An et al., 2013) except the adapter sequence and chicken genome-specific index files. The adapter sequence used in this analysis was TGG AAT TCT CGG GTG CCA AGG AAC TCC AGT CAC (Illumina); and the chicken genome build index (build_bwt_idx) files were constructed based on the chromosome information of the galGal 5.0 genome assembly along with the latest version of "gga.gff3" file downloaded from the miRbase website (Kozomara et al., 2019).
Target genes of the differentially expressed miRNAs were predicted using the built-in target gene prediction function in miRDeep*, which employees the most commonly used target gene prediction tool, TargetScan, to predict target genes of known and novel miRNAs (Lewis et al., 2005).

Droplet Digital PCR Analysis to Validate the RNA Sequence Data
A small random sample of identified miRNAs was subjected to Droplet Digital PCR (ddPCR) analysis to validate the levels of small RNA_Seq reads. Primers were customarily designed for ddPCR analysis of each of the sampled miRNAs using mature sequences following the procedures described by Balcells et al. (2011). The cDNA samples used in ddPCR validation were reversely transcribed from individual total RNA samples using the iScript  and following the manufacturer's instruction (Bio-Rad Laboratories, Inc., Hercules, CA). A ddPCR reaction of 25 mL in final volume was initially prepared per miRNA per biological sample containing 2 mL of cDNA (60 ng/mL), 12.5 mL of EvaGreen Supermix (final concentration: 1X; Cat No. 1864034), 0.5 mL of each forward and reverse primers (200 nM; synthesized by Eurofins Genomics, Huntsville, AL), and 9.5 mL of nuclease-free water. Of which, 20 mL were loaded into one of 8 sample channels of a DG8 cartridge (Cat No. 1864008, Bio-Rad). Each oil well was loaded with 70 mL of droplet generating oil (Cat No. 1864006, Bio-Rad). The loaded DG8 cartridges were placed on a QX200 droplet generator (Bio-Rad) to generate the digital droplets. Forty mL of the generated droplet emulsion per sample were transferred to a well in a 96-well PCR plate followed by polymerase chain reaction with EvaGreen on a C1000 Thermal Cycler (Bio-Rad). The cycling conditions were 95°C for 5 min, followed by 40 cycles of 95°C for 15 s, 58°C for 60 s, and a final extension step of 90°C for 5 min. The droplets post PCR were read well by well on a QX200 droplet reader (Bio-Rad). PCRpositive and PCR-negative droplets in each of the wells were counted and analyzed with the QuantaSoft Software (Version 1.7, Bio-Rad). The primers designed for ddPCR validation of the selected miRNAs are given in Table 2.

Differentially Expressed miRNA Identification and GO Terms Enrichment Analysis
The number of reads per microRNA for each biological sample was counted using HTSeq (Anders et al., 2015). In each of the pairwise comparisons (between the primary and the 21st passage CEF; between the primary and the ZS-1 line of cells), differentially expressed micro-RNAs were identified by use of a custom R script encompassing the DESeq R package (2.1.0). A filter criterion of FDR < 0.05 and FC > 2 was enforced. For some of the differentially expressed miRNAs that ended up with a zero-statistic estimate for a normalized average TPM (baseMeanA or baseMeanB) in a contrast, an arbitrary small value of 5 was assigned to substitute the zero in order to computer a numeric fold change, and then a log 2 fold change value for easier comprehension of the estimates. To better understand the functional involvements of the identified microRNAs differentially expressed between the primary and the ZS-1 line of cells, the predicted target gene lists of the upregulated micro-RNAs and the downloaded miRNAs that were identified only in the ZS-1 cell samples in contrast to the primary CEF samples were subjected to GO Terms and pathway analysis independently using the g:Proflier2 (https:// biit.cs.ut.ee/gprofiler/gost) online tools with the following options: Organism: Gallus gallus; Statistical domain scope: All known genes; Significance threshold: Bonferroni correction; User threshold: 0.01 (Reimand et al., 2007).

RESULTS AND DISCUSSION
Small RNA Sequence Reads of the Total RNA Samples Small RNA sequencing generated an average of 49.98 million PF reads per type of the primary, 21st passage, and 48th passage (ZS-1 cell line) cell samples, with a range of 13.9 to 20.9 and 45.9 to 54.9 million reads per biological sample and per type of cell-samples, respectively (Table 1). The raw sequence datasets (accession numbers: SAMN11676622 -SAMN11676630) are available at the Sequence Read Archive (SRA) under the National Center for Biotechnology Information (NCBI) website with an assigned BioProject SRA accession: PRJNA543529 (https://www.ncbi.nlm.nih.gov/sra/ PRJNA543529).

Validation of the Small RNA_Seq Expression Data by ddPCR
A randomly selected small set of miRNAs was subjected to ddPCR analysis using customer designed primers ( Table 2). The absolute quantification of the randomly selected four miRNAs in the same sets of total RNA samples by ddPCR was analyzed against the corresponding normalized PF reads (TPM) of the small RNA sequence data with a bivariate model. The correlation coefficients between the two sets of expression data for the examined miRNAs ranged from r = 0.83 to r = 0.89 with a range of statistic P < 0.01 to P < 0.002 (Figure 1), which provided a highly positive support in validation of the small RNA sequence data generated in this study and the estimates of the miRNA expression derived from the small RNA sequence datasets.
MicroRNA Profiles of the Primary CEFs, 21st Passage CEFs, and the ADOL ZS-1 Line of Cells A total of 531 miRNAs was identified at least in 3 of the nine bio-samples of the primary, 21st CEFs, and the ZS-1 line of cells. Two hundred and twenty-nine of those were known miRNAs and the rest of 302 were novel miRNAs (Supplementary Table S1). Of which, 504, 459, and 422 miRNAs were identified in the primary, 21st passage, and ZS-1 cell samples, respectively (Supplementary Tables S2−S4). A set of 50 miRNAs was reportedly identified in the DF-1 cells by small RNA sequencing (Chen et al., 2019). Thirty-one of those were overlapped with miRNAs identified in this study (Supplementary  Table S4). However, those 31 miRNAs were not only identified in the ZS-1 samples, but also in the primary and the 21st passage CEF samples. Another study reported that gga-miR-375 was identified in the DF-1 cells (Zhang et al., 2020b), which was identified in the primary CEF and the ZS-1 samples (Supplementary Table S4), but not in the 21st passage CEF samples. By comparing the 3 lists of identified miRNAs of the 3 types of cell samples, a total of 42 miRNAs (18 known and 24 novel miRNAs) was found present only on the list of the primary CEF samples; 5 novel miRNAs were identified only in the 21st passage CEF samples; and a set of different 5 novel miRNAs was detected only in the ZS-1 cells. The identification of different unique subsets of miRNAs in each of the 3 types of cell samples indicated that miR-NAs and miRNA expression may have been involved with cell passaging and even were associated with cellular senescence and immortalization processes. Furthermore, there were 62, 25, and 17 miRNAs identified in both the primary and 21st CEFs, the primary CEFs and ZS-1 cells, and the 21st passage CEFs and the ZS-1 samples, respectively (Supplementary Table S5). It was noticeable that there were over 5 times of specific miR-NAs detected only in the primary CEF samples in contrast to both the 21st passage CEF and the ZS-1 cell samples. There were 5 miRNAs detected either only in the 21st passage CEF samples or the ZS-1 cell samples, and all were novel miRNAs. Based on the reported observations during the development of the very first spontaneously immortalized chicken embryo fibroblastic cell line, the DF-1, and subsequent studies for characteristics of the DF-1 cells, normal primary CEFs entered a highly senescent phase post 10 passages with a slow growth rate in tissue culture; by the 40th passage, a survived homogeneous population of cells grew rapidly (Himly et al., 1998;Lassiter et al., 2014). Therefore, it is speculated that the unique subsets of miRNAs that detected only in the primary CEFs, 21st passage CEFs, and the ZS-1 cells may be more or less associated with the maintenance of normal cell growth, reduced normal cell growth and cell senescence, and spontaneous transformation in addition to other factors and other miR-NAs that were detected in more than one type of cell samples but with differential expressions between the different types of cell samples.

Differentially Expressed miRNAs of the ADOL ZS-1 Cells and the 21st Passage Chicken Embryo Fibroblasts in Contrast to the Primary Chicken Embryo Fibroblasts
A total of 188 (96 known and 92 novel) miRNAs was identified with significantly differential expression (2.30E-194 ≤ P ≤ 2.10E-02 and 6.06E-192 ≤ FDR ≤ 4.50E-02) between the ADOL ZS-1 line of cells and the primary CEFs. Fifty-five of those miRNAs were upregulated and the rest of 133 were downregulated in expression in the ADOL ZS-1 line of cells (Table 3). There were 175 (100 known and 75 novel) miRNAs identified with significantly differential expression (3.59E-268 ≤ P ≤ 2.05E-02 and 9.45E-266 ≤ FDR ≤ 4.95E-02) between the 21st passage and the primary CEFs. Sixty-two of those miRNAs were upregulated and the rest of 113 miRNAs were downregulated when the CEFs reached the 21st passage (Table 4). With a close look at the 2 lists of differentially expressed miRNAs derived from the 2 contrasts, between the ZS-1 cells and the primary CEFs and between the 21st passage and the primary CEFs, eighteen miRNAs were upregulated and 45 miR-NAs were downregulated in expression only in the ZS-1 cells; 24 miRNAs were upregulated and 26 miRNAs were downregulated in expression only in the 21st passage CEFs; the rest of differentially expressed miRNAs were all directionally consistent more or less upregulated or downregulated in expression between the 2 contrasts (Supplementary Table S6). Increasing evidence indicates that miRNAs play vital roles in a variety of biological processes including maintaining cellular mechanism and proliferation (Yin et al., 2021). The miR-21 is reported as an oncogene and plays a key role in resisting programmed cell death in cancer and targets apoptosis (Buscaglia and Li, 2011) The gga-mir-21 was identified in all 3 types of samples in this study, and it was overexpressed significantly in both the 21st passage CEF samples and the ZS-1 samples in contrast to the primary CEF samples (Supplementary Table S6), suggesting that gga-mir-21 may have exerted influence in overcoming cellular senescence in the 21st passage CEFs and in maintaining the immortalization status of the ZS-1 cells. Other differentially expressed miRNAs, including novel microRNAs, identified in the contrast between the ZS-1 samples and the primary CEF samples might play role (s) similar to the gga-mir-21, which remain to be further investigated in the future. Figure 1. A bivariate plot illustrating the moderately high associations between the droplet digital PCR generated absolute quantifications (ddPCR) and the normalized small RNA sequence data (TPM) of a set of randomly selected miRNAs. This result supports the validity of the small RNA sequence data generated in this study.    Predicted Target Genes of the Differentially Expressed miRNAs Observed in the Contrast between the ZS-1 Cells and the Primary CEFs A total of 11,499 target genes was predicted for the 188 differentially expressed miRNAs identified in the contrast between the ADOL ZS-1 cell samples and the primary CEFs (Supplementary Table S7). A further close look of the target gene list led to 2 short lists of target genes, one was a list of 5,436 target genes predicted for 18 upregulated miRNAs in expression, which were only identified in the ZS-1 cell samples in contrast to the primary CEF samples (Supplementary Tables S6, S8); the other one was a list of 4,757 target genes predicted for 45 downregulated miRNAs in expression, which were only identified in the ZS-1 cell samples in contrast to the primary CEF samples (Supplementary Tables S6, S9). To better understand how, if any, the miRNAs would be involved in overcoming the process of cell senescence and/or maintaining a stable spontaneous immortalization state, the latter two relatively short lists of target genes were independently subjected to GO Terms analysis to explore potential roles through which the differentially expressed miRNAs may play in modulating cellular senescence and immortalization processes.

Target-Gene-Set Enrichment in GO Terms and Pathways
Two sets of target genes of the upregulated and downregulated miRNAs identified only in the ZS-1 cell samples in contrast to the primary CEF samples were independently subjected to GO Terms enrichment analyses and the outputs of the analyses showed that both sets of the target genes were highly enriched in a variety of GO Terms and pathways (Supplementary Tables S10 −S11). The target genes of the upregulated miRNAs were enriched in over thirty molecular function (MF), 353 biological process (BP), and 48 cellular component (CC) GO Terms, in addition to KEGG and Reactome (REAC) pathways, and 57 transcription factors (Figure 2). The target genes of the downregulated miR-NAs were enriched in 24 MF, 330 BP, and 30 CC GO Terms, in addition to KEGG (MAPK, ErbB, and Endocytosis) and Reactome pathways, and 60 transcription factors ( Figure 3). Studies in human diploid somatic cell population doublings revealed that the cells enter a status known as terminal proliferation arrest state of senescence due to an intrinsic mechanism involving p53 and pRB/p16 INK4 mediated pathways. Telomere shortening is reportedly a major factor attributing to the proliferation arrest state of senescence. No GO Terms nor Figure 3. A Manhattan plot illustrating GO Terms enrichments of target genes of the downregulated miRNAs. These differentially expressed miRNAs were detected only in the ZS-1 cell samples in contrast to the primary chicken embryo fibroblast samples. The enrichments of the GO Terms MF, BP, CC, and the KEGG pathways, REAC pathways, WP pathways, as well as TF and MIRNA are depicted, respectively, in the plot, which also suggested strong and wide involvements of the downregulated set of miRNAs in biological processes through their set of target genes in the ZS-1 cells. Figure 2. A Manhattan plot illustrating GO Terms enrichments of target genes of the upregulated miRNAs. These differentially expressed miR-NAs were detected only in the ZS-1 cell samples in contrast to the primary chicken embryo fibroblast samples. The enrichment of the GO Terms MF (molecular function), BP (biological processes), CC (cellular component), and the KEGG (Kyoto Encyclopedia of Genes and Genomics) pathways, REAC (reactome) pathways, WP (WiKi) pathways, as well as TF (transcription factor) and MIRNA (microRNA target base) is depicted, respectively, in the plot, suggesting strong and wide involvements of the upregulated set of miRNAs in biological processes in the ZS-1 cells. pathways (Supplementary Tables S10−S11) identified for the target gene lists of the upregulated and downregulated miRNAs in the ZS-1 cell samples in contrast to the primary CEF samples were involved with the p53 or the pRB/p16 INK4 . Other evidence, however, suggests that some genes are involved in immortalization but have nothing to do with the telomere status (Duncan et al., 2000). The target gene list of the overexpressed miRNAs (Supplementary Table S8) is associated with the MAPK signaling pathway (Supplementary Table  S10). MAPK (mitogen-activated protein kinase) signaling pathway is a highly conserved module that is involved in various cellular functions, including cell differentiation, migration, and proliferation (Zhang and Liu, 2002;Mebratu and Tesfaigzi, 2009;Ables et al., 2011). The ErbB family of proteins contains 4 receptor tyrosine kinase receptors, which include EGFR/ ERBB1/HER1, ERBB2/HER2, ERBB3/HER3, and EFBB4.HER4. Transphosphorylation of the ErbB dimer partners by activation of HER2 and EGFR stimulate intracellular pathways and STAT transcription factors. Of which, the EGFR and HER2 can activate STAT3 via phosphorylation through SRC to promote cell survival (Kavarthapu et al., 2021).

CONCLUSIONS
In summary, a total of 422 miRNAs was identified in samples of the ADOL ZS-1 cell line and 188 of those miRNAs were differentially expressed between the ZS-1 and the primary chicken embryo fibroblast (CEF) samples. Bioinformatic analysis of the target gene list of the differentially expressed miRNAs that were only identified in the ZS-1 samples (not in the 21st passage CEFs) in contrast to the primary CEF samples suggested that those miRNAs might have modulated the cell immortalization processes and contributed to the maintenance of the immortalization state by targeting genes involved in various biological processes and key pathways. This study provided a piece of experimental and informatics evidence that the epigenetic factor, miRNAs, plays a role in modulating cell proliferation, senescence process and immortalization state.