Transcriptome Analysis of Long Noncoding RNAs and mRNAs in Granulosa Cells of Jinghai Yellow Chickens Illuminated With Red Light

Jinghai Yellow chickens are a new indigenous breed with a dual purpose in China, but their egg laying performance is limited. Compared with white light (WL), exposure to red light (RL) can improve the egg laying performance of hens. Herein, to elucidate the molecular mechanism by which RL affects the egg laying performance, RNA sequencing was used to analyze long noncoding RNAs (lncRNAs) and mRNAs from granulosa cells of small yellow follicles from Jinghai Yellow chickens in RL and WL groups. A total of 12,466 lncRNAs were identified among the assembled transcripts, of which 168 lncRNAs were significantly different between the RL and WL groups (101 downregulated and 67 upregulated). Additionally, 1182 differentially expressed mRNAs were identified (958 downregulated and 224 upregulated). Integrated network analysis demonstrated that numerous differential mRNAs were involved in follicular development through steroid hormone synthesis, oocyte meiosis, and the PI3K-Akt signaling pathway. The impact of lncRNAs on cis and trans target mRNAs indicates that some lncRNAs play important roles in follicular development of small yellow follicles. The results provide a starting point for studies aimed at understanding the molecular mechanisms by which monochromatic light affects follicular development and egg production in hens.


INTRODUCTION
Artificial illumination is widely used to promote reproductive performance in birds. Casey et al. (1969) concluded that shorter wavelengths stimulated reproductive development of pullets, but short wavelength illumination has a minimal effect on the rate of egg production in laying hens and breeder turkeys (Schumaier et al., 1968;Wells, 1971;Jones et al., 1982). Red light (RL) increases egg production in brown egg laying hens (Foss and White, 1983), turkeys (Pyrzak and Siopes, 1986), quails (Woodard et al., 1969), and pigeons (Wang et al., 2018). Furthermore, Pyrzak et al. (1984) suggested that hens produce fewer but larger eggs than earlier maturing hens exposed to red and white light illumination. However, Gongruttananun (2011) indicated that RL has no effect on egg production in hens, and Li et al. (2014) concluded that the wavelength of light does not affect egg production in hens. Thus, research on the effects in poultry is complicated, and few studies have focused on broiler chickens. Five or six hierarchical follicles simultaneously reside in the ovary to support daily ovulation in hens. One single follicle from a cohort of 8-13 small yellow follicles (SYFs) 6-8 mm in diameter is selected each day to enter the preovulatory hierarchy following ovulation of the largest follicle (Johnson and Woods, 2009;Onagbesan et al., 2009). Follicles consist of an outer layer of theca inerna cells that encircle inner layers of granulosa cells (GCs) (Hsueh et al., 1984). Follicular selection is associated with GCs (Hernandez and Bahr, 2003), during which GCs become capable of producing progesterone, resulting in a surge in luteinizing hormone (Robinson and Etches, 1986). Onagbesan et al. (2009) concluded that although the oocyte is the major source of epidermal growth factor-like peptides, the granulosa and theca are the sites of action of the ligands. Elucidating the molecular mechanisms mediating GC development in hens under monochromatic light could help to improve the reproductive performance of broilers.
Long noncoding RNAs (lncRNAs) are non-translated RNAs longer than 200 bp that play an important role in preand post-translational mRNA processing (Cheetham et al., 2013). Increasing evidence indicates that lncRNA-mediated gene expression is critical in reproduction in both male and female animals (Liu et al., 2018;La et al., 2020). Peng et al. (2019) identified 160 mRNAs and 550 lncRNAs that differ in follicles between two different chicken breeds, many involved in oocyte meiosis, progesterone-mediated oocyte maturation, and cell cycle pathways. Meanwhile, Ren et al. (2017) revealed that 52 lncRNAs were closely correlated with divergent reproductive mRNAs in the different phases of duck ovaries. The functions of lncRNAs are closely related to the development of follicles, but data associated with the functions of lncRNAs in follicle development in chicken under monochromatic light remains limited.
In the current study, Illumina sequencing technology was employed to identify lncRNAs and mRNAs in the GCs of SYFs under red light that are related to follicle development. The results could prove useful for exploring the molecular mechanisms mediating the development of GCs under monochromatic light, and help to improve the egg laying performance of broilers.

Chicken Rearing and Sample Preparation
Jinghai Yellow chickens were raised by Jiangsu Jinghai Poultry Industry Group Co., Ltd. (Nantong, Jiangsu, China). After transfer to the laying house, hens were caged individually, and 300 hens were divided into RL and WL groups, with five subgroups in each. All hens were provided with water ad libitum and restricted food. Experimental birds were exposed to red light (RL, 660 nm) while control birds were exposed to white light (WL, 400 to 760 nm) using light-emitting diodes (Shenzhen Hongda Technology Co., Ltd, Shenzhen, China) for 16 h each day (16 h light, 8 h dark). The light intensity was 15.0 lux as measured with a TES-1336A light meter (TES Electrical Electronic Crop., Taipei, China). Egg production and egg weight were measured, and based on the pedigree record, six hens at 300 days of age with an average body weigh were selected. All efforts were made to minimize distress. SYFs with a diameter of 6-8 mm were washed carefully in cold phosphate-buffered saline (PBS; Gibco, USA) and collected using the method for collecting GCs reported previously by Gilbert et al. (1977), flash-frozen in liquid nitrogen, and stored at −80 • C.

RNA Sequencing (RNA-Seq) Sample Preparation and Sequencing
Total RNA from each sample was isolated using TRIzol reagent (Invitrogen, USA). An RNA Nano 6000 Assay Kit and a Bioanalyzer 2100 system (Agilent Technologies, USA) were employed to determine the integrity of RNA, a Nanodrop instrument (Thermo Fisher Scientific, USA) was used to assess the purity and quantity of RNA, and the RIN ranged from 9.0 to 10.0. Six lncRNA libraries were constructed from SYFs of hens raised under RL (R1, R2, R3) or WL (W1, W2, W3). A total of 3 µg RNA from each sample was used as input material for RNA sample preparation. An Epicentre Ribo-zero rRNA Removal Kit (Epicentre, USA) was used to remove ribosomal RNA, and ethanol precipitation was applied to clean up the rRNA-free samples. An NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB, USA) was employed to generate sequencing libraries. Random hexamer primer and M-MuLV Reverse Transcriptase were used to synthesize first-strand cDNA, DNA Polymerase I and RNase H were applied to synthesize second-strand cDNA, and dNTPs with dTTP were replaced by dUTP. The 3 ′ ends of DNA were adenylated, and NEBNext Adaptor with a hairpin loop structure was ligated to prepare for hybridization. An AMPure XP system (Beckman Coulter, USA) was employed to purify the library fragments. The quality of the library was then measured by an Agilent Bioanalyzer 2100 System. Paired-end reads were sequenced on an Illumina Hiseq 4000 platform (30×) at Shanghai Oebiotech Co., Ltd (Shanghai, China).

Co-expression (trans) and Co-location (cis) Analyses
Pairwise significantly differentially expressed mRNAs and lncRNAs were estimated using the Pearson's correlation coefficient (r). The small number of samples (three each from WL and RL groups) used for the co-expression analysis is a limitation of our study. Those mRNAs with a p-value < 0.01 and |r-value| > 0.9 were considered co-expressed genes of their respective lncRNAs. FEELnc (v 0.1.1) (Wucher et al., 2017) was employed to screen the coding genes located within 100 kb upstream and downstream of lncRNAs for potential cis-regulation. RIsearch-2.0 software (Alkan et al., 2017) was used to identify target genes in trans, with parameters set as the base number of direct interactions between lncRNAs and mRNAs ≥10 and free energy ≤-100. Differentially expressed lncRNAs and their corresponding differentially expressed cisand trans-target genes were used to construct lncRNA-gene interaction networks using Cytoscape v3.2.1 (Smoot et al., 2011).

Target Gene Prediction
To gain further insight into the functions and classifications of the identified lncRNA targets, we performed Gene Ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotation of predicted lncRNA targets using the DAVID gene annotation tool (http://david.abcc.ncifcrf.gov/). We used KOBAS software to test the statistical enrichment of differentially expressed genes and lncRNA target genes in KEGG pathways (Peng et al., 2019).

Real-Time Quantitative PCR (RT-qPCR) Analysis
Samples were isolated from GCs of SYFs and RT-qPCR was used to validate DE lncRNAs and mRNAs identified by RNA-Seq. RT-qPCR was performed using a LightCycler 480 II Real-time PCR Instrument (Roche, Swiss) with ChamQ SYBR qPCR Master Mix (Vazyme, China). Each 10 µl PCR mixture contained 1 µl of cDNA, 5 µl of 2× ChamQ SYBR qPCR Master Mix, 0.2 µl of forward primer, 0.2 µl of reverse primer, and 3.6 µl of nucleasefree water. Reactions were incubated in a 384-well optical plate (Roche, Switzerland) at 95 • C for 30 s, followed by 40 cycles at 95 • C for 10 s, and 60 • C for 30 s. Each sample was run in triplicate for analysis. At the end of each PCR cycle, melting curve analysis was performed to validate the specific generation of the expected PCR product. Specific primers for mRNAs and lncRNAs are listed in Supplementary Table 1. Using ACTB as a reference, relative expression levels of mRNAs and lncRNAs were quantified using the 2 − CT method (Livak and Schmittgen, 2001).

Statistical Analysis
Data are expressed as mean ± standard error, and one-way analysis of variance was performed with SPSS 13.0 software (SPSS Inc., Chicago, IL, USA). The statistical significance of differences among the various groups was evaluated by least significant difference post-hoc multiple comparisons tests, and p < 0.05 was considered statistically significant.

Reproductive Performance of Jinghai Yellow Chickens
Egg production by hens aged 300 days in the RL and WL groups is shown in Table 1. Egg production in the RL group was significantly greater than that in the control group (p = 0.023), but there was not statistically significant difference in egg weight between groups (p = 0.667). Although the body weight of hens was slightly higher in the RL group (p = 0.925), the difference was not statistically significant. These results suggest that RL enhanced the reproductive performance of hens.

Identification of lncRNAs and mRNAs in Hen Ovaries
Six cDNA libraries were built from the RL (n = 3) and WL (n = 3) groups to identify lncRNAs and mRNAs expressed in GCs of SYFs. We obtained 97.97-99.10 million raw reads after filtering out contaminated reads, low-quality reads, and those with unknown bases accounting for >5% of reads, resulting in 90.45-95.06 million clean reads (Supplementary Table 2). Next, 87.66-91.81% of clean reads from each library were mapped to the chicken reference genome. The average GC content was 47.81%, and Circos analysis showed that lncRNAs in GCs were distributed on almost all chromosomes, with the fewest on chromosome 32 and the most on chromosome 1 (Figure 1). A stringent filtering pipeline was applied to discard transcripts lacking all lncRNA characteristics, transcripts < 200 bp in length, and those with only two exons and three reads of coverage. The lncRNA genes had an average length of 1,408 bp and 2.5 exons. A total of 12,466 lncRNAs were included in the assembled transcripts, comprising 10,969 and 1,497 known and unknown lncRNAs (Supplementary Table 3). The majority of lncRNAs were from the genic intronic region (Supplementary Table 3). Expression levels, transcript lengths, and the number of exons between lncRNAs and mRNAs generated from six individual chicken samples are shown in Figure 2. The length of mRNA transcripts was greater than the length of lncRNAs, and most mRNAs included more than 20 exons, compared with only two or three exons in most lncRNAs. Furthermore, the average expression level measured for lncRNAs was significantly lower than that of mRNAs.

Analysis of Differentially Expressed lncRNAs
Based on the cutoff criteria for distinguishing differentially expressed lncRNAs, 168 lncRNAs were significantly FIGURE 1 | Circos plot overview of lncRNA sequencing data. The inner plots show chromosomes 1 to 32 shared between GCs from SYFs in the RL and WL groups. different between RL and WL groups (p < 0.05), of which 101 were downregulated and 67 were upregulated in the RL group (Figure 3, Supplementary Table 4). Furthermore, 1,182 differentially expressed mRNAs were identified (958 downregulated and 224 upregulated). RT-qPCR was applied to validate the RNA-Seq results using six candidate lncRNAs and eight mRNAs; cytochrome P450 family 11 subfamily A member 1 (CYP11A1), cytochrome P450 17A1 (CYP17A1), 7β-hydroxysteroid dehydrogenase type 7 (HSD17B7), bone morphogenetic protein 15 (BMP15), bone morphogenetic protein receptor 2 (BMPR2), luteinizing hormone/choriogonadotropin receptor (LHCGR), follicle stimulating hormone receptor (FSHR), and insulin-like growth factor 1 receptor (IGF1R). The results showed that the expression patterns of these mRNAs and lncRNAs were similar to the results of high-throughput sequencing (Figures 4, 5), indicating that the RNA-Seq data were reliable.

Interactions Between lncRNAs and mRNAs Involved in Reproduction
Co-expression of different lncRNAs and mRNAs was analyzed using Pearson correlation tests to calculate expression correlations between differentially expressed lncRNAs (length <6000 bp) and differentially expressed mRNAs. Correlation coefficient ≥ 0.9 and p-value ≤ 0.01 were selected as the cutoff criteria, and 120 lncRNAs targeting 1,175 mRNAs were identified (Supplementary Table 5). Potential cis and trans targets of lncRNAs were predicted to explore the functions of lncRNAs. Nine lncRNAs were identified with cis functions and 39 lncRNAs were predicted with trans functions. Four lncRNAs (ENSGALT00000099549, ENSGALT00000103341, ENSGALT00000097648, and XR_001464460.2) were predicted to target 344 mRNAs (Figure 6, Supplementary Table 6). ENSGALT00000099549, ENSGALT00000097648, and XR_001464460.2 target a set of mRNAs related to ovarian steroidogenesis, including CYP11A1, BMP15, HSD17B7, and LHCGR. ENSGALT00000103341 targets a set of mRNAs related to sphingolipid metabolism, including ABC transporters such as UGT8, UGT8L, ABCB5, ABCA4, and ABCG2. XR_003077640.1 regulates genes such as CYP17A1, IGF1R, and BMPR2, while XR_001468386.1 regulates genes including PIK3R1, PAK5, SEMA3C, and BMPR2. Furthermore, the target mRNAs of these lncRNAs were enriched in various GO functions including negative regulation of cell proliferation, cell differentiation, and the Notch signaling pathway (Figure 7). Meanwhile, pathway enrichment analyses identified ovarian steroidogenesis, cell adhesion molecules, and ABC transporters (Figure 8).

DISCUSSION
Our results demonstrate that monochromatic light affects the egg laying performance of Jinghai Yellow chickens. Compared with WL, RL improved egg production, in agreement with a previous study on hens . Foster and Follett (1985) reported that RL elicits a greater photosexual response than the same photon flux of other wavelengths, because longer wavelengths penetrate more easily to reach the hypothalamus. Pyrzak and Siopes (1986) suggested that RL increases egg weight, but our results showed that RL had only a minor effect on the egg weight of chickens. Specifically, the body weight of chickens at 300 days of age under RL was similar to that under WL, consistent with Li et al. (2014). RL can clearly improve the egg laying performance of chickens, but the molecular mechanisms associated with the effects of monochromatic light remain obscure.
We compared the lncRNA expression profiles of GCs from SYFs of Jinghai Yellow hens under RL and WL, and identified a number of lncRNA target mRNAs related to egg production. Compared with the WL group, 168 lncRNAs were differentially expressed in the RL group (101 downregulated and 67 upregulated). Target prediction and functional analysis of these genes, including ENSGALT00000099549, ENSGALT00000103341, ENSGALT00000097648, XR_001468386.1, XR_003077640.1, and XR_001464460.2, showed that many lncRNA target mRNAs were associated with steroid hormone biosynthesis, ovarian steroidogenesis, TGF-β signaling pathway, prolactin signaling pathway, MAPK signaling pathway, estrogen signaling pathway, PI3K-Akt signaling pathway, Hippo signaling pathway, MAPK signaling pathway, and cell adhesion molecules. In our previous studies, steroid hormone biosynthesis and the PI3K-Akt signaling pathway were also altered in pigeons under monochromatic light (Wang et al., 2015. Huang et al. (1979) suggested that progesterone is synthesized by GCs, and Hu et al. (2017) concluded that PI3K and AMPK signaling pathways converge to modulate ERK activity, and thereby regulate GC differentiation. Meanwhile, Kim et al. (2013) showed that the TGF-β signaling pathway is involved in initiating follicle selection, in accordance with our current results, suggesting that these pathways play a key role in follicular development under monochromatic light.
Eight overlapping targets identified by RNA-Seq were selected for verification by RT-qPCR, including CYP11A1, CYP17A1, and 17β-HSD, which are important for steroidogenesis and the synthesis of hormones such as testosterone, progesterone, and estrogen (Zhou et al., 2020). P450scc, the protein product of the CYP11A1 gene, plays a major role in the control of steroidogenesis (Simpson, 1979;Hara et al., 2014). In chicken follicular GCs, expression of CYP11A1 is a prerequisite for progesterone synthesis, and is related to follicle selection. Nebert et al. (2013) reported that CYP11A1 knockout mice display a variety of aberrant phenotypes associated with various steroid hormone deficiency syndromes. CYP11A1 and CYP17A1 have the capacity to synthesize androgenic substrates that diffuse into adjacent pre-granulosa cells (Lagaly et al., 2008). HSD17B7 converts estrone to estradiol, HSD17B7 is involved in cholesterol biosynthesis and has 3-ketosteroid reductase activity (Nokelainen et al., 1998;Shehu et al., 2008), and mouse HSD17B7 −/− fetuses lack proteins involved in cholesterol synthesis (Saher et al., 2005).
Compared with large follicles, GCs of small follicles have lower CYP11A1 levels. Hatzirodos et al. (2014) also pointed out that early in follicle development, the hormone inhibin is secreted, and oestradiol is later secreted at the preovulatory stage, which may explain why these genes were down-regulated in the RL group in the present work.
Follicular development and differentiation of GCs are dependent on regulation of LH and FSH, and their specific receptors (Richards et al., 1976). LH induces progesterone secretion through LHCGR, and an increase in granulosa LHCGR causes a progressive increase in the responsiveness of GCs to LH in maturing follicles (Bahr and Johnson, 1984). LHCGR is differentially regulated between small and large follicles (Peng  , 1991), and is downregulated in small follicles (Hatzirodos et al., 2014), in accordance with our current results. Recent studies have shown that all SYFs isolated from the same ovary in a laying hen can express FSHR, and respond to stimulation by FSH (Webb et al., 2003;Johnson et al., 2015). One of the earliest markers for differentiating GCs is elevated expression of FSHR mRNA, specifically within the granulosa layer (Johnson and Woods, 2009). In hens, SYFs expressing the highest levels of FSHR are recruited into the preovulatory hierarchy during ovarian follicle development (Woods and Johnson, 2005;Wang et al., 2017). Herein, GCs of SYFs in the RL group had higher FSHR levels, indicating more preovulatory hierarchy follicles.
As members of the TGF-β superfamily, BMPs and their receptors have been shown to play important roles during folliculogenesis (Juengel and McNatty, 2005). BMP15 promotes follicle selection in hens (Stephens and Johnson, 2016). Moore et al. (2003) demonstrated that BMP15 promotes mitosis of GCs and suppresses FSHR expression, leading to the suppression of FSH-induced progesterone synthesis, and stimulation of kit ligand expression, in accordance with the results of the present study. FSHR was upregulated and BMP15 was downregulated in the RL group. BMPR2 binds GDF-9 and BMP-15 (Hatzirodos et al., 2014), and dysregulation of BMPR2 gene expression has been associated with abnormalities in follicular development (Andreas et al., 2016). In vivo, GCs require IGF1R to undergo differentiation upon FSH stimulation (Zhou et al., 2013), and upregulation of IGF1R in the RL group is in accordance with FSH. Baumgarten et al. (2017) concluded that IGFR is necessary for the formation of preovulatory follicles and the differentiation of GCs to the preovulatory stage.
In conclusion, we found that RL enhanced the egg production of hens. High-throughput sequencing identified a set of differentially expressed lncRNAs and mRNAs in GCs of SYFs between RL and WL groups. The results provide a starting point for studies aimed at understanding the molecular mechanisms by which red light increases the egg production of Jinghai Yellow chickens. Furthermore, the findings will facilitate research on the effects of red light on follicular development in Jinghai Yellow chickens and other poultry.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI Sequence Read Archive (SRP278036), BioProject (PRJNA657681).

ETHICS STATEMENT
The animal study was reviewed and approved by Institutional Animal Care and Use Committee of the Department of Animal Science and Technology, Yangzhou University, China.

AUTHOR CONTRIBUTIONS
YW, JW, and GZ contributed to the overall design of the study. YW, HS, and YG collected data. YW, HS, MS, TL, XL, and LC contributed to sample collection. YW participated in manuscript writing and revision. All authors approved the final version of the manuscript.