Full-length transcriptome analysis of Adiantum flabellulatum gametophyte

Plant materials

The spores of A. flabellulatum were collected from the rubber forest in Ma’an mountain, Danzhou City, Hainan Province, China (109.519°E, 19.501°N). The mature spores were sterilized and inoculated in 1/4 MS medium to form young gametophytes after germination. Young gametophytes with diameters of 1–2 mm were selected and transferred into new 1/4 MS medium. These young gametophytes were then placed in illuminations of 0, 10, 500, 1,000 and 10,000 lx, respectively, and cultured at 25 °C for 12 days. Three biological replicates were set in the above 5 illuminations, respectively, with a total of 15 samples. Approximately 200 µg of gametophytes from each sample were selected for RNA extraction.

RNA extraction, Library construction and Sequencing

Total RNA (more than 40 ng/µL) was extracted from the above 15 samples using the CTAB method, and the mRNAs were purified and mixed to construct a cDNA library (Yu et al., 2021). Raw reads were obtained on the PacBio Sequel platform. After reads of insert (ROI) were generated from raw reads, clustering and correction (Chin et al., 2013) were performed using the Quiver algorithm to obtain high-quality full-length consensus reads and then the transcripts were obtained after redundancy removal (using the Cd-hit program, Li & Godzik, 2006). RNA extraction, library construction and sequencing were performed by the Beijing Genomics Institute (BGI).

Coding sequences prediction

The Transdecoder software (v3.0.1, Haas et al., 2013) was used to identify the candidate coding regions in the transcripts. The longest open reading frames were extracted, and the coding sequences (CDSs) were predicted using the SwissProt (http://ftp.ebi.ac.uk/pub/databases/swissprot, Boutet et al., 2007) and protein families (Pfam, http://pfam.xfam.org, Finn et al., 2016) databases.

Simple sequence repeats detection

The MISA software (v1.0) (Thiel et al., 2003) was used to detect the simple sequence repeats (SSRs). The detection parameters were mono-nucleotide, di-nucleotide, tri-nucleotide, tetra-nucleotide, penta-nucleotide and hexa-nucleotide repeats at least 12, 6, 5, 5, 4 and 4 times, respectively. The Primer3 software (v2.2.2) (Untergasser et al., 2012) was used to design primers for the SSRs.

Functional annotation

BLASTN was used for NT (https://blast.ncbi.nlm.nih.gov/Blast) annotations; Blastx (BLAST, v2.2.23, Altschul et al., 1990) or Diamond (v0.8.31) (Buchfink, Xie & Huson, 2015) were used for NR (https://blast.ncbi.nlm.nih.gov/Blast.cgi), KEGG (http://www.genome.jp/kegg), KOG (http://www.ncbi.nlm.nih.gov/KOG) and SwissProt annotations; and Blast2GO (v2.5.0) (Conesa et al., 2005) was used for GO (http://geneontology.org) annotations. After translating CDSs into protein sequences, the protein sequences were uploaded to iTAK (http://itak.feilab.net/cgi-bin/itak/index.cgi) for the prediction of transcription factors (TFs), transcription regulators (TRs) and protein kinases (PKs).

Long non-coding RNAs prediction

The Coding Potential Calculator (CPC) (v0.9-r2) (Kong et al., 2007), txCdsPredict (http://hgdownload.soe.ucsc.edu/admin/jksrc.zip) (Zweig et al., 2008) and the Coding-Non-Coding Index (CNCI) (https://github.com/www-bioinfo-org/CNCI) (Sun et al., 2013) were used to score transcripts other than CDSs, and the scoring thresholds (CPC_threshold = 0, txCdsPredict_threshold = 500, CNCI_threshold = 0) were set to determine long non-coding RNAs (lncRNAs). Pfam was used to predict the protein-coding potential of a transcript. The transcript that could not be predicted in Pfam is potential lncRNA. If the results with three of the above four methods were consistent, the transcript was finally confirmed as lncRNA.

Quantitative real-time polymerase chain reaction

In order to verify the reliability of the full-length transcript sequences and the annotation information, we selected nine chlorophyll synthase genes and used the method applied by Cai et al. (2021) for qRT-PCR (Table S1). Each of the four treatments (0 lx, 10 lx, 500 lx and 10,000 lx) were detected with three biological replicates for a total of 12 samples. Three parallel tests were performed for each sample, and the isoform_55882 (Ubiquitin C) was used as a reference gene. Relative quantification was performed using the 2^−ΔΔCt method (Pfaffl, 2001; Cai et al., 2021).

PacBio SMRT sequencing data

In this study, a PacBio ISOSeq library was established, and the full-length transcriptome sequencing was performed on the PacBio Sequel platform. The total base of polymerase reads was 69.87 Gb, and the total number of reads was 923,898. The average length, the longest sequence and the N50 of polymerase reads were 75,625.26 bp, 303,495 bp and 129,782 bp, respectively (Fig. 1A). After removing the adapters, 16,221,787 subreads were obtained, with the total base of 68.72 Gb. The average length, the longest sequence and the N50 of subreads were 4,236.18 bp, 218 714 bp and 4,762 bp, respectively (Fig. 1B).

After subreads from the same circular molecule were combined into circular consensus sequences (CCSs), a total of 809,367 CCSs were obtained. The average length and average quality were 4,934 bp and 0.97, respectively (Fig. 1C). The CCSs were split into full-length, non-full-length, chimeric and non-chimeric sequences according to whether 5′ primer, 3′ primer and ploy(A) tail were detected. There were 2,762,613 full-length non-chimeric (FLNC) sequences, with an average length of 998 bp and an average quality of 0.98 (Fig. S1).

The FLNC sequences from the same copy were grouped into clusters, and the Arrow algorithm (Hong et al., 2020) was used to correct the error. After removing the low-quality sequences, 2,608,705 consensus reads were obtained with an average quality of 0.99 (Fig. S2). After removing redundancy, 354,228 transcripts were finally obtained (Table 1, Fig. 1D). The CCS and transcript data have been uploaded to the NCBI Sequence Read Archive (SRA) database (https://www.ncbi.nlm.nih.gov/sra), and the accession numbers are PRJNA774117 and PRJNA733457, respectively.

CDS prediction

CDS prediction was performed on transcripts, and a total of 231,705 CDSs were obtained. The total length was 194,628,012 bp, and the N50 was 1,056 bp (Table 2). There were 164,535 CDSs with a length of 100–1,000 bp, accounting for 71.0% of the total. There were 57,009 CDSs of 1,000–2,000 bp, 8,722 CDSs of 2,000–3,000 bp and 1,439 CDSs of ≥ 3,000 bp, accounting for 24.6%, 3.8% and 0.6% of the total, respectively (SRA database; Fig. 2A). The CDSs have been uploaded to the SRA database, and the accession number is PRJNA777740.

SSR detection

Through SSR detection of 354,228 transcripts, 138,995 SSRs were obtained from 81,304 transcripts, with a total length of 3,875,607 bp (Table S2). Of these, 50,482 transcripts contained one SSR locus and 30,822 transcripts contained two or more SSR loci (Table 3). Among SSR repeat types, di-nucleotide repeats (90,741) accounted for the highest proportion (65.28%). The second was the mono-nucleotide repeats (35,856) accounting for 25.8% of total repeats and the third was tri-nucleotide repeats (10,528) accounting for 7.56% of repeats. Tetra-nucleotide, penta-nucleotide and hexa-nucleotide repeats accounted for the minority at 0.40%, 0.25% and 0.69% of repeats, respectively. Forty-one thousand three hundred sixty-eight SSRs were present in compound formation. SSRs with six repeat units (19,697, 14 17%) were the most common, followed by seven repeat units (12,575, 9.05%), eight repeat units (9,806, 7.05%) and 12 repeat units (9,272, 6.67%) (Table 4). Among the SSRs with di-nucleotide repeats, AG/CT (68,289) was the most common type, and AT/AT (824) was the least. Among the SSRs with tri-nucleotide repeats, ATC/GAT (2,375) was the most common type and AAT/ATT (51) was the least common (Fig. 2B). Twenty-nine thousand five hundred eighty-nine pairs of primers were designed and screened by Primer3. The primer lengths were 18–28 bp, the GC contents were 20.00–72.22%, the annealing temperatures were 57–63 °C, and the product lengths were 80–276 bp (Table S3).

Functional annotation

Three hundred fifty-four thousand two hundred twenty-eight transcripts were annotated by seven databases. There were 251,501, 197,474, 193,630, 194,639, 195,956, 113,069 and 197,883 transcripts annotated by the NR, KEGG, KOG, Swissprot, Pfma, NT and GO databases (Table S4), respectively. Most of them were annotated by the NR database, accounting for 71.00%. The NT database annotated the least, accounting for 31.92% (Table 5). Fifty-five thousand seven hundred ninety-four transcripts were annotated by all the seven databases combined, accounting for 15.75% of the transcripts. There were 125,336 transcripts annotated by NR, KEGG, KOG, SwissProt and Pfam databases together, accounting for 35.38% (Fig. 3A).

NR annotation

The proportion of transcripts annotating to homologous species was counted from the NR annotation results. Among them, the proportions of transcripts annotating to Marchantia polymorpha, Selaginella moellendorffii, Physcomitrella patens and Picea sitchensis were 59.31%, 10.98%, 9.45% and 8.70%, respectively, and the proportion of transcripts annotating to other known species was 11.56% (Table S5; Fig. 3B).

KEGG annotation

A total of 197,474 transcripts were annotated by the KEGG database, of which 112,839 were annotated to 137 pathways (Table S6). The KEGG pathways were divided into five categories: “cell processing,” “environmental information processing,” “genetic information processing,” “metabolism” and “organismal system.” For “metabolism,” “global and overview maps” had the greatest number of transcripts (55,321), accounting for 28.01%, followed by “carbohydrate metabolism,” which had 21,432 transcripts, accounting for 10.85% (Fig. 4A).

KOG annotation

There were 193,630 transcripts annotated by the KOG database, which could be divided into 25 groups according to KOG functional classification. Among them, the transcripts belonging to “general function prediction only” were the most common with a total of 48,175, accounting for 24.88%. There were 26,065 “signal transduction mechanisms” and 23,151 “posttranslational modifications protein turnover chaperones,” accounting for 13.46% and 11.96%, respectively (Fig. 4B).

GO annotation

A total of 197,883 transcripts were annotated by the GO database. GO domains can be divided into three categories: “biological processes,” “cell components” and “molecular functions.” A large number of transcripts in “biological processes” were mainly involved in the “cellular process” (94,069) and the “metabolic process” (81,073). The category “cell components” mainly consisted of transcripts involved in “cellular anatomical entity” (129,811) and “intracellular” (82,674). The “Cell components” category mainly consisted of transcripts involved in “binding” (95,413) and “catalytic activity” (95,176) (Fig. 4C).

Prediction of TFs, TRs and PKs

TFs play an important role in gene expression regulation. A total of 5,749 CDSs were predicted belonging to 64 TF families. Among them, the predicted sequences of C3H family were the most common (614), followed by bHLH (449), STAT (1) was the least common (Table 6; Table S7; Fig. 5). In addition, 2,214 sequences belonged to 24 TR families. Among them, the SET family had the most sequences (251), followed by the GNAT family (248), and the two families that had the least number of sequences were RB (4) and MED7 (4) (Table 6; Table S8; Fig. 6).

PKs regulate protein activity and amplify signals to induce a cellular response during plant signal transduction (Ardito et al., 2017). A total of 4,950 CDSs were predicted belonging to the PK families. Among them, 833 sequences belonged to leucine rich-repeat receptor-like kinases (LRR-RLKs) families and 746 sequences belonged to receptor-like cytoplasmic kinase (RLCK) families, accounting for 16.83% and 15.07% of the total PKs, respectively (Table 6; Table S9; Fig. 7).

lncRNAs prediction

106,931, 118,406, 99,672 and 115,438 lncRNAs were predicted using CPC, txCdsPredict, CNCI tools and the pfam database, respectively (Table S10). Among them, 85,487 were predicted using the 4 methods together, and 2,105, 7,652, 1,839 and 14,710 were predicted by three of the four methods alone (Fig. 8). Finally, 111,793 lncRNAs were obtained, accounting for 31.56% and 91.24% of all transcripts and non-CDSs, respectively.

Quantitative real-time polymerase chain reaction

Chlorophyll is an important pigment in plant photosynthesis, and its synthesis enzyme gene expressions are regulated by light (Lee et al., 2005; Masuda et al., 2003; Matsumoto et al., 2004; Okamoto et al., 2020; Tzvetkova-Chevolleau et al., 2007). In order to verify the reliability of the full-length transcript sequences and the annotation information, we selected nine chlorophyll synthase genes according to KEGG annotation for qRT-PCR, including HemA, HemL, HemB, CHLH, CHLD, CRD, DVR, POR and CAO. The results showed that the expressions of the above nine genes increased at first and then decreased with the increase of illumination (Table S11; Fig. 9). It can be seen that the chlorophyll synthase genes in gametophyte of A. flabellulatum are also regulated by light. This also shows that the data of the full-length transcripts has strong reliability.