All experimental procedures involving the use of animals were conducted according to the Guidelines for Experiment established by the Science Research Department of the Institute of Animal Sciences, Chinese Academy of Agriculture Sciences (IAS-CAAS), and ethical approval was received from the Animal Ethics Committee of the IAS-CAAS (no. IAS 2022-35).

In the present study, White Leghorn and Beijing You chicken were employed as parents to generate purebreds (WW and YY) and reciprocal crossbreeds (WY and YW). Briefly, 30 males and 150 unrelated females from White Leghorn, and 30 males and 150 unrelated females from Beijing You chicken were selected for mating to generate WW and YY and for reciprocal mating to generate WY and YW based on laying consistency, semen quality and body weight. Chickens from a single hatch were raised under the same environment with free access to feed and water in the experimental farm of IAS-CAAS.

Egg was collected for each individual daily from the age at first egg (AFE). When the egg laying ratio of population reached 5%, thirty chickens from each group were randomly selected to measure pubic space. Meanwhile, six individuals from WW, YY, WY and YW were selected based on the average of pubic space of each group, respectively. Subsequently, all visible follicles were carefully excluded when ovary tissue was collected as previously reported (23). Then, the ovary was frozen immediately in liquid nitrogen for RNA sequencing. In addition, the oviduct was captured for length measurement by image J (24). Heterosis of phenotype was calculated according to the following equation:

Where F 1¯ , PM¯ and PF¯ are the mean phenotype of reciprocal crossbreeds, maternal and purebreds, respectively. In addition, Student’s t-value was estimated to evaluate the significance of H% based on the formula of Wu and Zhang (25):

where F_1i is the phenotype of individual i in crossbreeds; N is the number of chickens in WY or YW. P value was obtained according to the student’s t-test and degrees of freedom. H% was considered statistically significant if P < 0.05.

Total RNA was extracted using TRIzol^® Reagent (Invitrogen, USA) following the manufacturer’s instructions. The RNA integrity and quality were determined by agarose gel electrophoresis and NanoPhotometer^® spectrophotometer (IMPL EN, CA, USA). Three micrograms of total RNA were used for the construction of lncRNA library, and then the ribosomal RNA was removed from total RNA with TruSeq Stranded Total RNA Library Prep Gold Kits (Illumina, United States). Finally, 24 RNA sequencing libraries were constructed for paired-end sequencing according to the instructions of NEB Next^® Ultra Directional RNA Library Prep Kit for Illumina^® (NEB, Ipswich, MA, USA). RNA-Seq were performed using Novaseq6000 (Illumina, San Diego, United States) to generate 150bp paired-end reads.

The reads, including adapter contamination, low-quality reads (Phred Quality Score < 5%), reads with poly-N > 5% and reads matched rRNA were filtered out using in-house perl scripts to generate clean reads, which were aligned to the chicken reference genome (http://ftp.ensembl.org/pub/release-106/fasta/gallus_gallus/) using HiSAT2 (v2.1.0) (26) with default parameters, and then the mapped reads was assembled by StringTie (v2.1.5) (26) with gene transfer format (GTF) file of Ensembl gene annotation (http://ftp.ensembl.org/pub/release-106/gtf/gallus_gallus/). According to the chicken gene annotation, the number and proportion of three lncRNA functional elements (exon, intron and intergenic) were calculated based on the unique mapping sequences.

The potential lncRNAs were winnowed with the following criteria. Firstly, the transcripts less than 200 bp and without strand information were removed. Secondly, transcripts with class code “i”, “u” and “x” were retained. Subsequently, the sequence of known lncRNAs were download from two lncRNA databases, Domestic-Animal LncRNA Database (ALDB) (12) and NONCODE (27). We then used BLAST to align the unannotated transcripts to known lncRNAs. The transcripts perfectly aligned with sequences in either ALDB or NONCODE were regarded as known lncRNAs. Then, the protein coding potential of the remaining transcripts was predicted using CNCI (28), CPC (29), PLEK (30) and Pfam (31) software. Only transcripts with CNCI score < 0, CPC score < 0, PLEK score < 0 and Pfam E-value < 1e-5 were established as putative lncRNAs. Finally, the cis- and trans- acting relationship between lncRNAs and potential protein-coding genes was predicted according to their distances and expression correlations (32).

Differentially expressed lncRNAs (DELs) and differentially expressed genes (DEGs) between two different groups were analyzed based on gene counts by DESeq2 (v1.16.1) (33) in R project. Meanwhile, |Log2(fold change) | >1 and adjusted P value < 0.05 were taken as criteria to identify the DELs and DEGs in the corresponding comparison. These genes were further classified into three inheritance patterns: additivity, dominance and overdominance, based on the gene expression level of purebreds and crossbreeds as previously reported (34). Briefly, additivity (IV and X) occurred when gene expression was significantly different between two purebreds, and gene expression of crossbreeds was higher than one purebred but lower than the another. Gene expression within crossbreeds that was not significantly different from one purebred but significantly higher (or lower) than the another was regarded as dominance (II, IV, IX, and XI). Gene expression within crossbreeds that was significantly higher (or lower) than both purebreds was considered as overdominance (V, VI, VIII, III, VII, and X).

The weighted gene co-expression network analysis (WGCNA) was performed to identify co-expressed modules of highly correlated genes and summarized such modules based on the correlation between module eigengene and phenotypes (35). Briefly, an expression matrix of all lncRNAs and mRNAs was constructed. The top 40% most variant genes were used for subsequent analysis. After the weighted adjacency matrix established, an unsigned weighted correlation network was constructed by creating a matrix of Pearson correlation coefficients of each pair genes. Then, each topological matrix was used as input for linkage hierarchical clustering analysis, and primary gene modules was detected using a dynamic tree cutting algorithm (deepSplit = 2, minModuleSize = 50, mergeCutHeight = 0.25). Subsequently, module eigengenes were correlated with different traits and searched against the most significant correlations. Finally, the lncRNAs-mRNAs networks were visualized with Cytoscape (v.3.8.0), which was run with layout “attribute circle layout”.

To investigate biological function of nonadditive genes, we performed functional enrichment analysis, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways via KOBAS 3.0 platform (36). The GO terms and KEGG pathways with P value < 0.05 were considered significantly enriched.

GnRH and MT content was detected using the GnRH and MT assay kit (Beijing Laibotairui Technology Co. Ltd, Beijing, China) following the manufacture guidelines, respectively. Briefly, serum (40 μL) from each chicken was added to each well of 96-well plate. Then, the 96-well plate was incubated for one hour at 37°C. Following washed five times with washing liquid, color liquid A and B was added successively. After incubating for 10 min at 37°C, 50 μL of termination solution was added. Finally, the OD value was detected at 450 nm by Thermo Multiskan Ascent.

One microgram total RNA of each sample was reverse transcribed into cDNA using the PrimeScript RT Reagent Kit (TaKaRa, Shiga, Japan). The specific primer sequences were designed using Primer Premier5.0. The RT-qPCR reaction mixture (10 µL) consisted of cDNA 1.5 µL, forward primer 0.5 µL (10µM), reverse primer 0.5 µL (10µM), SYBR Green Master (Mix) 5 µL and ddH₂O 2.3 µL. The amplification condition was as follows: 95°C for 3 min, followed by 40 cycles at 95°C for 3 s, 60 °C for 32s, finally, a single melt cycle was 95°C for 15 s and 65°C for 1 min. Triplicate was performed for each sample using PrimeScript One Step RT-PCR Kit (TaKaRa, Shiga, Japan). The results were normalized to GAPDH expression.

Fold changes in gene expression were calculated using the 2^−ΔΔCt method. All data are presented as the mean ± standard deviation (SD). The correlation between results of RT-qPCR and RNA-seq was calculated using Pearson^’s correlation method in R (v.4.0.2).

As shown in Figure 1A , the pubic space of crossbreeds was larger than the average value of purebreds, and significant positive heterosis (25.79%, 32.45%) was observed in both WY and YW (P < 0.05). Similarly, the larger oviduct length was observed in WY and YW compared with the average length of purebreds ( Figure 1B ), and heterosis of oviduct length was 30.55% in WY and 20.15% in YW, respectively. Moreover, AFE in crossbreeds was earlier than the average age of purebreds exhibiting negative heterosis ( Figure 1C ).

In the present study, an average of 36,579,837 raw reads was generated from each library, with average Q30 bases higher than 94.56% ( Supplementary Table S1 ). The percentage of uniquely mapped reads and multi-mapped reads was 82.13% and 2.73%, respectively. In total, we characterized 3164 lncRNAs in chicken ovary (details in Supplementary Table S1 ). Among these lncRNAs, 1170 known lncRNA were identified through aligning lncRNAs sequence to the NONCODE or ALDB database. Besides known lncRNAs, a total of 1994 transcripts were identified as putative lncRNAs by CNCI, CPC, PLEK and Pfam software ( Figure 2A ). The majority of putative lncRNA were intergenic (50.82%), and 25.76% were intron, and 23.42% were antisense lncRNAs ( Figure 2B ). Furthermore, the length of the putative lncRNAs was mainly ranging from 2,800 to 3,000bp, similar to that of the known lncRNAs in ovary tissue ( Figure 2C ). The distributions of putative and known lncRNAs on the chromosomes were varied ( Supplementary Figure S1 ). Moreover, the number of lncRNAs located on chromosome 1 was the largest, which accounted for 8.93% in putative lncRNAs and 13.08% in known lncRNAs. The most putative and known lncRNAs possessed two exons. The exon numbers of putative lncRNAs were greater than known lncRNAs ( Figure 2D ).

A principal component analysis (PCA) was performed to visualize the group differences of gene expression. The purebreds and crossbreeds were obviously separated from each other ( Supplementary Figure S2 ), indicating significant differences of gene expression between two purebreds and between purebreds and crossbreeds. With the criteria of |log2 (fold change) | > 1 and adjusted P value < 0.05, we identified 686 DELs between two purebreds and between purebreds and crossbreeds ( Figure 3A ), including 630 in YY vs. WW, 225 in WY vs. WW, 83 in WY vs. YY, 246 in YW vs. WW and 68 YW vs. YY ( Supplementary Table S2 ). These DELs were clustered into 12 types (I, II, III, IV, V, VI, VII, VIII, IX, X, XI and XII, details in Supplementary Table S3 ) according to their expression in purebreds and crossbreeds. The 12 types were further clustered into three main inheritance patterns: additivity (IV, X), dominance (III, V, IX, XI) and overdominance (I, II, VI, VII, VIII, XII). The number of dominant lncRNAs was 288 in WY and 294 in YW, respectively. The number of overdominant lncRNAs was 2 and 2 in WY and YW, respectively. The nonadditive lncRNAs accounted for 44.24% in WY and 45.16% in YW, respectively ( Figures 3B, C ). Additionally, a total of 200 nonadditive lncRNAs were shared by WY and YW ( Supplementary Figure S3A ). The number of unique nonadditive lncRNAs was 90 in WY and 96 in YW, respectively.

There were 2023 DEGs between two purebreds and between crossbreeds and purebreds identified ( Figure 3D ), including 1918 in YY vs. WW, 899 in WY vs. WW, 131 in WY vs. YY, 781 in YW vs. WW and 220 YW vs. YY ( Supplementary Table S2 ). These DEGs were clustered into 12 types, and the gene number of different patterns in reciprocal crossbreeds was shown in Figure 3E . The dominant genes were 956 in WY and 868 in YW, respectively. The number of overdominant genes was 3 and 1 in WY and YW, respectively. The nonadditive genes, accounted for 49.30% in WY and 45.42% in YW ( Figure 3F ). Among the nonadditive genes, a total of 682 genes were shared by WY and YW, and the number of unique genes was 274 in WY and 186 in YW ( Supplementary Figure S3B ), respectively. As shown in Figure 4A , gene ontology (GO) analysis showed that the nonadditive genes were mainly enriched in three categories, including growth and development related GO terms, such as animal organ development, gland morphogenesis, female gonad development and female sex differentiation, ion transport related GO terms, such as calcium ion transport, calcium ion homeostasis and metal ion transport, and other GO terms, such as response to hormone, steroid binding, response to endogenous stimulus and sterol transport (details in Supplementary Table S4 ). Interestingly, the majority of GO terms were also enriched by shared nonadditive genes. Given that the nonadditive genes of WY and YW were also significantly enriched in the biological process of female gonad development, we further characterized the common nonadditive genes harbored in the GO terms shared by WY and YW. The four common nonadditive genes, including transforming growth factor beta 1 (TGFB1), collagen type VI alpha 1 chain (COL6A1), CCAAT/enhancer binding protein beta (CEBPB) and receptor tyrosine kinase (KIT) genes, were involved in the growth and development of ovarian cells ( Supplementary Table S4 ).

KEGG pathway analysis showed that nonadditive genes were significantly enriched in the focal adhesion, ECM-receptor interaction, vascular smooth muscle contraction, GnRH signaling pathway, calcium signaling pathway and folate biosynthesis ( Figure 4B ). These pathways were also enriched by shared nonadditive genes. Regarding the GnRH signaling pathway, we found eight common nonadditive genes including G protein subunit alpha q (GNAQ), calcium voltage-gated channel subunit alpha1 C (CACNA1C), ENSGALG00000005884, adenylate cyclase 6 (ADCY6), calcium voltage-gated channel subunit alpha1 D (CACNA1D), ENSGALG00000008727, early growth response 1 (EGR1) and matrix metallopeptidase 2 (MMP2) were enriched in the pathway ( Supplementary Table S4 ).

We constructed a co-expression network of lncRNAs and mRNAs to infer the underlying regulatory function and potential co-expressed genes and lncRNAs. Correlated lncRNAs and mRNAs were clustered into 20 modules by WGCNA, and each module contained at least 50 genes ( Figure 5A ). Based on the correlation between the clusters and phenotypes, seven modules were significantly correlated with pubic space or oviduct length (P < 0.05) ( Figure 5B ). Among these modules (ME), the midnightblue ME (r = -0.57, p = 0.01) and greenyellow ME (r = -0.56, p = 0.01) were negatively correlated with pubic space. The purple ME was negatively correlated with pubic space (r = -0.61, p = 0.005) and oviduct length (r = -0.57, p = 0.01). The green ME was positively correlated with the oviduct length (r = 0.53, p = 0.02). The turquoise ME (r = -0.59, p = 0.008), the brown ME (r = -0.54, p = 0.02) and the red ME (r = -0.64, p = 0.003) was negatively correlated with oviduct length. Additionally, for all lncRNAs and protein coding genes harbored in the co-expression network, the top 1% of highly correlated gene pairs were chosen to visualize the correlation of modules ( Figure 5C ). Function enrichment results of each cluster suggested that the highly correlated lncRNAs and genes were mainly enriched in the developmental related process. Notably, the annotated genes in cluster green ME, were significantly enriched in oocyte meiosis, ECM-receptor interaction, calcium ion binding and GnRH signaling pathway, suggesting that lncRNAs in this cluster is essential for those biological processes. In addition, genes in turquoise ME were found closely correlated with fatty acid biosynthetic process, cell cycle and peptide biosynthetic process ( Figure 5D ).

We found 43 nonadditive genes were interactive and harbored in the seven associated modules, and these genes were regulated by 40 nonadditive lncRNAs in cis or trans ( Figure 6A ). Interestingly, function of these nonadditive genes were associated with reproductive structure development, gonad development, female gonad development, calcium signaling pathway and GnRH signaling pathway. The number of common nonadditive genes enriched in female gonad development and GnRH signaling pathway was four and eight in crossbreeds ( Figure 6B ), respectively. Moreover, CNCAN1C and TGFB1, which were regulated by MSTRG.17017.1 and MSTRG.6475.20 in trans, were enriched in GnRH signaling pathway and female gonad development, respectively. Accordingly, we found the expression of FSHR and ESR2 in crossbreeds was both higher than that of purebreds. But these two genes were not nonadditive genes ( Supplementary Figure S4 ).

The nonadditive genes were related to GnRH signaling pathway and female gonad development, implying that hormone secretion plays a vital role in positive heterosis of sexual maturation. Then, the serum GnRH and MT concentration were detected. GnRH content in crossbreeds was significantly higher than the average value of purebreds, and showed positive heterosis (P < 0.05) ( Figure 7A ). MT content was lower in reciprocal crossbreeds compared with the average value of purebreds exhibiting negative heterosis ( Figure 7B ).

To validate the RNA sequencing results, the expression of 13 DEGs and 6 DELs was detected using RT-qPCR. The specific primers for genes and lncRNAs were listed in Table 1 . The expression of DEGs and DELs between any two groups showed consistent trend with the results of RNA-seq. A high correlation coefficient of gene expression level was detected (R² = 0.91), indicating that RNA-seq data were reliable ( Figures 7C, D ).