To create F₁ hybrids, gametophytes were cultivated from the parental diploid species, Phegopteris koreana and P. taiwaniana, on agar media in petri dishes. The spores sown were collected from two samples of P. koreana and two samples of P. taiwaniana ( Supplementary Table S1 ). The culture conditions were primarily based on those outlined by Watano and Masuyama (1991). From each of the two species, one gametophyte (total of two) was transplanted into the same well of a 12-well petri dish. Sporophytic outcrossing or gametophytic selfing (Haufler et al., 2016) was then promoted with the addition of two to three drops of sterilized water to each well once a week. The resulting sporophytes were transferred into pots individually and cultivated in the experimental field at the Nishi-chiba campus of Chiba University. To evaluate the hybridity of the produced sporophytes, Sanger sequencing was conducted for the plastid DNA marker, rbcL, and the single-copy nuclear gene, PgiC. Subsequently, phylogenetic analyses were performed for each gene by incorporating the sequences into the previous sequence matrix used by Fujiwara et al. (2021). The experimental procedures, from PCR to phylogenetic analysis, followed the methods described by Fujiwara et al. (2021). As a result, three F₁ hybrid individuals with P. koreana as a maternal parent (referred as to F_1K) and two F₁ hybrid individuals with P. taiwaniana as a maternal parent (referred as to F_1T) were obtained. The maternal diploid ancestor of the allotetraploid P. decursivepinnata has been proven to be P. koreana (Fujiwara et al., 2021), and thus the F_1K hybrids were verified to be the products of crossing in the same direction as P. decursivepinnata, while the F_1T hybrids were in the opposite direction. The sporophytes genetically identified as F₁ hybrids were used for subsequent analyses ( Supplementary T able S2 ).

Sporophytes of the parental diploid species, Phegopteris koreana and P. taiwaniana, were also produced by artificial fertilization of the gametophytes. This was necessary as the parental sporophytes of the F₁ hybrids ( Supplementary Table S1 ) were not cultured. Two gametophytes of the same species were transplanted into the same well of a 12-well petri dish, and sporophytic selfing or gametophytic selfing was promoted as previously described. Genotypes of the resulting sporophytes were determined using the same methods as described for F₁ hybrid identification. Ultimately, three P. koreana and three P. taiwaniana sporophytes were selected for subsequent analyses ( Supplementary Table S2 ).

For the allotetraploid species, Phegopteris decursivepinnata, three sporophyte samples were collected ( Supplementary Table S1 ), and then cultivated in an experimental field alongside those of the F₁ hybrids and parental species, for the period from April 2021 to April 2022 prior to RNA extraction.

Approximately 2 × 2 cm squares of the unfolding young leaves were collected in April 2021 from three Phegopteris koreana, three P. taiwaniana, and three F_1K individuals, and in April to May 2022 from two F_1T and three P. decursivepinnata individuals. Additionally, individuals representing each parental species (K1 and T1) were resampled in April 2022 to assess the reproducibility of results across the two years ( Supplementary Ta ble S2 ). The collected leaf samples were rapidly frozen in liquid nitrogen and subsequently crushed using zirconia beads. Total RNA was extracted using the PureLink™ Plant RNA Reagent (ThermoFisher Scientific). The RNA-seq libraries were constructed using NEBNext^® Poly(A) mRNA Magnetic Isolation Module (for PolyA selection) and NEBNext^® Ultra™ II Directional RNA Library Prep Kit according to the manufacturer’s instructions. Illumina 150 bp paired-end sequencing was performed in NovaSeq6000. The raw read data was deposited in the DDBJ Sequence Read Archive under the accession numbers (DRR498230–DRR498244).

For all RNA-seq reads, adapter sequences and low-quality bases (Q < 20) were trimmed using Trimmomatic v.0.39 (Bolger et al., 2014). After trimming, the filtered reads from biological replicates of each species ( Supplementary T able S2 ) were assembled using Trinity v.2.9.1 (Grabherr et al., 2011). To remove isoform redundancy in the Trinity assembly outputs, the highest expressed isoforms within genes estimated by RSEM (Li and Dewey, 2011) were selected using the “filter_low_expr_transcripts.pl” provided in Trinity v.2.9.1. TransDecoder v.5.5.0 was used to recognize the ORFs, and only the longest ORF per transcript was retained. The obtained sequences were annotated with TAIR IDs based on the best BLASTP hit with a cutoff E-value of 10^–5 against the Arabidopsis protein sequence database (Araport11_pep_20220914_representative_gene_model). Organelle-encoded genes were removed from subsequent analyses. After the above filtering steps, orthogroups were constructed among parental species, F_1K, F_1T, and Phegopteris decursivepinnata using OrthoFinder v.2.5.2 (Emms and Kelly, 2019). The venn diagrams showing the numbers of coexpressed orthologs within orthogroups were drawn using the venn package v1.11 in R software (R Core Team, 2022). The orthogroups containing a single gene each in parental species were used as reference gene sets ( Supplementary Table S3 ).

For the differentially expressed gene (DEG) analyses of the homoeologs in the F₁ hybrids and the allotetraploids, the reads of the parental species were mapped to their own reference gene sets and the reads from the F₁ hybrids and the allotetraploids were then mapped to the reference genes in both P. koreana and P. taiwaniana using bowtie2 v.2.4.5 (Langmead and Salzberg, 2012). Considering that the F₁ hybrids and the allotetraploids retain homoeologs from both parental species, the output bam files from bowtie2 were used as input for the file for EAGLE-RC (Kuo T. et al., 2018), a tool that identifies the parental origin of the reads from the F₁ hybrids and the allotetraploids. The number of mapped reads was counted using eXpress v.1.5.1 (Roberts and Pachter, 2013). The edgeR package v.3.40.2 (Robinson et al., 2010) in R, was used to determine the read count data for the replicates and this was normalized using the trimmed mean of M values (TMM) methods, and differential expression analyses were performed between parental species, between the F₁ hybrids and each parent, and between the allotetraploids and each parent. In the present study, the DEGs were identified using Fisher’s exact test, with the following criteria: |log₂(fold change (FC))| >1 and false discovery rate (FDR) < 0.05. To evaluate the consistency of the trends of the expression patterns of the entire genomes/subgenomes, we conducted the three patterns of transcripts per million (TPM) filtering, TPM > 0.0, TPM > 0.5, and TPM > 1.0, in all five groups.

Genes identified as DEGs in the F₁ hybrids and the allotetraploids relative to their diploid parents were classified into 12 expression-level categories ( Figure 1 ), using the differential expression classification outlined by Rapp et al. (2009). These categories include the following: additivity (I and XII), Phegopteris taiwaniana ELD (T-dominance; II and XI), P. koreana ELD (K-dominance; IV and IX), transgressive expression lower than either parent (Transgressive Down-regulation; III, VII, and X), and transgressive expression higher than either parent (Transgressive Up-regulation; V, VI, and VIII). Genes that did not show significant differential expression levels were categorized as “No Change.” When comparing the F₁ hybrids to the allotetraploids, the 12 categories were grouped into five broader categories (K-dominance, T-dominance, Additivity, Transgressive Up-regulation, and Transgressive Down-regulation). The changes among the five differential expression categories from the F₁ hybrids to the allotetraploids were visualized using the qgraph package (v.1.9.5) in R.

Using the edgeR package in R, the expression level ratio, log₂(FC), of the orthologous genes from Phegopteris taiwaniana and P. koreana was calculated. This calculation involved comparing the mean count-per-million (CPM) of the alleles from P. taiwaniana (Pt) replicates with the mean CPM of those from the P. koreana replicates (Pk), as follows: log₂(Pt_{_ortholog}/Pk_{_ortholog}). The same method was used to determine the expression level ratios of the parental alleles in the F₁ hybrids, as follows: log₂(F_{1_T-homoeolog}/F_{1_K-homoeolog}), and those of the homoeologs in the allotetraploid P. decursivepinnata (Pd): log₂(Pd_{_T-homoeolog}/Pd_{_K-homoeolog}). Subsequently, the distribution of the log₂(FC) was compared among the parental species, the F₁ hybrids, and the allotetraploids.

The expression bias of the homoeologs in the F₁ hybrids and the allotetraploids was analyzed using the existing difference in gene expression levels between the parental species. To achieve this, the ratio of the mean CPM of the orthologous genes in the parental species and the ratio of the mean CPM of homoeologous genes in each of the F₁ hybrids or the allotetraploids was compared using Fisher’s exact test. The ratio of the log₂ (FC) for the orthologous genes in the parental species and that of the homoeologs in the F₁ hybrids was then determined as follows: log₂((F_{1_T-homoeolog}/F_{1_K-homoeolog})/(Pt_{_ortholog}/Pk_{_ortholog})). Similar comparisons were conducted between the parental species and the allotetraploids, as follows: log₂((Pd_{_T-homoeolog}/Pd_{_K-homoeolog})/(Pt_{_ortholog}/Pk_{_ortholog})). The ratio of log₂(FC) was referred to as the “magnitude” of the bias. The HEB was detected for each of the homoeologous gene pairs using the following criteria: P-value <0.05 and FDR <0.05 with the Fisher’s exact test, and the absolute value of the “magnitude” >1. The homoeologous gene pairs with a “magnitude” >1 and <−1 are considered as HEB toward P. taiwaniana (T-bias) and P. koreana (K-bias), respectively. The relationship of the gene sets exhibiting the HEB between the F₁ hybrids and the allotetraploids are shown using Circular layout graphs using the qgraph packages in R. The plots that show the relationship between the log₂(FC) values of parental orthologs and those of homoeologs and the histograms showing the difference in the expression bias between the F₁ hybrids or the allotetraploids and the parents were generated using the ggplot package in R.

The artificial crossing experiment successfully produced several sporophyte individuals. The nuclear PgiC tree confirmed that five of these were interspecific F₁ hybrids between Phegopteris koreana and P. taiwaniana ( Supplementary Figure S1 ). The plastid rbcL tree showed that there were three and two hybrids with P. koreana and P. taiwaniana as their mother, respectively, given that maternal inheritance of organelle genomes in ferns (reviewed in Kuo L.-Y. et al., 2018) ( Supplementary Figure S1 ; Supplementary Table S2 ).

Approximately, 45–72 millions of 150 bp pair-end raw reads were sequenced in each studied sample. After filtering the RNA-seq reads for quality, the initial transcripts were assembled and filtered to 26,663–31,796 nuclear genes containing complete or partial ORFs across five groups containing three species and reciprocal F₁ hybrids of the parental species ( Supplementary Table S3 ). The ortholog identification using these filtered genes yielded 31,936 orthogroups among the five groups. Of these, 12,380 genes that were detected in both parental species and in which no gene duplication was identified for each parental species, were retained as mapping references for subsequent gene expression analysis ( Supplementary F igure S2 ). To reveal the effects of the low-expressed genes, we performed three levels of TPM filtering for all five groups, and 12,380 orthologous genes were retained with TPM > 0, 9,540 with TPM > 0.5, and 9,003 with TPM > 1.0 ( Supplementary Table S3 ). MDS plot of the expression matrix ( Supplementary Figure S3 ) showed that the replicates of each parental species, F_1K hybrids, F_1T hybrids, and allotetraploids, Phegoteris decursivepinnnata are closely clustered together, suggesting that the expression pattern is highly consistent among the replicate samples of each group. The samples from each parental species include those collected in different years (2021 and 2022), but all are closely clustered, and there is no apparent influence due to the difference in years. Furthermore, the clustering result is congruent among the datasets with different TPM filtering. In the comparison of expression levels across species, genes showing stable expression should be used because genes that respond to transient environmental changes should be excluded from the analyses. Therefore, we mainly showed the results obtained from the dataset with TPM > 0.5 for the subsequent analysis.

To identify additivity, transgressive expression, and ELD in the F₁ hybrids and the allotetraploids, in comparison to the expression levels in the parental species, the 9,540 orthologous genes were classified into 12 categories, following the categorization system described by Rapp et al. (2009) ( Figure 1 ). The analysis revealed that most genes from the F₁ hybrids and the allotetraploids did not display significant differential expression levels when compared to the orthologs of their parents (denoted as “No Change” in Figure 1 ; 80.1% in the F_1K hybrids, 81.8% in the F_1T hybrids, and 77.4% in the allotetraploids), with a lesser proportion being observed in the allotetraploids (2.7% and 4.4% decreases compared to the F₁ hybrids). Few differences were found in the count for “No Change” genes between the hybrids with Phegopteris koreana (F_1K) and P. taiwaniana (F_1T) as their mothers, respectively. For “Additivity,” the category I and XII showed one of the lowest proportions among the all of categories in the F_1T hybrids, the F_1K hybrids, and the allotetraploids, and the proportion of the categorized genes was almost consistent between the F₁ hybrids and the allotetraploids.

For transgressive expression, a substantial number of genes were found to display lower and higher expression levels when compared to both parents (Transgressive Down- and Up-regulation) in the F₁ hybrids and the allotetraploids. The number of genes in “Transgressive Down-regulation” always outnumbers the number in “Transgressive Up-regulation” across both the F₁ hybrids and the allotetraploids. Interestingly, the number of genes exhibiting transgressive expression was increased in the allotetraploids, and also the difference in the number of genes between “Transgressive Down-regulation” and “Transgressive Up-regulation” became more pronounced in the allotetraploids: from 574 in the F_1K hybrids and 587 in the F_1T hybrids to 807 genes in the allotetraploids across “Transgressive Down-regulation” categories (III, VII and X), and from 352 in the F_1K hybrids and 198 in the F_1T hybrids to 381 genes in the allotetraploids across “Transgressive Up-regulation” categories (V, VIII, and VI).

The analyses have also revealed significant trends in ELD. There were significantly more instances where the F₁ hybrids and the allotetraploids exhibited the same expression level as the parent with higher expression, regardless of either “T-dominance” or “K-dominance” categories. For instance, for the category “T-dominance” in the F_1K hybrids, there were 453 instances where Phegoteirs taiwaniana showed higher expression than P. koreana (category II), compared to only 79 instances where the opposite was observed (category XI). A similar trend was observed in the F_1T hybrids and the allotetraploids ( Figure 1 ). In addition, while a higher number of the genes exhibited “T-dominance” when compared to “K-dominance” in both the F_1K and F_1T hybrids, the allotetraploids showed a similar proportion of the genes exhibiting “T-dominance” and “K-dominance”. In the allotetraploids, 395 and 87 genes in categories II and XI, respectively, exhibited “T-dominance”, and 348 and 108 genes in categories IV and IX, respectively, exhibited “K-dominance”.

The changes in the expression level categories for each gene from the F_1K hybrids to the allotetraploids were then examined ( Figure 2 ). The F_1K hybrids are the product of crossing in the same direction as the natural allotetraploids. Most genes that were categorized as “No Change” (89.4%), “T-dominance” (63.3%), and “K-dominance” (67.3%) in the F₁ hybrids also belonged to the same categories in the allotetraploids. The transition between categories occurred the most frequently between “No Change” and “Transgressive-Up-regulation” or “Transgressive Down-regulation” (208–403 genes). Especially, the transition from “No Change” to “Transgressive Down-regulation” detected almost twice as many cases (403 genes) compared to other transitions between “No Change” and transgressive expressions (208–215 genes). In the genes showing ELD, the number of genes that transitioned between “No Change” and “T-dominance” are almost the same in both directions (71 and 79 genes), while the number of genes that transitioned from “No-Change” to “K-dominance” (92 genes) was almost double that of the other direction (47 genes).

To test the subgenome dominance of Phegopteris koreana or P. taiwaniana in the F₁ hybrids and the allotetraploids, the relative homoeolog contributions and distributions of the expression level ratios for all 9,540 orthologous gene pairs were examined ( Figure 3 ). Although significant deviations from the mean of 0 were detected in the parental pairs, the F_1K hybrids, and the allotetraploids, an apparent expression bias toward either subgenome was not observed in all cases, as the mean and median values were close to 0 ( Table 1 ; Figure 3 ). However, the kurtosis of the distribution was highest in the parental species, and lowest in the allotetraploids. This suggests that the genes with biased expression increased in number after the acquisition of hybridization and polyploidization.

To understand the influence of the expression bias of the homoeologous gene pairs on gene expression, while considering the pre-existing differences in the parental species, the expression level ratios between the homoeologs in the F₁ hybrids and the allotetraploids were compared with those between the parental species. Among the 9,540 orthologous gene pairs, in both the F₁ hybrids and the allotetraploids, approximately 1,000–1,600 pairs exhibited different expression level ratios when compared with the parental ortholog pairs, a phenomenon referred to as HEB ( Table 2 ). Notably, there is no substantial difference between the proportions of T- or K-biased gene pairs observed in both the F₁ hybrids and the allotetraploids (1,400 “T-bias” vs 1,339 “K-bias” in the F_1K hybrids, and 1,020 vs 991 in the F_1T hybrids, 1,619 vs 1,558 in the allotetraploids). Of the gene pairs showing HEB, 60.5% (1658 gene pairs) were shared between the F_1K hybrids and the allotetraploids ( Figure 4 ). The majority of the genes in the other 40% varied between the F_1K hybrids and the allotetraploids by transitioning to or from the “Non-bias” category, while a limited number of gene pairs, only 40 or 41 gene pairs, showed a direct switch between the “T-bias” and “K-bias” categories ( Figure 4 ). More gene pairs altered their expression ratios from “Non-bias” in the F_1K hybrids to “T-bias” or “K-bias” in the allotetraploids than from “T-bias” or “K-bias” to “Non-bias.”

To examine the differences in strength of the HEB between the F_1K hybrids and the allotetraploids, the relative expression levels of the homoeolog pairs in the F₁ hybrids and the allotetraploids were compared with the orthologous gene pairs between the parental species ( Figure 5 ). The F₁ hybrids and the allotetraploids showed a similar trend whereas both the T-biased and K-biased gene pairs showed a similar distribution in the scatter plot and histogram ( Figure 5 ). In the allotetraploids, along with increases in both subgenome-biased gene pairs, the strength of both biases also increased ( Table 3 ; Figure 5 ). The number of genes with a “magnitude” > |10| for T- and K-biased gene pairs increased from nine in the F_1K hybrids and none in the F_1T hybrids to 11 in the allotetraploids for T-biased genes, and five in the F_1K hybrids and two in the F_1T hybrids to 13 in the allotetraploids for K-biased genes, respectively ( Table 3 ). Additionally, the mean values of “magnitude” for T- and K-biased genes were also amplified from 2.72 in the F_1K hybrids and 2.39 in the F_1T hybrids to 2.83 in the allotetraploids for T-biased genes, and −2.36 in the F_1K hybrids and −2.19 in the F_1T hybrids to −2.64 in the allotetraploids for K-biased genes ( Table 3 ). Finally, we explored the functional differences associated with HEB, showing both T-biased and K-biased expression in the F_1K hybrids and the allotetraploids. Among the five K-biased genes and nine T-biased genes in the the F_1K hybrids, and 13 K-biased genes and 11 T-biased genes in the allotetraploids (a “magnitude” >|10|), four, nine, 11 and nine genes were successfully identified using blast analysis against the Arabidopsis protein sequences for K- and T-bias of the F_1K hybrids and the allotetraploids, respectively ( Supplementary Tables S4–S7 ).

Remarkable changes were observed in the gene expression levels of the F₁ hybrids and the allotetraploids when compared to the parental species ( Figure 1 ). Notably, more than half of genes (71.14%) showing non-additive expression in the F_1K hybrids were also detected as non-additive in allotetraploids ( Figure 2 ). These findings suggest that the merger of divergent genomes through hybridization itself had a substantial impact on gene expression changes, which is consistent with previous studies (Hegarty et al., 2006; Chaudhary et al., 2009; Flagel and Wendel, 2010; Buggs et al., 2011; Duan et al., 2023). However, we also observed substantial differences in expression levels between the F₁ hybrids and the allotetraploids. In the allotetraploids, the proportion of genes categorized as “No Change” was decreased compared to those in the F₁ hybrids (77.4% in the allotetraploid vs. 80.1% and 81.8% in the F₁ hybrids).

The gene showing transgressive expression was more frequent in the allotetraploids than in the F₁ hybrids (1,188 genes in the allopolyploids vs. 926 and 785 genes in the F₁ hybrids) ( Figure 1 ). Interestingly, although very low rates of transgressive expression patterns in allopolyploid fern species (approximately 1% or less in both “Transgressive Up-regulation” and “Transgressive Down-regulation”) have been reported in Polypodium hesperium (Sigel et al., 2019), the proportion in this study is rather consistent with the reports from allopolyploids of flowering plants: 12.5% in this study, 8.9–10.9% in Achillea (Chen et al. 2021), 16–18% in Coffea (Bardil et al., 2011), 13.2–30.4% in Gossypium (Flagel and Wendel, 2010) and 5.5–8.9% in Tragopogon (Shan et al., 2020). Tracing the category transition of genes between the F_1K hybrids and the allotetraploids, the transition from the category “No-Change” to “Transgressive Down-regulation” mostly accounted for the increase in the number of genes showing transgressive expression ( Figure 2 ). This result indicates that additional expression changes that had occurred after polyploidization contributed to the higher rate of transgressive expression genes in the allotetraploids. A similar phenomenon including the category transition from “No change” to transgressive expression has been reported in angiosperms, and was likely to be caused by epigenetic modification or cis- and trans-regulatory evolution (Flagel and Wendel, 2010; Hu and Wendel, 2019).

ELD is a well-known phenomenon observed in numerous allopolyploids as demonstrated in several previous studies (Rapp et al., 2009; Flagel and Wendel, 2010; Bardil et al., 2011; Sigel et al., 2019; Chen et al., 2021). In this study, a substantial number of genes showing “K-dominance” or “T-dominance” in the F_1K hybrids was also highly conserved in the allotetraploids: 67.3% in “K-dominance” and 63.3% in “T-dominance” ( Figure 2 ). Therefore, hybridization plays a major role in shaping the pattern of ELD in this study. In many allopolyploid species, unbalanced ELD has been reported, with a tendency for expression levels of more genes to resemble those of one parental species (Rapp et al., 2009; Wu et al., 2018; Chen et al., 2021). Similarly, in the allopolyploid fern Polypodium, almost all genes showing ELD are reportedly biased toward one parental species (Sigel et al., 2019). In this study, however, although the number of genes showing “T-dominance” was slightly higher than those showing “K-dominance” in both the F_1K hybrids (5.6% vs. 4.3%) and the F_1T hybrids (6.0% vs. 3.8%), the allotetraploids showed the same extent in the proportion of genes exhibiting either “T-dominance” and “K-dominance” (5.1% vs. 4.8%). This pattern of ELD in the allotetraploids was shaped by larger outflows from “T-dominance” to other categories than inflows and larger inflows from other categories into “K-dominance” than outflows ( Figure 2 ). This suggests that the hybridization between Phegopteris koreana and P. taiwaniana contributed to “T-dominance”, and the alleviation of “T-dominance” and the activation of “K-dominance” occurred over evolutionary time after polyploidization.

HEB is the phenomenon where F₁ or allopolyploid homoeologous gene pairs exhibit expression ratios that differ from those of orthologous gene pairs of their parent species (Yoo et al., 2013). HEB has been frequently documented in allopolyploid species (Chaudhary et al., 2009; Koh et al., 2010; Edger et al., 2017; Boatwright et al., 2018; Sigel et al., 2019; Boatwright et al., 2021; Chen et al., 2021; Vasudevan et al., 2023). In this study, 21.1–33.3% of the 9,540 gene pairs showed HEB in both the F₁ hybrids and the allotetraploids ( Table 2 ). Importantly, a significant number of genes displaying HEB were shared between the F₁ hybrids and the allotetraploids. More than half of the genes showing HEB in the F₁ hybrids (60.5%) were also identified with HEB in the allotetraploids. The transition from “K-bias” to “T-bias” or vice versa between the F₁ hybrids and the allotetraploids was infrequently observed ( Figure 4 ). These findings demonstrate that HEB primarily originated in the response to hybridization before the occurrence of polyploidization. However, there was also a substantial change in the pattern between the F₁ hybrids and the allotetraploids—the intensification of HEB in the allotetraploids compared to the F₁ hybrids. This intensification was evident not only in the number of genes showing HEB but also in the “magnitude” of bias for both “K-bias” and “T-bias” in the allotetraploids (see Tables 2 , 3 ; Figure 5 ). These changes resulted in the lowest kurtosis in the distribution of the expression level ratio between homoeologs in the allotetraploids ( Table 1 ; Figure 3 ). Therefore, this observation suggests that HEB gradually intensified in the number of genes and the “magnitude” of bias after polyploidization. Unbalanced HEB is commonly reported across allopolyploids (Chaudhary et al., 2009; Koh et al., 2010; Edger et al., 2017; Sigel et al., 2019; Chen et al., 2021; Vasudevan et al., 2023). In this study, there was no apparent bias toward either parent in the F₁ hybrids, and the allotetraploids (1,339 “K-bias” vs. 1,400 “T-bias” in the F_1K hybrids, 1,020 vs. 991 in the F_1T hybrids, and 1,619 vs. 1,558 in the allotetraploids) ( Table 2 ). This pattern of HEB aligns with the previous reports of allopolyploids exhibiting a balanced HEB (Boatwright et al., 2018; Boatwright et al., 2021; Burns et al., 2021; Wu et al., 2021). The balanced HEB obtained in this study contrasts with the previous study focusing on allotetraploid fern, which reported an exceptionally strong HEB towards one parent (744–923 genes vs. 39–100 genes) (Sigel et al., 2019). Therefore, the strong unbalanced HEB observed in the previous study is not a common pattern in transcriptomic change after allopolyploidization among ferns.

The HEB is expected to be partly attributed to cytonuclear gene coordination (Sharbrough et al., 2017), although its significance for gene expression patterns in allopolyploids remains inconsistent among studies (Gong et al., 2012; Gong et al., 2014; Sehrish et al., 2015; Ferreira De Carvalho et al., 2019). In allopolyploid plants, the nuclear genome is contributed from both parental species, but organellar genomes are inherited only from maternal parents, including ferns (Gastony and Yatskievych, 1992; reviewed in Kuo L.-Y. et al., 2018). Therefore, cytonuclear coordination should occur following hybridization, where the expression of nuclear-encoded organelle-targeted genes is biased toward the maternal parent, or the function of the paternal homoeolog is lost (Sharbrough et al., 2017). Given that the maternal parent of Phegopteris decursivepinnta is P. koreana (Fujiwara et al., 2021), it is suggested that cytonuclear interactions might favor higher expressions of organelle-targeted homoeologs derived from P. koreana. In this study, among genes exhibiting extremely strong HEB (“magnitude” > |10|), whereas one out of five K-biased genes and two out of nine T-biased genes in the F_1K hybrids were associated with chloroplast and mitochondria-related functions ( Supplementary Tables S6, S7 ), four out of 13 K-biases genes and four out of 11 T-biased genes in the allotetraploids were detected to be functionally related to organelle ( Supplementary Tables S4, S5 ). Two of four organelle-related K-biased genes detected in the allotetraploids are orthologue to AtRBCX1 (AT4G04330), which functions in RuBisCO assembly (Kolesiński et al., 2011), and to NDU9/B14.5b, a subunit of mitochondrial Complex I (Subrahmanian et al., 2016) ( Supplementary T able S4 ). Consistent with our initial prediction, our results support the idea that cytonuclear coordination favors the higher expression of the maternal homologue in organelle-targeted genes, particularly those directly related to organelle function and structure. Therefore, this finding indicates that the increased HEB in the allotetraploids observed in this study can be partly explained by the concerted action between the duplicated nuclear genome and the increased copy of organelle genomes in cells.

Our results showed that neither ELD nor HEB showed a pronounced bias toward one of the subgenomes in the F₁ hybrids and the natural allotetraploids, Phegopteris decursivepinnata. The pattern observed in P. decursivepinnata can be interpreted as the absence of subgenome dominance when compared with the previous studies (Flagel and Wendel, 2010; Edger et al., 2017; Bird et al., 2018; Sigel et al., 2019; Chen et al., 2021; Glombik et al., 2021; Vasudevan et al., 2023). The absence of subgenome dominance in P. decursivepinnata contrasts with the extreme case of subgenome dominance observed in Polypodium hesperium, the only fern species for which DEG analysis has been conducted to date (Sigel et al., 2019). The difference in transcriptomic consequence of allopolyploidy between the two ferns may be explained by two possible reasons.

Firstly, the age of allopolyploid ferns may reflect the degree of subgenome dominance. In other words, it is likely that subgenome dominance has not extensively proceeded in Phgopteris decursivepinnata due to its relatively recent origin. The previous phylogenetic analysis showed no genetic differentiation from parental species for cpDNA gene (rbcL) and nuclear DNA (pgiC), indicating the recent origin of the allotetraploid species (Fujiwara et al., 2021). However, because the parental species currently do not co-occur due to their ecological differentiation, the allopolyploid speciation is not an ongoing process, but rather more likely to have occurred during the climatic oscillations in the Quaternary, as suggested for the other allopolyploid ferns in Japan (Fujiwara and Watano, 2020; Fujiwara et al., 2022). However, the allopolyploid fern showing strong subgenome dominance, Polypodium hesperium is also considered to have originated very recently during glacial cycles in the Pleistocene (Sigel et al., 2014). Therefore, the difference in the degree of subgenome dominance between the two allopolyploid ferns with recent origins cannot be attributed to the difference in their age since polyploidization. Secondly, the degree of subgenome dominance may be linked to genomic differentiation between parental species, known as parental legacy (Buggs et al., 2014). In other words, the genomic properties inherited from the parental species may predetermine the tendency of subgenome dominance. Indeed, in allopolyploid species that exhibit strong subgenome dominance, the bias toward one subgenome is already evident in the F₁ hybrids and the first generation of synthetic allopolyploids and subsequently tends to be strengthened in the later generations of allopolyploids (Flagel and Wendel, 2010; Edger et al., 2017; Bird et al., 2021; Vasudevan et al., 2023). Increasing evidence suggests that one significant aspect of the genomic difference between parental species that contributes to subgenome dominance is the variation in the abundance of transposable elements (TEs) between them (Bird et al., 2018; Alger and Edger, 2020). TEs can be activated during allopolyploidization, but subsequently, most TEs are rapidly silenced. This process affects the expression of the neighboring coding genes, resulting in a lower expression or silencing in the homoeologous genes of a subgenome possessing more TEs when compared with another subgenome. According to the hypothesis, the observed difference in subgenome dominance between the two allopolyploid ferns may be explained by the difference between the two species in genomic differentiation in TE abundance between the parental species. Future studies will require an investigation into the differences in the genomic composition of these species by focusing on the abundance of TEs.

Although we found the contrasting result of subgenome dominance between two allopolyploid ferns, several lines of genomic evidence can predict that the absence of subgenome dominance may be more common, particularly in homosporous ferns. Unlike angiosperms, ferns are hypothesized to exhibit low TE activity, and this is supported by the old insertion time of LTR retrotransposons (Baniaga and Barker, 2019) and genome size stability over long periods (Schneider et al., 2015; Clark et al., 2016; Fujiwara et al., 2023). Therefore, ferns may maintain a highly similar genomic landscape (Szövényi et al., 2021; Huang et al., 2022), even between species from different genera (Rothfels et al., 2015), and the relatively less divergent genome in ferns could potentially generate subtle or absence of subgenome dominance in most homosporous fern allopolyploids, as observed in this study. Furthermore, subgenome dominance is hypothesized to trigger a biased fractionation of the recessive subgenome, leading to rapid diploidization and genome-downsizing in polyploid genomes (Cheng et al., 2018; Wendel et al., 2018; Alger and Edger, 2020). In ferns, however, it is traditionally believed that the diploidization process is slower when compared with angiosperms (Barker, 2013). Instead of the physical loss of genes, the diploidization in ferns is hypothesized to involve gene silencing, leading to slow rates of chromosome reduction and genome downsizing (Haufler, 2014). The limited empirical evidence based on genome sequencing also supports this hypothesis on slow genome evolution in homosporous ferns (Huang et al., 2022; Zhong et al., 2022), with one exception (Marchant et al., 2022). Therefore, the observed absence of subgenome dominance in Phegopteris decursivepinnata is consistent with the slow genome downsizing in ferns and could partially explain the process of post-polyploid evolution seen in the broad range of homosporous ferns. Also in homosporous lycophyte, a lineage of land plants considered to exhibit a similar genome evolution as homosporous ferns (Sessa and Der, 2016), a recent finding from ancient allotetraploid Huperzia showed a limited subgenome dominance and slow genome evolution following polyploidization (Li et al., 2023), which support our prediction. To test the validity of our prediction, we still need to examine transcriptomic change and the pattern of genomic differentiation including TE abundance in more diverse allopolyploid fern species.

In this study, we successfully demonstrated the dynamics of gene expression change between the F₁ hybrids and the allotetraploids of Phegopteris fern, highlighting the absence of subgenome dominance in the allotetraploid. Nevertheless, it is essential to note several limitations in our study design and technical approaches that need to be addressed in future investigations.

Firstly, one of the limitations is the relatively limited number of genes examined in this study. We targeted 9,540 genes co-expressed among the parental species, the F₁ hybrids, and the allotetraploids. This represents less than half of the total number of protein-coding genes reported in homosporous ferns (36,857 genes in Ceratopteris, Marchant et al., 2022; 31,244 genes in Adiantum, Fang et al., 2022). However, this is because the necessity of focusing on co-expressed genes arose from the absence of available genome sequence data for Phegopteris fern and the specific goal of tracing the transition of genes exhibiting ELD and HEB from the F₁ hybrids to the allotetraploids. Although our study successfully provided a robust overview of the transcriptomic changes in the allotetraploids as evidenced by the consistency across three datasets with TPM filtering thresholds of > 0, > 0.5, and > 1.0 ( Supplementary Figures S4–S7 ), it should be noted that our results are solely based on co-expressed genes and thus might have led to the oversight of critical changes, such as gene silencing, and reactivation of genes suppressed in the parental species. To comprehensively examine genome-wide expression changes following allopolyploidization, future studies should consider analyzing transcriptome data using parental genome sequences as references.

Secondly, this study is solely based on gene expression in leaf tissues, possibly overlooking tissue-specific expression differences. Previous studies comparing homoeolog gene expression among different tissues reported substantial differentiation of homoeolog expression bias in direction and expression level for each gene among the examined tissues (Adams et al., 2003; Buggs et al., 2010; Huang et al., 2022). On the other hand, recent studies focusing on the whole transcriptome-wide pattern of homoeolog gene expression suggest that the overall trend of homoeolog expression bias is highly conserved among different tissues (Edger et al., 2017; Griffiths et al., 2019; Chen et al., 2021; Vasudevan et al., 2023). Therefore, while the trend obtained in this study could be consistent among other tissues, we might have missed key changes in homoeolog gene expression occurring in different tissues, such as tissue-specific silencing and tissue-specific homoeolog bias, as reported in other studies. Thus, in future studies, it is required to examine gene expression patterns using different tissue types of Phegopteris species.

Lastly, the pattern of gene expression in the allotetraploids might change in response to environmental conditions (Shimizu, 2022). In Cardamine, studies have shown that allopolyploids in each dry and humid environment tend to exhibit gene expression patterns resembling those of each parental species specialized for each environment, which enables the allopolyploids to obtain adaptive plasticity to diverse niches (Shimizu-Inatsugi et al., 2016; Akiyama et al., 2021). Phegopteris decursivepinnata has a wider distribution than those of its parental species (Fujiwara et al., 2021) and thus is expected to have a higher tolerance to a wide range of environmental conditions than its parental species. Given that, it is likely that P. decursivepinnata could also show different patterns of homoeolog gene expression from that observed in this study according to different environmental conditions. Thus, it is essential to trace the change in the pattern of gene expression in response to different environmental conditions in this species.