Spores of A. flabellulatum were collected from the rubber forest at Ma’an Mountain (109.519°E, 19.501°N) in Danzhou, Hainan Province, China, as described by Cai et al. (2022b). The spores were sterilized with 0.1% HgCl₂ for 5 min, washed with sterile water five times, and then plated on 1/4 MS solid medium (3% w/v sucrose and 0.7% w/v agar). The plated spores were exposed to meticulously controlled illumination, with a photosynthetic photon flux density (PPFD) level of 7.1 μmol m^-2 s^-1, and maintained under an uninterrupted 24-hour photoperiod. Upon spore germination, young gametophytes with diameters ranging from 1 to 2 mm were transferred to fresh 1/4 MS solid medium, also containing 3% w/v sucrose and 0.7% w/v agar. Subsequently, the gametophytes were divided into seven distinct experimental groups, each cultivated under varying PPFD conditions: 0, 0.1, 1.4, 7.1, 14.4, 73.6, and 145.3 μmol m^-2 s^-1, respectively ( Supplementary Table 1 , Supplementary Figure 1 ). These cultivations were carried out over a continuous 40-day period under an uninterrupted 24-hour photoperiod. The perimeter and projected area of the gametophytes were quantified using Image J software (v1.8.0). Each treatment was replicated at least eight times, and the cultivation temperature was maintained at 25°C. Statistical analysis was conducted using the least significant difference (LSD) test.

The gametophytes of A. flabellulatum were cultivated under diverse light conditions: 0, 0.1, 7.1, and 145.3 μmol m^-2 s^-1, for a period of 12 days, maintaining an uninterrupted 24-hour photoperiod. Subsequently, the four sample groups, denoted as Af0, Af0.1, Af7.1, and Af145.3, were collected for transcriptome sequencing. Three biological replicates were established for each treatment, resulting in a total of 12 samples. Approximately 200 mg of gametophytes were collected from each sample, and total RNA was extracted using the CTAB method. mRNA was then enriched using oligo dT magnetic beads, fragmented, and reverse-transcribed into cDNA for library construction. Sequencing was performed on the DNBSEQ platform. The raw reads were filtered using SOAPnuke software (v1.5.2) to remove reads with adapters, N content exceeding 10%, and more than 50% of bases with a quality score below 15. The clean reads were aligned to the full-length transcriptome of A. flabellulatum gametophytes using Bowtie2 (v2.2.5), and gene expression levels were calculated using RSEM (v1.2.8). These steps were executed by the Beijing Genomics Institute (BGI). Finally, we conducted sequencing saturation analysis, correlation analysis, and principal component analysis using the method described by Cai et al. (2021).

The study involved three comparison groups: Af0 vs Af0.1, Af0.1 vs Af7.1, and Af7.1 vs Af145.3. Genes with average expression levels less than 0.5 in each comparison group were removed, and differentially expressed genes (DEGs) were screened using the criteria of |log₂ FoldChange| ≥ 1 and Q-Values < 0.001. The expression levels of all DEGs in the 12 samples were used to construct gene co-expression modules using the WGCNA package in R software, with “sft$powerEstimate=13” and “mergeCutHeight=0.25”. Modules that were strongly correlated with the phenotypes and light conditions of A. flabellulatum gametophyte were selected for further analysis.

The module genes were annotated and classified using Blast2GO (v2.5.0) for Gene Ontology (GO, http://www.geneontology.org/) and Blastx (v2.2.23) and Diamond (v0.8.31) tools for KEGG (http://www.genome.jp/kegg/) annotation and classification, respectively (Cai et al., 2022b). Enrichment analysis was performed using the Phyper function package in the R software. BGI completed the GO and KEGG annotations. Enrichment pathways related to the gametophyte phenotype and light conditions of A. flabellulatum were selected for further analysis and discussion.

In this study, the identification of open reading frames (ORFs) from gene sequences was conducted using getorf (EMBOSS: 6.5.7.0). The resulting ORFs were aligned to transcription factor protein domains (http://plantregmap.gao-lab.org/) using hmmsearch (v3.0) (Tian et al., 2020). The genes’ potential for encoding transcription factors were then identified based on the transcription factor family characteristics outlined in PlantTFDB (http://planttfdb.gao-lab.org/) (Jin et al., 2017) and classified according to their respective transcription factor families.

We used Cytoscape software (v3.9.1) to visualize the co-expression networks of flavonoid biosynthesis genes, as well as MYBs and bHLHs in the green, brown, yellow, and turquoise modules, following the method described by Shannon et al. (2003). In the visualization, nodes represent genes, and the size of each node is positively correlated with its degree. Lines represent the weight of co-expression between genes, with thicker and darker lines indicating higher weights.

Fifteen genes were chosen for validation by quantitative real-time polymerase chain reaction (qRT-PCR) ( Supplementary Table 2 ). A total of 12 samples were used, with three technical replicates per sample. To achieve relative quantification, we selected isoform_78364 (Elongation Factor Ts), which exhibited stable expression across all samples, as the internal control gene, following the method described by Pfaffl (2001). Relative quantification was carried out using the 2^-ΔΔCt method, and significance analysis was performed using the LSD method.

The perimeters and projected areas of gametophytes grown under different PPFD levels were measured. The results showed that gametophytes are highly sensitive to low light, with their perimeter increasing significantly under a PPFD as low as 0.1 μmol m^-2 s^-1. Below 7.1 μmol m^-2 s^-1, the perimeter and/or projected area of gametophytes increased with increasing PPFD levels and reached a maximum at 7.1 μmol m^-2 s^-1. When the PPFD levels ranged between 7.1 μmol m^-2 s^-1 and 73.6 μmol m^-2 s^-1, the average value of the perimeter and projected area of gametophytes decreased slightly, but no significant differences were observed. However, when the PPFD reached 145.3 μmol m^-2 s^-1, the perimeter and projected area of gametophytes decreased significantly ( Figure 1 ). Therefore, we conclude that 7.1 μmol m^-2 s^-1 is the optimal PPFD for the growth of A. flabellulatum gametophytes, while both PPFD levels below 7.1 μmol m^-2 s^-1 and at 145.3 μmol m^-2 s^-1 can induce light stress on the gametophytes.

We selected four treatments with PPFD levels of 0, 0.1, 7.1, and 145.3 μmol m^-2 s^-1 (represented by Af0, Af0.1, Af7.1, and Af145.3, respectively) to explore the genes involved in the response of A. flabellulatum gametophytes to light using transcriptome sequencing, based on previous experimental results. Three biological replicates were performed for each treatment, resulting in a total of 12 samples that were sequenced using the DNBSEQ platform. The full-length transcript sequences of A. flabellulatum gametophytes reported by Cai et al. (2022a) were used as reference genes for relative quantification.

The raw read counts of each sample were greater than 42 M, and after filtering, each sample had at least 41 M clean reads, with a ratio greater than 96.17% and Q20 and Q30 values higher than 94% and 87%, respectively ( Supplementary Table 3 ). Sequencing saturation analysis indicated that the gene identification ratio plateaued when the reads number exceeded 50×100 K, suggesting that the sequencing data had reached saturation ( Supplementary Figure 2A ). The gene expression level distributions of the 12 samples were analyzed, and the median gene expression levels were all ≥0.97 ( Supplementary Figure 2B ). Pearson correlation coefficients between biological replicates were ≥0.91 ( Figure 2A ), and the four treatments were separated from each other on the plane composed of PC1 and PC3 in the principal component analysis ( Figure 2B ). These results confirm the reliability and high quality of the sequencing data, which can be used for subsequent DEGs screening and analysis. The sequencing data have been deposited in the NCBI SRA database under accession number PRJNA868202.

A total of 248,028 genes were detected among the 12 samples. Three comparison groups were set: Af0 vs Af0.1, Af0.1 vs Af7.1, and Af7.1 vs Af145.3. In each comparison group, we filtered out genes with an average expression level of less than 0.5 for both treatments and then selected DEGs with |log₂ FoldChange| ≥ 1 and Q-Values < 0.001. This resulted in 1,059, 5,536, and 10,152 DEGs, respectively. The Af7.1 vs Af145.3 comparison had the most DEGs, with 7,692 upregulated and 2,460 downregulated genes, followed by Af0.1 vs Af7.1, which had 1,864 upregulated and 3,672 downregulated genes. The Af0 vs Af0.1 comparison had the fewest DEGs ( Supplementary Tables 4 – 6 ; Figures 3A–C ). After removing duplicates, a total of 15,874 DEGs were obtained ( Figure 3D ).

For these 15,874 DEGs, we performed weighted gene co-expression analysis and obtained 12 gene co-expression modules ( Figure 3E ). The turquoise module contained 8,796 DEGs and was positively correlated with “PPFD levels,” while the yellow (746 DEGs) and brown (875 DEGs) modules were positively correlated with “gametophyte area.” The green module (387 DEGs) was positively correlated with “presence or absence of light” ( Supplementary Table 7 ). In contrast, the black module (191 DEGs) was negatively correlated with “presence or absence of light” ( Supplementary Table 7 , Supplementary Figure 3 ). Based on these results, we selected DEGs from these five modules for further analysis.

We performed GO annotation on a total of 10,995 DEGs from the five gene modules mentioned earlier. Of these, 8,033 were successfully annotated, accounting for 73.06% of the total. These annotations covered three major functional categories: “Biological process,” “Cellular component,” and “Molecular function,” which contained 4,911, 5,451, and 6,424 genes, respectively. In the “Biological process” category, all five modules had the highest number of genes annotated in the “cellular process” term, while in the “Cellular component” category, it was the “cellular anatomical entity.” In the “Molecular function” category, the black and brown modules had the highest number of genes annotated in the “binding” term, while the green, yellow, and turquoise modules were most enriched in the “catalytic activity” term ( Supplementary Tables 8 , 9 ; Figure 4A ).

We present the top 21 enriched GO terms with Q-values < 0.01. Firstly, the green module was enriched in at least 8 terms related to light signaling (8/21). Among them, “photoreceptor activity” had the smallest Q-value (Q=1.67E-07) and the most genes (8) ( Supplementary Table 10 , Supplementary Figure 4A ). This indicates that the DEGs in the green module play an important role in regulating light signaling during the development of A. flabellulatum’s gametophyte. Secondly, most of the terms in the black and brown modules were related to photosynthesis. For example, “thylakoid” had the smallest Q-value (Q=1.63E-28) and the most genes (52) in the black module, while “regulation of photosynthesis, dark reaction” had the highest enrichment rate (0.11) ( Supplementary Table 10 , Supplementary Figure 4B ). In the brown module, “photosystem II” had the smallest Q-value (Q=2.35E-115) and the most genes (166), while “photosynthetic acclimation” had the highest enrichment rate (0.125) ( Supplementary Table 10 , Supplementary Figure 4C ). This indicates that the DEGs in the black and brown modules mainly participate in regulating photosynthesis during the development of A. flabellulatum’s gametophytes. Additionally, the two terms with the smallest Q-value in the yellow module were “flavonoid biosynthetic process” (Q=3.61E-09) and “flavonoid metabolic process” (Q=3.87E-09), and they both had the most genes (12) ( Supplementary Table 10 , Supplementary Figure 4D ). The terms enriched in the turquoise module were mainly related to cell damage, such as “cell killing” and “killing of cells of other organism,” and they had the most genes (36) ( Supplementary Table 10 , Supplementary Figure 4E ). This indicates that the DEGs in the yellow and turquoise modules mainly participate in resisting light stress caused by strong light.

We performed KEGG annotation on 10,995 DEGs, of which 5,405 genes were successfully annotated with an annotation rate of 49.16%. These genes were classified into five major categories, including “Cellular processes”, “Environmental information processing”, “Genetic Information Processing”, “Metabolism”, and “Organismal Systems”. Among these categories, the “Metabolism” category had the most DEGs, with a total of 3,363 genes. The five gene modules were annotated to the most genes in the “Transport and catabolism” (354 genes) of “Cellular processes”, “Signal transduction” (235 genes) of “Environmental information processing”, “Global and overview maps” (2,858 genes) of “Metabolism”, and “Environmental adaptation” (207 genes) of “Organismal Systems”. In the “Genetic Information Processing” category, the black and turquoise modules had the most genes annotated in “Translation”, while the green, brown, and yellow modules were involved in “Folding, sorting and degradation”. Importantly, in the “Metabolism” category, all five gene modules had a relatively high number of genes annotated in “Energy metabolism” (1,088 genes) and “Carbohydrate metabolism” (956 genes), indicating that the regulation of these two metabolisms by light is also crucial ( Supplementary Tables 11 , 12 ; Figure 4B ).

We presented the top 10 KEGG pathways with the smallest Q-values in enrichment analysis ( Supplementary Figure 5 ; Figure 5 ). Firstly, the green, black, and yellow gene modules were all enriched in the “Circadian rhythm-plant” pathway, with the green module having the smallest Q-value and the largest number of genes (14) in this pathway ( Supplementary Table 13 , Supplementary Figure 5A ). Secondly, the black, yellow, brown, and turquoise gene modules were all enriched in multiple pathways related to photosynthesis, such as “Photosynthesis-antenna proteins” and “Photosynthesis” ( Supplementary Table 13 ). The black module had the smallest Q-value (1.63E-20), highest enrichment rate (0.009240924), and the largest number of genes (28) in the “Photosynthesis” pathway ( Supplementary Table 13 , Supplementary Figure 5B ). In the brown module, the “Photosynthesis-antenna proteins” pathway had the smallest Q-value (1.08E-77), highest enrichment rate (0.041363336), and the largest number of genes (125) ( Supplementary Table 13 ; Figure 5A ). Additionally, the brown, yellow, and turquoise gene modules were enriched in the “Flavonoid biosynthesis” pathway, which is involved in photoprotection ( Supplementary Table 13 , Supplementary Figure 5C ; Figure 5 ). In summary, our KEGG analysis results also indicate that light plays an essential regulatory role in light signal transduction, photosynthesis, and photoprotection during gametophyte development in A. flabellulatum.

We predicted transcription factor families for the 10,995 DEGs and identified 34 families, with mTERF, MYB, and bHLH families having the highest number of DEGs at 13 each. The turquoise module had the most predicted transcription factor families (30) and individual transcription factors (95), followed by the yellow module with 13 families and 17 individual factors ( Supplementary Table 14 , Supplementary Figure 6 ).

As the “Photosynthesis-antenna proteins” and “Photosynthesis” pathways are significantly enriched in multiple gene modules, we first analyzed the expression of DEGs in these two pathways. The “Photosynthesis-antenna proteins” pathway includes Light-harvesting complex I (LHC I) and Light-harvesting complex II (LHC II), with LHC I containing five subunits (Lhca1-Lhca5) and LHC II containing seven subunits (Lhcb1-Lhcb7). Results showed that among the five gene modules, the “Photosynthesis-antenna proteins” pathway had a total of 327 DEGs, including the 12 subunits mentioned above, mainly belonging to the brown (125) and turquoise (168) modules. These genes were generally expressed at the highest levels under a PPFD of 7.1 μmol m^-2 s^-1 (brown module) or 145.3 μmol m^-2 s^-1 (turquoise module) ( Supplementary Table 15 ; Figure 6 ). In addition, the “Photosynthesis” pathway includes the Photosystem II complex (PSII complex), cytochrome b₆/f complex, Photosystem I complex (PSI complex), photosynthetic electron transfer chain, and Adenosine triphosphate synthase (ATP synthase). Results showed that among the five gene modules, the “Photosynthesis” pathway had a total of 370 DEGs, including 30 subunits, such as 11 subunits of the PSII complex (PsbD, PsbO, PsbP, PsbQ, PsbR, PsbS, PsbW, PsbY, PsbZ, Psb27, Psb28), one subunit of the cytochrome b₆/f complex (PetC), 10 subunits of the PSI complex (PsaB, PsaD, PsaE, PsaF, PsaG, PsaH, PsaK, PsaL, PsaN, PsaO), 4 subunits of the photosynthetic electron transfer chain (PetE, PetF, PetH, PetJ), and 4 subunits of ATP synthase (gamma, delta, a, b). These DEGs were also mainly distributed in the brown (76) and turquoise (231) modules ( Supplementary Table 16 , Figure 7 ).

We also analyzed the DEGs involved in Calvin cycle (141), chlorophyll (61), and carotenoid (28) biosynthesis, and found that these DEGs were mainly classified into the turquoise module ( Supplementary Tables 17 - 19 , Supplementary Figures 7 - 9 ).

Moreover, mitochondrial transcription termination factors (mTERFs) located in plastids play important regulatory roles in the development and function of chloroplasts. AfmTERFs were detected in the brown, yellow, and turquoise modules, with 1, 1, and 11 members, respectively. Therefore, we constructed a phylogenetic tree of mTERFs from A. flabellulatum and Arabidopsis. The results showed that four AfmTERFs (isoform_43473, isoform_47623, isoform_50012, and isoform_240354) in A. flabellulatum were homologous to AtmTERF6, one (isoform_20996) was homologous to AtmTERF9, and two (isoform_83065 and isoform_91611) were homologous to AtmTERF12 ( Supplementary Table 20 , Supplementary Figure 10 ). Interestingly, all seven AfmTERFs belonged to the turquoise module.

Based on the analysis in the previous sections, we found that the flavonoid biosynthesis pathway is significantly enriched in multiple gene modules, particularly in the brown, yellow, and turquoise modules. Therefore, we analyzed the expression of DEGs in this pathway. Our results show that there are 210 DEGs encoding 13 enzymes, including CHS, CHI, F3H, and FLS, which are key enzymes in the flavonoid biosynthesis pathway. These DEGs mainly belong to the turquoise module ( Supplementary Table 21 ; Figure 8 ). Thus, we suggest that light, especially strong light, strongly upregulates the expression of genes encoding flavonoid biosynthesis enzymes in the gametophyte of A. flabellulatum.

The MYB and bHLH families of transcription factors are important in plant responses to abiotic stress. These two types of transcription factors were detected in four modules, with the highest number found in the turquoise module. Therefore, we performed co-expression analysis of flavonoid biosynthesis genes and MYB/bHLH transcription factor genes in A. flabellulatum. Our results showed that the encoding genes of flavonoid biosynthesis enzymes, CHS, DFR, HCT, and CHI, have a wide range of co-expression relationships with MYBs and bHLHs, with 39, 20, 20, and 19 genes in the turquoise module, respectively ( Supplementary Table 22 ; Figure 9 ). Furthermore, we also found co-expression relationships between MYBs and bHLHs in the turquoise module.

To validate the accuracy of the transcriptomic data, we selected 15 DEGs for qRT-PCR analysis. These included two genes from the “Photosynthesis-Antenna Proteins” pathway (Lhca3 and Lhcb1/Lhcb2), three genes from the “Photosynthesis” pathway (PsaG, PsaL, and PetE), nine genes from the “Flavonoid Biosynthesis” pathway (CHSs, CHIs, F3H, ANS, DFR, and HCTs), and Cryptochrome (CRY). Our qRT-PCR results demonstrated that the expression patterns of these 15 genes were consistent with the transcriptomic data, confirming the reliability of our transcriptomic analysis ( Supplementary Table 23 ; Figure 10 ).