Co-expression network analyses of anthocyanin biosynthesis genes in Ruellia (Wild Petunias; Acanthaceae)

Background Anthocyanins are major pigments contributing to flower coloration and as such knowledge of molecular architecture underlying the anthocyanin biosynthetic pathway (ABP) is key to understanding flower color diversification. To identify ABP structural genes and associated regulatory networks, we sequenced 16 transcriptomes generated from 10 species of Ruellia and then conducted co-expression analyses among resulting data. Results Complete coding sequences for 12 candidate structural loci representing eight genes plus nine candidate regulatory loci were assembled. Analysis of non-synonymous/synonymous (dn/ds) mutation rates indicated all identified loci are under purifying selection, suggesting overall selection to prevent the accumulation of deleterious mutations. Additionally, upstream enzymes have lower rates of molecular evolution compared to downstream enzymes. However, site-specific tests of selection yielded evidence for positive selection at several sites, including four in F3'H2 and five in DFR3, and these sites are located in protein binding regions. A species-level phylogenetic tree constructed using a newly implemented hybrid transcriptome–RADseq approach implicates several flower color transitions among the 10 species. We found evidence of both regulatory and structural mutations to F3′5'H in helping to explain the evolution of red flowers from purple-flowered ancestors. Conclusions Sequence comparisons and co-expression analyses of ABP loci revealed that mutations in regulatory loci are likely to play a greater role in flower color transitions in Ruellia compared to mutations in underlying structural genes. Supplementary Information The online version contains supplementary material available at 10.1186/s12862-021-01955-x.

Background
Across angiosperms, flower color is determined by a number of different factors including the presence of one or more pigment classes (e.g., anthocyanins, non-anthocyanin flavonoids, carotenoids including xanthophylls, betalains, and chlorophyll) and environmental factors including pH and UV exposure [93,95]. Among these different classes, anthocyanins are widely appreciated as one of the most important and evolutionarily widespread plant pigment pathways [26,52,72,81]. Anthocyanins are water-soluble pigments that belong to the flavonoid biosynthetic pathway. The three primary aglycone forms of anthocyanin plus three derivatives-pelargonidin, cyanidin (including peonidin), and delphinidin (including petunidin and malvidin)-are together responsible for red, pink, and purple to blue pigmentation, respectively [92].
The Anthocyanin Biosynthetic Pathway (ABP) has been extensively studied from a molecular perspective and is highly conserved across flowering plants. Studies in model species such as petunia, snapdragon, and Page 2 of 17 Zhuang and Manzitto-Tripp BMC Ecology and Evolution (2022) 22:27 Arabidopsis have identified eight structural genes encoding the following core ABP enzymes: chalcone synthase (CHS), chalcone isomerase (CHI), flavonoid 3-hydroxylase (F3H), flavonoid 3'-hydroxylase (F3'H), flavonoid 3′5'-hydroxylase (F3′5'H), dihydroflavonol 4-reductase (DFR), anthocyanidin synthase (ANS), and UDP-3-Oglucosyltransferases (UF3GT). The temporal and spatial expression of these eight genes as well as any structural variations due to mutations in these genes, together with differential binding of pathway intermediates to enzymes depending on substrate specificity, yields flower color variation [26,72]. Anthocyanin biosynthesis is regulated primarily at the transcriptional level by a complex of three classes of transcription factors-MYB, basic helix-loop-helix (bHLH), and WD-40 repeat (WDR) proteins-collectively referred to as the MBW complex [50]. Compared to bHLH and WDR, MYBs play a central role in recognizing specific target genes [6,45,94]. A plethora of studies have identified MYB homologs and their functions in various plant families. An R2R3-MYB encoding gene named VvmybA1 is a key regulator of UFGT in grape skin and other parts of the grape plant [88]. In Antirrhinum majus, expression of F3H, DFR, ANS and UF3GT in flowers is regulated by three closely related R2R3-MYBs named Rosea1, Rosea2 and Venosa [23,62]. Three MYBs isolated from apple, MdMYB1, MdMYB10 and MdMYB, and 'MYB gene Ruby' from Citrus, facilitate the accumulation of anthocyanins in the fruit [14,20]. MYBs furthermore play tissue-and location-specific roles in flower coloration. In Mimulus lewisii, two R2R3-MYBs named PELAN and NEGAN regulate anthocyanin biosynthesis in petal lobes and nectar guide spots, respectively [90].
Significant progress made in elucidating ABP regulatory mechanisms in non-model plants has been possible owing to the conserved nature of ABP enzymes and their MBW regulators. However, study of ABP regulation through homolog-based approaches limits potential identification of novel, regulatory elements because most identified MYBs are highly conserved and belong to the same clade in phylogenetic analyses. Somatic mutations of identified MYBs generally affect the expression of one or a few key structural genes [17], suggesting the existence of other unidentified regulators. In the present study, in addition to identification of genetic elements involved in ABP expression through a homolog-based search method, we undertake co-expression network association analyses (Additional file 1: Fig. S1) using 16 newly generated transcriptomes from 10 closely related species of Ruellia (Wild Petunias: Acanthaceae) to (1) yield a first, comprehensive estimate of candidate structural and regulatory loci underlying flower coloration in Ruellia, (2) examine molecular evolution of these loci, (3) assess whether mutations to structural vs. regulatory elements have likely been of greater importance during flower color evolution in the group, (4) identify blocks of the ABP likely under co-regulation, and (5) identify potentially novel candidate transcription factors involved in the ABP.

Plant materials
This study sampled flower petal (P) and/or leaf (L) tissues from the following species, all of which are in cultivation in glasshouses at the University of Colorado-Boulder: Ruellia bourgaei (P), R. breedlovei (P), R. brevifolia (P), R. elegans (P, L), R. fulgida (P, L), R. hirsuto-glandulosa (P, L), R. longepetiolata (P), R. lutea (P, L), R. simplex (P, L), and R. speciosa (P, L). Flower petal tissues were removed from the S3 stage ( Figure 3) and leaf tissues were collected at a mature stage, approximately 2-3 nodes below the uppermost meristem. All species except R. longepetiolata and R. elegans were derived from wild field populations; the former two species were acquired and grown from cuttings sent by colleagues.

Anthocyanidin identification
We followed anthocyanidin identification and quantification methods described in Zhuang and Tripp [97]. Briefly, 25 mg of silica-dried leaf or petal tissue was placed in 1.5 ml tubes filled with 2 N HCl. Sugars were cleaved from anthocyanin molecules by placing tubes in a 103 °C heat block for 90 min then centrifuged for 5 min at 10,000 rpm at room temperature to obtain supernatants. Retained supernatants were washed twice with 400 uL of ethyl acetate and centrifuged for 1 min at 10,000 rpm to restore phase separation. The pigmented bottom layer was further washed twice with 200 uL of isoamyl alcohol and injected into an Agilent 1260 Infinity system (Thermo Scientific) for anthocyanidin identification and quantification.

cDNA library construction and sequencing
Fresh leaf and petal tissues were collected from plants growing in the glasshouses and placed immediately into liquid N 2 for RNA extraction. Total RNA was extracted with MasterPure ™ RNA Purification Kits (Epicentre) following the manufacturer's protocols. Integrity of total RNA was determined by running samples on a 1.2% agarose gel. RNA-seq library preparation was conducted using an Illumina ScriptSeq Complete Kit, Low Input (BL1224, Illumina) according to the manufacturer's instructions. The sequence-ready libraries were quantified using a Qubit (Invitrogen) and quality was assessed using a Bioanalyzer. Libraries were sent to the Genomics and Microarray Core, University of

Hybrid RAD loci-transcriptome approach to reconstructing a species phylogeny
To understand evolutionary relationships among the ten species used in our study, we developed a simple pipeline that uses unassembled transcriptome reads for phylogenetic reconstruction. First, all error-corrected reads were screened for the presence of the sequence tag "GAA TTC ", which represents the EcoRI enzyme recognition sequence. The sequence tag 5'-AATTC as well as its downstream sequences were extracted. The resulting reads were used for de novo locus assembly implemented in pyRAD v3.0.6 (Eaton [19]-mindepth 15,-clustering threshold 0.85, -minsamples: 3). Resulting biallelic loci were extracted from the final alignment using an inhouse Perl script. We used maximum likelihood methods to reconstruct phylogenetic relationships in RAxML v8 (Stamatakis [68],-f a; -m GTRGAMMA; -N 200).

PCA analysis
To generate a uniform read count table for downstream analysis, identified ORF sequence files from all samples were pooled together and clustered using CD-HIT-EST with a similarity threshold of 95%, and each resulting transcript was considered as a unigene then used as a reference for read mapping. Reads from each tissue-specific library were mapped to the constructed reference using bowtie2 (-sensitive), and eXpress v1.5.1 [58] was used to extract raw read counts and FPKM values for each sample. Transcripts were included in downstream analysis when > 10 reads were recovered in a minimum of five samples. Transcripts failing to meet this threshold were removed from downstream analyses. We then conducted a regularized logarithm (rlog) transformation using the rlog function in DESeq2 [39] on the remaining data. A PCA analysis was conducted on these transformed data using the plotPCA function in DESeq2 and results were visualized using the R package ggplot2 [85].

Identification of ABP structural genes
To identify candidate structural genes involved in the ABP in all species of Ruellia included in this study, six core pathway genes from Arabidopsis (CHS: AT5G13930, CHI: AT3G55120, F3H: AT3G51240, DFR: AT5G42800, ANS: AT4G22880, and UF3GT: AT5G54060) and two branching genes from Penstemon neomexicanus (F3'H: KM388824 and F3′5'H: KM388828) were used as references for homolog identification using NCBI-BLAST [2], with a blast e-value threshold of 1 × 10 −100 for both blastn and blastp, in all instances. To further improve the quality of this dataset, all identified transcripts were further blasted against NCBI's nr Database; candidate transcripts were retained only if the top blast hit was annotated as a structural gene belonging to the ABP. cDNA sequences of transcripts obtained from all species were aligned to their corresponding reference genes to assess sequence integrity using MAFFT v7.3.10 (Katoh et al. [30]; mafft-maxiterate 1000-localpair).
For species for which we were unable to assemble fulllength cDNAs for all eight ABP loci, a second round of de novo assembly was conducted. First, all error-corrected reads were mapped to the reference gene of its most closely related species (see Fig. 1A) using bowtie2 v2.2.6 (Langmead and Salzberg [33],-sensitive). All mapped reads were then extracted using samtools v1.4.1 [35] and served as input for de novo assembly in SPAdes v3.9.0 [5] with default parameters.

Co-expression analysis of candidate ABP genes
To evaluate the expression correlations among identified candidate ABP genes, their FPKM values were used for testing the significance of correlations (an FDR adjusted p-value threshold of < 0.05) using the corr.test() function in the R package "psych" [56] and p-values were adjusted using the Holm-Bonferroni method. All statistical analyses were conducted within the R environment. The correlation was visualized as a matrix using the R package "corrplot" [80]. For candidate ABP structural genes associated with MYB10L1, correlations were predicted using a linear model with a 95% confidence level interval (geom_smooth(method = "lm", level = 0.95)) and visualized using ggplot2 [85].

Non-synonymous to synonymous mutation rates and tests for selection
To test for evidence of positive selection on ABP genes, MAFFT-aligned cDNA sequences and their corresponding coding sequences were used to guide coding sequence (CDS) alignments in Pal2Nal (Suyama et al. [69], http:// Fig. 1 Phylogenetic tree depicting evolutionary relationships among 10 sampled species of Ruellia and principle component analysis (PCA) of transcriptome data from these species. A. Maximum likelihood tree was inferred using a hybrid transcriptome-RADloci approach newly implemented in this study. Briefly, raw, unassembled transcriptome reads were mined for the subset of reads that contained the GAA TTC EcoR1 overhang; these yielded a total 110,263 biallelic sites following data processing in pyRAD. The ML tree was inferred using a GTRGAMMA model of sequence evolution implemented. Shown to the right of taxon labels are: photographs of flowers of the 10 study species as well as a shorthand coding of the major anthocyanin pigments present in petal tissues of these species (C = cyanidin + derivatives; D = delphinin + derivatives; P = pelargonidins+ derivatives). B. Gene expression profiles of ten species projected onto the first two principal components, showing two clear clusters of leaf (triangle) vs. corolla (circle) transcripts as well as two subclusters of corolla transcripts separated by PC2 that correspond to Clade I vs. Clade II in Fig. 1A Page 5 of 17 Zhuang and Manzitto-Tripp BMC Ecology and Evolution (2022) 22:27 www. bork. embl. de/ pal2n al/). The ratio of non-synonymous to synonymous mutations (dN/dS, or ω) was estimated using codon-based maximum likelihood methods implemented in codeml, a module in the PAML package (version 4.8; Yang [89]). To test for positive selection across an entire gene, we used the M0 model, which assumes the same ω ratio for all sites in the gene. To test for positive selection in a site-specific manner (for nuances and best practices of site-specific dN/dS analyses, see, for example, [8,67,78), we first identified sites under positive selection through Naive Empirical Bayes (NEB) analysis (P > 95%) using the site models M1, M2, M7 and M8, which allow the ω ratio to vary among sites. Models M1 + M2 and M7 + M8 form two pairs of models that can be used to test for positive selection using a likelihood ratio test (LRT). Model M1and M2 allow the site classes ω = 1 and 0 < ω < 1, and model M2 additionally includes a third class: ω > 1 [89]. Models M7 and M8 include several site classes with ω ratios that follow a beta distribution, and model M8 additionally includes a ω > 1 class. Two sets of model comparisons were conducted for each gene: the first compares M1/M2 whereas the second compares M7/M8. Codon sites that are under positive selection were identified using the Bayes Empirical Bayes (BEB) method [89].

Secondary structure analysis and functional predictions
To examine potential effects of amino substitutions on protein function, we predicted the secondary structure of proteins as well as functional effects of mutations at amino acid sites under positive selection. Secondary structure of proteins (helices, sheets, and coils) was predicted using the PredictProtein online server (Yachdav et al. [87], https:// www. predi ctpro tein. org, last accessed 07/21/2017). Functional effects of sequence variants were predicted using snap2, which is integrated into the pre-dictProtein server [87].

Identification of ABP interactive pathways through co-expression network analysis
To test for correlations in expression among assembled transcripts and ABP genes, transcripts with at least 10 read counts and transcripts that were found in at least three samples were selected. Pearson correlation coefficients were calculated using the corr.test function in the R package "physic" [56]. To calculate adjusted P values, expression levels of each ABP gene were permuted 1000 times and calculated as (r + 1)/(n + 1), where n = 1000 and r is the number of these replicates that produce a test statistic greater than or equal to that calculated for the actual data [47]. We deemed correlations to be significant when Pearson's r ≥ 0.65 and P < 0.05. All transcripts with expression levels significantly associated with the 12 identified candidate structural genes were annotated using Trinotate pipeline release version 3.0.0 (Grabherr et al. [24], http:// trino tate. sourc eforge. net/) with a blast e-value threshold of 1 × 10 −5 . We additionally built a custom plant transcription factor database for ABP loci using data downloaded from the Plant Transcription Factor Database v4.0 (Pérez-Rodríguez et al. [51]; http:// plant tfdb. cbi. pku. edu. cn, last accessed 04/18/2017). KEGG pathway analysis was conducted on all resulting transcripts using BlastKOALA (Kanehisa et al. [29], http:// www. kegg. jp/ blast koala/ Last accessed 4/01/2017) and was further processed using an in-house Perl script then visualized using the R package ggplot2 [85].

Phylogenetic analysis and sequence comparison
We aligned translated protein sequences of ABP genes using MAFFT v7.3.10 (Katoh et al. [30]; mafft-maxiterate 1000-localpair) and then identified conserved regions of resulting alignments using Gblocks [70]. Phylogenetic analyses were conducted on final alignments using RAxML and the GTRGAMMA substitution model, which is more computationally intensive compared to GTR CAT but is suitable for datasets with smaller numbers of taxa, such as ours (see RAxML documentation: https:// cme.h-its. org/ exeli xis/ resou rce/ downl oad/ NewMa nual. pdf ) with 100 bootstrap replicates.

Real time qPCR
To confirm s transcript's abundance measured from the RNA-Seq data and better characterize assembled homologs of F3'H and DFR, we extracted total RNA from petal tissues of Ruellia simplex collected at different growth stages using MasterPure ™ RNA Purification Kit (Epicentre) following the manufacturer's protocols. For reverse transcription, ~ 500 ng of DNase I-treated total RNA was used for cDNA synthesis using the GoScript ™ Reverse Transcriptase system (Promega) with the oligo(dT) 15

HPLC analysis of anthocyanin accumulation in leaf and petal tissues
To phenotype ABP pigments that underlie flower color variation in Ruellia, 10 species with distinct petal colors (four purple, three red, and three yellow) were selected based on visual inspection coupled with floral color delimitation into one of several discrete categories. The latter was accomplished via floral reflectance measured using an Ocean Optics JAZ-COMBO Spectrometer ( Fig. 1A; Additional file 15: Data S1

De novo assembly of tissue-specific Ruellia transcriptomes
To help elucidate molecular mechanisms regulating anthocyanin accumulation, we analyzed data from 12 newly sequenced transcriptome libraries (eight from flower petals, four from leaf tissues) together with data from four previously published libraries (two from flower petals, two from leaf tissues; [97]. Quality statistics for tissue-specific assemblies using the best TransRate score are shown in Table 1 nearly all (98.78 ± 1.04%) could be mapped concordantly as read pairs. We also evaluated transcriptome completeness by searching for single-copy orthologs of highly conserved plant genes in each assembly using BUSCO v2 [64]. On average, assembled transcriptomes, regardless of whether they derived from petal or leaf tissue, contained a high percentage (74.06 ± 5.91%) of the 1440 groups of conserved plant genes deposited in the BUSCO database; this high percentage further reflects the completeness of our transcriptome assemblies.

Reconstruction of phylogenetic relationships in Ruellia
To resolve phylogenetic relationships among the 10 species of Ruellia, we employed a new approach that involves retrieving RAD loci from transcriptome data. Phylogenies constructed using assembled transcriptome data sometimes face challenges due to complications associated with pseudogenes and paralogous loci, resulting in potentially misassembled contigs [82]. We thus implemented a hybrid transcriptome-RADseq phylogenetic tree building approach that made use of a subset of error-corrected reads having the EcoR1 GAA TTC overhang. These reads were treated in a manner similar to RADtags and were analyzed using pyRAD [19]. The minimum locus length was set to 35, but the resulting average read length was 78nt, which is suitable for phylogenetic analysis as it has been shown that large phylogenetic tree can be built accurately from even short (50 bp average) sequences [27]. A total of 110,263 loci were concatenated and used for phylogenetic inference via maximum likelihood implemented in RAxML v8 [68]. Branch support as assessed via 100 bootstrap replicates was ≥ 75% for all nodes (most nodes with ≥ 95% support), and the resulting phylogeny composed of two major clades (Fig. 1A) mirrors relationships recovered in prior phylogenetic studies in the genus [75]. We specifically did not conduct ancestral state reconstructions of flower color and estimations of color transitions due to very few species terminals [59], especially in a taxonomically species-rich genus such as Ruellia (ca. 350 species worldwide). However, the 10-tip phylogeny suggests four potential color transitions: in Clade I, one from purple to red and one from red to yellow flowers, and in Clade II, one from red to purple and one from purple to yellow flowers.

PCA analysis of leaf and petal transcriptomes
To assess overall expression level similarities among transcriptomes, we conducted a principal component analysis (PCA) on normalized gene count data based on combined petal and leaf transcriptome data from R. simplex, which had the best assembly statistics among all species (Table 1). After removing transcripts with low expression levels (see "Methods"), a total of 60,626 transcripts were used in PCA analysis. Results demonstrate that flower petal vs. leaf transcripts were separated into two clusters (Fig. 1B). Additionally, two subclusters among the flower petal transcripts can be identified and correspond to the two clades resolved in our phylogeny (Fig. 1A

Identification and co-expression analysis of candidate ABP structural genes
Our BLAST + search of Trinity-assembled Ruellia transcriptomes against known ABP genes identified 12 orthologous candidate ABP structural genes in Ruellia. Among these, one copy each of CHS, CHI, F3H, F3′5'H, ANS and UF3GT was identified whereas three copies each of F3'H and DFR were recovered. Subsequent reference-based gene assembly using SPAdes v3.10.1 [5] recovered full-length cDNA of ABP genes in most of the species (Additional file 12: Table S1, Additional file 2: Fig.  S2). We then measured expression levels of the 12 candidate genes in FPKM (fragments per kilobase of transcript per million fragments sequenced) and compared results in a tissue-specific manner. Results indicate that expression of DFR1 and F3′5'H was restricted to petal tissues whereas expression of all other genes occurred in one or more leaf samples (Fig. 2). We further characterized the multiple copies of F3'H and DFR identified among our Ruellia transcriptomes by comparing translated amino acid sequence across samples (Additional file 2: Fig. S2) as well as expression patterns measured at five flower developmental stages in Ruellia simplex using transcript-specific primers and qPCR (Fig. 3, Additional files 13, 14: Table S2, S3). Two of the three F3'H copies (F3'H1, F3'H2) were highly homologous (sequence divergence: 23.44% across all samples; Additional file 2: Fig. S2) whereas F3'H3 was markedly divergent from the former two (Additional file 2: Fig.  S2). Similarly, DFR3 is highly divergent from the other two copies (i.e., DFR1 and DFR2). This marked divergence in F3'H3 and DFR3 coupled with more widespread expression of these two copies in leaf tissues across species (compared to expression of F3'H1, DFR1, DFR2, and to a lesser extent F3'H2 in petal tissues; Fig. 2) indicates potential functional diversity of these two enzymes.
Using qRT-PCR, developmental stage-specific expression of multiple copies of F3'H and DFR found in R. simplex petals (i.e., DFR1, DFR2, DFR3 and F3'H2 and F3'H3 [F3'H1 not expressed in petals of this species; Fig. 2]) confirmed transcriptome-based expression patterns of ABP enzymes in this species (Fig. 3). These data further support the conclusion that our transcriptome assemblies are relatively complete and are of high quality. Analyses of R. simplex indicated that expression of F3'H2, F3'H3 and DFR3 all peaked at stage two of five developmental stages, also indicating potentially coordinated regulation of these three genes (Fig. 3).
Network based co-expression analysis across tissue types in all samples recovered strong, correlated expression among CHS, F3H, DFR1, and ANS as well as between CHI and DFR3 (r ≥ 0.65), indicating possible blocks of structural genes under regulatory control  Fig. 4). Taken together, results indicate that regulation of ABP is highly coordinated and likely involves different regulatory blocks.

Identification and co-expression analysis of candidate ABP regulatory genes
Our homology-based searches for candidate MBW regulatory elements in Ruellia recovered four clades of MYBs (Fig. 5), two of WDRs (Additional file 3: Fig. S3), and three of bHLHs (Additional file 4: Fig. S4). Transcript expression indicated that the three families of transcription factors exhibited different expression patterns: expression of all four clades of MYBs was specific to petals (Fig. 5) whereas WDRs and bHLHs were expressed in both petal and leaf tissues (Additional files 3, 4: Figs. S3 and S4). Of the four MYB clades, copies of MYB10L1 were relatively highly expressed compared to copies in the MYB10L2, MYB10L3, and MYB10L4 clades (Fig. 5B).
MYB10L3 and MYB10L4 were expressed only in a few species and at relatively low levels (Fig. 5). Network based co-expression analysis across tissue types in all samples recovered strong, correlated expression of MYB10L1 with five candidate structural genes (CHS, F3H, DFR1, DFR2 and ANS; r ≥ 0.65; Fig. 4) and a weak correlation with F3′5'H (r = 0.56). No correlations among expression of MYB10L2, WDRs, bHLHs, and candidate structural genes were found; co-expression analyses were not conducted for MYB10L3 and MYB10L4 due to limited data.
We compared amino acid-predicted sequences of MYBs from our four identified clades in Ruellia to amino acid translations of functionally validated MYBs known to regulate ABP in other flowering plants: Arabidopsis thaliana, Petunia x hybrida, Pyrus pyrifolia, and Mimulus guttatus. Specific attention was given to the amino acid sequence residue ' ANDV' because this residue has been reported to distinguish anthocyanin from non-anthocyanin MYBs [37]. We found this residue to be present in candidate Ruellia MYBs from the MYB10L1 and MYB10L4 clades ( Fig. 5B; note that expression levels of MYB10L4 copies were markedly lower than were expression levels of MYB10L1 copies). In contrast, replacement of the amino acid V with an I (resulting in the residue ' ANDI' and considered to be a conserved signature for non-anthocyanin gene regulation) characterized most of the MYB10L2s and MYB10L3s (as well as two of the four Mimulus guttatus sequences; Fig. 5B). We detected an additional non-synonymous mutation in this residue in Ruellia lutea, which resulted in an A to an S amino acid substitution (yielding 'SNDV' ; Fig. 5B).

Estimating amino acid substitutions in ABP loci and potential effects on protein function
To test whether codon sites of identified ABP loci have potentially undergone adaptive evolution, we conducted tests of non-synonymous to synonymous substitution rates (dN/dS ratios, or ω) across all candidate structural genes plus one candidate MYB predicted to be significantly associated with the ABP in Ruellia (see Fig. 4). We used genomic data from the closely related species Ruellia matudae to serve as an outgroup for dN/dS analyses. First, we implemented the site test model M0 in codeml, which assumes a single ω ratio across all sites in an alignment. Results yielded dN/dS ratios that ranged between 0.20 and 0.72 (Table 2). Among all loci tested, DFR3 had the highest dN/dS ratio (0.72) followed by F3'H2 (0.54) and MYB10L (0.52). The fact that no gene was found to have dN/dS ratio > 1 indicates that all ABP genes are predominantly under purifying selection. However, because positive selection is unlikely to affect all sites in a given protein over extended evolutionary time, we further explored the potential of positive selection on individual amino acid sites. Results indicated that several sites in both DFR3 and F3'H2 are likely to be under positive selection (Fig. 6). Subsequent protein secondary structure analysis indicated that four of the sites that appear to be under positive selection in F3'H2 and five sites in DFR3 are located inside protein binding regions (PBRs). We further examined the potential effects of amino acid substitutions on protein function at the above sites likely to be under positive selection. Using algorithms that predict both secondary structure and functional effects, we found one substitution in F3'H2 and three substitutions in DFR3 for which a moderate functional effect of the mutation is predicted (colored as pink in Fig. 6). These results indicate that, at these sites, observed amino acid substitutions are more likely to influence the stability or substrate binding affinity of the resulting protein or even cause protein dysfunction. At all other sites, observed amino acid substitutions colored as green indicate a strong signal for neutral effect.

Co-expression analysis of transcripts and associations with ABP loci
After removing transcripts that either had low expression levels or were expressed only in a small fraction of our study species (i.e., fewer than 3), a total of 60,613 transcripts were tested for correlated expression with the 13 ABP structural genes identified in Ruellia. Significant positive correlations were found for 1748 transcripts (2.89% of total transcripts present in 3 or more species, Pearson's r ≥ 0.65, P < 0.05). We used the KEGG classification system to predict and annotate functionality of these ABP-associated transcripts by assigning them to known biochemical pathways. Of the 1748 transcripts, 532 (30.43%) were assigned KEGG annotations. Among these, genes involved in a variety of metabolic pathways formed the largest group (Additional files 5, 6: Figs. S5 and S6), among which 10 are putatively involved in phenylpropanoid biosynthesis and two are involved in flavone and flavone biosynthesis.
To gain further insight into interactions between transcripts and ABP loci as well as identify potential novel regulatory elements, transcription factors were further characterized using a custom database derived from the Plant Transcription Factor Database. We four MYB orthologs from Mimulus guttatus were also included in phylogenetic analyses given the close evolutionary relationship of Mimulus to Ruellia, but these four copies have not been functionally validated in prior works). Four, strongly supported clades of candidate MYBs in Ruellia were recovered (highlighted in yellow). B. Protein sequence alignment of candidate ABP MYB regulators in Ruellia and orthologs in other species. Relative expression of each gene shown in FPKM (fragments per kilobase of transcript per million fragments sequenced). Results indicate that copies from the MYB10L1 clade in Ruellia were among the most highly expressed and also contain an amino acid residue (' ANDV') that has been reported to characterize MYBs that regulate the ABP (see text); this residue was also found copies from the MYB10L4 clade, but transcripts in this clade were in low abundance. In contrast, a non-synonymous substitution yielding an ' ANDI' residue characterized most of the candidate Ruellia MYBs in the MYB10L2 and MYB10L3 clades, and this residue is conserved across MYBs that function elsewhere besides in ABP production in other plant families identified 13 families of TFs that are positively correlated to ABP loci expression in Ruellia (Additional file 5: Fig. S5). Among these, MYBs were the dominant group of ABP-associated TFs. Among 19 specific MYBs identified, nine belonged to the R2R3 family including MYB10L1 that we earlier identified as likely to be especially relevant to ABP expression in Ruellia based on homology to functionally verified MYBs in other plant families (Fig. 5). To further characterize potential function roles of the remaining eight R2R3 MYBs, we used   [3,15,21,26,32,46,50,79]).

Discussion
In the present work, we have documented the effects of mutations in ABP genes-coding and regulatory-and the effects of expression of these genes on flower color variation. Our results demonstrate the importance of both structural and regulatory changes during ABP evolution in Ruellia. We identified 12 candidate ABP structural genes plus 9 candidate transcriptional regulators, including 4 MYBs, 2 WDRs, and 3 bHLHs among the 10 investigated species. Through co-expression network analyses, we found evidence for strong, correlated expression among most ABP structural genes as well as between most of these genes and one candidate regulatory MYB.
Phylogenetic relationships among the 10 species here investigated yield a hypothesis of four potential flower color transitions within the group: purple to red, red to yellow, red to purple, and purple to yellow (Fig. 1A). These four flower color transitions are consistent with patterns documented in a study with much more extensive phylogenetic sampling [74] and implicate both gain and loss of function mutations during flower color evolution in Ruellia. Prior studies on ABP have yielded evidence that four classes of mutations contribute importantly to color shifts across angiosperms: (1) mutations in ABP structural genes that cause either loss of function or changes to enzyme binding affinity, (2) mutations in cis-regulatory regions of ABP structural genes, and (3) mutations in coding regions or (4) cisregions of ABP regulatory genes that impact expression of structural genes [17,26,37,43,49,61,63,84,88,90,92]. The purple to red color transition that occurred between the R. breedlovei lineage and the clade containing R. fulgida likely involved a loss of function mutation in the F3′5'H gene in the latter: our data showed that all essential ABP structural genes including F3′5'H itself were expressed at detectable levels (Fig. 2), but sequence analysis of F3′5'H revealed a premature stop codon caused by a C to T point mutation at positon 373 in R. fulgida (Additional file 8: Fig. S8). By contrast, both our de novo and reference-based approaches failed to assemble F3′5'H in the red-flowered species R. brevifolia, which belongs to the same clade as R. fulgida, and expression analysis indicated F3′5'H was expressed at extremely low levels (Fig. 2). These data indicate that a regulatory mutation was likely involved in this purple to red color transition. In the yellowflowered species R. bourgaei, R. speciosa, and R. lutea, we observed down-regulation of multiple structural genes, including down-regulation of F3H in the closely related R. bourgaei and R. speciosa as well as down-regulation of ANS in all three species (Fig. 2). These results suggest either independent loss of functions in the cisregulatory elements of these structural genes in the two clades or, more likely, mutations in shared regulatory genes given that expression of F3H and ANS is strongly coordinately regulated (Fig. 4). Finally, the shift from red to purple flowers that occurred between Ruellia elegans and R. hirsuto-glandulosa likely involved reactivation of the F3′5'H pathway branch. Although restoration of pathway function is expected to be more difficult than degradation of function [54,66], putative gains in floral anthocyanins have been documented ( [4], 31). In R. elegans, malvidin-a derivative of delphinidin-was detected in leaf tissues albeit in low concentrations, indicating functionality of the F3′5'H pathway in this species (Additional file 16: Data S2). Restoration of the F3′5'H pathway in petals of R. hirsuto-glandulosa may have been achieved by mutation in the cisregulatory region of this enzyme, leading to decreased binding affinity of transcription inhibitors or increased binding affinity of transcription activators.
With exception of a few mutations in F3'H2 and DFR3 (Fig. 6) and two mutations yielding premature stop codons in F3′5'H (Additional file 8: Fig. S8), sequence comparisons and protein predictions of ABP structural genes among species with different flower colors indicated that the majority of amino acid mutations likely had neutral or minimally different functional effects, thus unlikely to result in enzyme dysfunction (Additional file 17: Data S3). Thus, much of the observed variation in flower colors and anthocyanin accumulation in Ruellia is more likely to be the result of mutations impacting the expression of ABP structural genes (Fig. 2). Regulatory mutations are considered to be among the most important factors in driving morphological evolution [7,84] given that such changes can potentially lead to modification in expression of entire pathways. Through coexpression analysis, we found that expression of one of the R2R3 MYBs identified in our study-MYB10L1, which is homologous to the known ABP regulator MdMYB10 in apples (Malus × domestica; 79)-was highly associated with several structural genes including CHS, F3H, DFR1, DFR2, and ANS (Fig. 4). Thus, mutations that affect MYB10L1, be they regulatory or coding, may potentially impact all five of these candidate structural genes. However, sequence comparison of MYB10L1 across all species of Ruellia sampled here recovered no lethal mutations (Fig. 5, Additional file 17: Data S3), indicating MYB10L1 itself is likely regulated by a higher order of regulators.
Notably, based on patterns of expression correlation with MYB10L1, the regulation of candidate ABP structural genes in our dataset can be separated into two distinct regulatory blocks (Fig. 4). Group 1 contains genes whose expression are likely under the control of MYB10L1: CHS, F3H, F3′5'H (significantly correlated with DFR2, r = 0.65, weakly associated with MYB10L1, r = 0.56), DFR1, DFR2, and ANS. Group 2 contains genes whose expression is independent of MYB10L1, among which CHI and DFR3 seem to be regulated in a coordinated manner (Fig. 4). In contrast to these two blocks, each copy of F3'H and the single copy UF3GT are regulated independently. While independent regulation of F3'H has been observed in other studies, co-expression of UF3GT with CHS or F3H in Arabidopsis has been reported in the past [1,63]. Thus, the observed expression pattern of UF3GT in the present study seems to be a unique feature of our dataset.
Prior research has demonstrated lower rates of nonsynonymous mutations among genes upstream in the ABP pathway compared to downstream genes [53], likely as a function of relaxed evolutionary constraint on downstream (vs. upstream) genes [40] or higher rates of positive selection on downstream genes, potentially playing a role in adaptive floral color evolution [28]. This pattern is expected given that upstream enzymes function in a much broader set of biosynthetic pathways (e.g., lignin and other flavonoids) than downstream enzymes, which are largely specific to anthocyanin production (Additional file 9: Fig. S9). Data from Ruellia reveal similar patterns of rate evolution of ABP genes compared to prior studies. Specifically, two of the lowest dn/ds substitution rates were documented for CHS and CHI (Table 2), consistent with a hypothesis that these two genes are undergoing stronger purifying selection. Furthermore, across all candidate ABP loci, loss of expression was observed only in instances of genes downstream from CHS and CHI (Fig. 2), likely reflecting the essential role of these two enzymes in other metabolic pathways besides the ABP [86].
Flower color often reflects adaptation to pollinators [11]. For example, floral anthocyanin content has a direct effect on pollinator behavior [60] and plant fitness [11]. That all candidate structural genes identified in the present study are under purifying selection adds to evidence that selection acts to maintain overall primary protein function of ABP genes. However, our data also yielded evidence that selection may drive diversification of ABP products in some lineages, specifically, we found signatures of positive selection in several sites in DFR and F3H across species (Fig. 6), and some of these sites are located in or nearby protein binding regions (PBR). This is especially the case for DFR3, in which five out of seven sites under positive selection are located within predicted PBR Mutations occurring at these sites are more likely to impact substrate binding affinities of the corresponding enzyme resulting in new interactions [55]. In this study, we found multiple, orthologous copies of both F3H and DFR to be present and here hypothesize that duplication of these genes may have facilitated functional diversification of the ABP in Ruellia [28].
Transcription factors (TFs) are key players in regulating flux through secondary metabolic pathways by controlling relative levels of gene expression [13]. To further investigate regulation of the ABP in the context of potential interactions with other metabolic pathways, we conducted a phylogenetic study of candidate MYBs as well as co-expression analyses of broader diversity of TFs across all species. We identified nine candidate R2R3 MYBs associated with ABP genes in Ruellia. Among these, orthologs of RsMYB10L1, RsMYB114L (i.e., AtMYB75, [10]), RsMYB305L1, and RsMYB305L2 (i.e., NtMYB305; Wang et al. [79]) have been reported to function in regulating anthocyanin accumulation. Our investigation demonstrated that RsMYB10L1 expression is highly correlated with five ABP structural genes (Fig. 4), and phylogenetic analyses demonstrated that RsMYB10L1 belongs to a strongly supported clade that contains AtMYB75 orthologs (Additional file 7: Fig. S7). This result in combination with sequence analysis and other expression patterns (Fig. 5) indicate specifically that RsMYB10L1 is a very strong candidate regulator of the ABP in Ruellia. Pilot experimentation with virus-induced gene silencing focusing on ABP knockdowns of RsMYB10L1 yielded additional support for this candidate (Zhuang and Manzitto-Tripp, unpub. data; Additional file 10: Fig. S10). Our phylogenetic analyses placed the remaining Ruellia R2R3 MYBs in a clade with Arabidopsis orthologs that function in lignin biosynthesis (e.g., AtMYB4, AtMYB31, and AtMYB46/AtMYB83, Additional file 7: Fig. S7; [3,15,32,46,50]. These results support a hypothesis of different functional roles of the R2R3 MYBs recovered in Ruellia-some in ABP production vs. others potentially in lignin production-albeit in pathways that interact via common upstream enzymes (Additional file 9: Fig. S9, Tohge et al. [73]). Functionally, anthocyanins and lignins are both proposed to protect plants against UV radiation [44], and redirection of metabolic fluxes between the lignin and flavonoid pathways has been observed in many plant species [57]. Over-expression of AmMYB308 in transgenic tobacco plants and several of its orthologs belonging to the same clade led to significant repression of lignin biosynthesis and increased total flavonoids [42,71]. Other studies have similarly demonstrated MYBs as negative regulators of lignin content [3] whereas others likely serve as activators of lignin biosynthesis [12,32,96]. Taken together, data suggest extensive interactions between anthocyanin and lignin biosynthesis through modification of expression of relevant TFs. As such, future attempts to understand regulatory impacts on flower color evolution in Ruellia would optimally incorporate information from other biosynthetic pathways that likely interact with anthocyanin production (e.g., [76]).
We also identified several TFs that are highly coexpressed with ABP-associated TFs in Ruellia. These include auxin responsive factors (ARFs), abscisic acid signaling genes (ABA), and other TFs, some of which function as activators and others as repressors based on prior studies in other plants (Additional file 5: Fig. S5; [18,25,48,77,83]). In the present study, we found positive co-expression between two auxin signaling pathway activators ARF6 and ARF19 [48,77] and ABP-associated TFs in Ruellia. We speculate that this positive association results from the role of auxin in flower development. Specifically, the accumulation of anthocyanins gradually increases during flower development in Ruellia (shown in Fig. 3A). Auxin is required for the initiation of floral primordia and disruption of auxin biosynthesis leads to the failure of flower formation [16]. Our hypothesis is further supported by the observed correlation of a group of MADS-box TFs with ABP structural genes, which is exclusively expressed only in petal tissues (Additional file 11: Fig. S11). Sequence analysis indicated one gene was homologous to AGL6, two genes were homologous to MADS1, and another two genes were homologous to MADS2, all of which impact flower development [41,61]. Taken together, our data reveal strong, coordinated regulation between anthocyanin accumulation and flower development.

Conclusions
The ABP is one of the most well-studied secondary metabolic pathways in plants. Although key structural genes have been extensively characterized and are highly conserved across unrelated species [52], the complex mechanisms regulating the ABP are still not fully understood largely as a result of prior emphasis on model plant systems. The present study lays the first steps towards understanding molecular evolution and expression of ABP genes in Ruellia. Specifically, it yields a first comprehensive survey of candidate structural and regulatory loci involved in the production of diverse flower colors in the genus. Data from our expression and co-expression analyses indicate that regulatory mutations have likely been more important in color shifts but that mutations to structural genes have also played a role. Finally, different portions of the ABP appear to be under coordinated regulation such that a single 'regulatory block' is unlikely to characterize the entire pathway, perhaps as a function of extensive interactions with other metabolic pathways.


s2fieldsofstudy: Biology Medicine

text:
Clinical Implications of Circulating Circular RNAs in Lung Cancer

Circular RNAs (circRNAs) are single-stranded RNAs with a covalently closed-loop structure that increases their stability; thus, they are more advantageous to use as liquid biopsy markers than linear RNAs. circRNAs are thought to be generated by back-splicing of pre-mRNA transcripts, which can be facilitated by reverse complementary sequences in the flanking introns and trans-acting factors, such as splicing regulatory factors and RNA-binding factors. circRNAs function as miRNA sponges, interact with target proteins, regulate the stability and translatability of other mRNAs, regulate gene expression, and produce microproteins. circRNAs are also found in the body fluids of cancer patients, including plasma, saliva, urine, and cerebrospinal fluid, and these “circulating circRNAs” can be used as cancer biomarkers. In lung cancer, some circulating circRNAs have been reported to regulate cancer progression and drug resistance. Circulating circRNAs have significant diagnostic value and are associated with the prognosis of lung cancer patients. Owing to their functional versatility, heightened stability, and practical applicability, circulating circRNAs represent promising biomarkers for lung cancer diagnosis, prognosis, and treatment monitoring.

Introduction
Despite advancements in breakthrough therapies such as targeted therapy and immunotherapy, the survival rate of lung cancer patients has failed to improve for decades. Lung cancer remains the leading cause of cancer-related death worldwide, while the 5-year survival rate of patients with distant metastatic lung cancer is only 6% [1]. Thus, these outcomes urgently require the development of an effective early diagnosis method for lung cancer. Currently, early screening methods for lung cancer in clinical practice include sputum cytology, low-dose chest computed tomography, and autofluorescence bronchoscopy [2]. For the minimal or non-invasive early diagnosis of lung cancer, blood-based biomarkers are rapidly emerging as new alternatives and include circulating tumor cells, exosomes, and circulating nucleic acids (DNAs, microRNAs, and non-coding RNAs) [3,4].
Circular RNAs (circRNAs) are single-stranded endogenous RNAs with a covalently closed-loop structure [5]. Since their discovery in the 1970s, the study of circRNAs has been limited; however, with the development of next-generation sequencing, over the past decade, there has been an increased interest and research in circRNAs. As the physiological and pathological functions of circRNAs have been discovered, many studies have attempted to use circRNAs as biomarkers in the diagnosis of cancer. Owing to their stable structure [5], circulating circRNAs are thought to be more advantageous than normal linear RNAs for application in liquid biopsies. In this review, we briefly discuss the biogenesis and functions of circRNAs and describe in depth the recent findings and prospects for the application of

Biogenesis of circRNAs
The formation of a circular structure via the covalent phosphodiester linkage between the 3 -hydroxyl group of a downstream exon and the 5 -phosphate group of an upstream exon is a seemingly sporadic event for RNA molecules. Although the exact mechanism through which circRNAs are generated has not been elucidated yet, circRNAs are considered to be generated with the help of specialized factors during the splicing process. An alternative "back-splicing" that joins together the 5 GU donor site at the downstream intron and the 3 AG acceptor site at the upstream exon promotes RNA circularization [20]. circRNAs can form from all precursor RNA regions including exons (protein-coding or untranslated regions) and introns.
During the typical mRNA splicing process, a lariat-structured intermediate molecule is produced [21], which can similarly generate circRNAs. The lariat formation catalyzes back-splicing by bringing the splicing sites closer together, thereby potentially providing an optimal situation for circRNA splicing [22]. However, circRNAs might also be produced through an alternative mechanism of direct back-splicing. Whereby, the back-splicing processes a pre-mRNA transcript into a circRNA along with an intermediate transcript composed of exons and introns, which is further processed into a linear RNA [23].
An alternative model for RNA circularization is through complementary sequencemediated RNA paring across flanking introns ( Figure 1A) [15]. Here, the production of normal linear mRNAs or circRNAs from a single precursor transcript could be determined by the competitive RNA paring between competitive sequences within individual introns or across flanking introns. All these models have been supported by evidence produced both in vitro and in vivo; however, it is still elusive how circRNAs are generated, although, importantly, circRNAs can be produced via different mechanisms depending on individual RNA transcripts or the surrounding molecular and cellular context [15]. circRNA biogenesis is precisely controlled by cis-regulatory elements and trans-acting factors, in addition to the spliceosome. Reverse complementary sequences in the flanking introns are particularly crucial in the production of circRNA [15,24]. The reverse complementary sequences are either repetitive sequences such as Alu elements, which are found in the flanking introns of 88% of circRNAs in humans [24], or non-repetitive sequences such as "GUUG" or "ACUU" regions participating in circLONP2 production [25]. The pairing between reverse complementary sequences brings the downstream splice donor and upstream splice acceptor sites closer together, increasing the likelihood of back-splicing events.
RNA circularization can be facilitated by various trans-acting factors, which help bring the circle-forming exons together via the interaction with intronic sequence motifs. Spliceosome complexes (small nuclear ribonucleoproteins) and splicing regulatory factors (heterogeneous nuclear ribonucleoproteins and serine/arginine-rich proteins) regulate the production of circRNAs during the splicing of mRNA precursors [26]. Intronic reverse complementary sequences encourage the initiation of RNA circularization, while the splicing factors exquisitely control the expression levels of specific circRNAs. Indeed, the U2AF2 splicing factor binds to and facilitates the production of circARF, ultimately promoting glioma tumorigenesis [27]. Another splicing factor, ESRP1 expedited the production of circUHRF1 through the interaction with intronic repeat sequences, which then promoted oral squamous cell carcinoma tumorigenesis [28]. In the mouse cortex, global circRNA biogenesis was suppressed by knocking out NOVA2, an RNA-binding protein regulating alternative splicing. NOVA2 was shown to bind to the flanking introns and facilitate backsplicing [29]. QKI is also known to bind to the flanking introns of circRNA-forming RNA precursors and promote circRNA formation during epithelial-to-mesenchymal transition in epithelial cancer cells [30]. Various cis-regulatory elements and trans-acting factors have been revealed to regulate the biogenesis of circRNAs [31]; however, the exact mechanism is still unclear, and more in-depth studies are required.

Functions of circRNAs
Like other regulatory non-coding RNAs, circRNAs interact with genomic DNA, mR-NAs/miRNAs, and proteins to perform a variety of activities within cells ( Figure 1B-F). As previously stated, the role of circRNAs as miRNA sponges has received the greatest attention. CDR1as is the first identified circRNA functioning as an efficient miRNA sponge that contains 74 miR-7-binding seed matches [14,16]. Upregulated in gastric cancer, circ-NRIP1 was shown to work as a sponge against miR-149-5p and promote gastric cancer tumorigenesis and metastasis through the regulation of AKT/mTOR signaling [32]. In non-small cell lung cancer, circPIP5K1A (circ_001430) upregulates BCL2 expression and thus inhibits apoptosis via sponging miR-136-3p [33].
circRNAs also bind to proteins and exert various functions. circTNPO3 binds to and decoys IGFBP3 protein to prevent it from stabilizing MYC mRNA. circTNPO3, therefore, downregulates MYC and SNAIL expression and suppresses proliferation and metastasis in gastric cancer [34]. During neuronal differentiation, circZNF827 was identified as a binding scaffold for a transcription-repressive complex composed of hnRNP-K/L proteins and the ZNF827 protein encoded from its host gene [35]. circPABPC1 inhibited the migration and adhesion of hepatocellular carcinoma cells. Intriguingly, circPABPC1 directly bound to and guided β1-integrin to the 26S proteasome, ultimately accelerating its degradation [36].
In addition, circZNF609 directly interacts with several mRNAs including CKAP5, UPF2, and SRRM1 mRNAs. circZNF609 enhances the stability and translatability of these mRNAs by recruiting the ELAVL1 RNA-binding protein, which controls microtubule dynamics and drug resistance in rhabdomyosarcoma cells [37]. circYAP also specifically recognizes YAP mRNA and then suppresses its translation by inhibiting the assembly of the translation initiation machinery [38]. Furthermore, circSMARCA5 directly binds to the genomic DNA of its host gene, SMARCA5, and suppresses its expression. The SMARCA5 downregulation induced by circSMARCA5 inhibited DNA damage repair and enhanced the drug sensitivity of breast cancer cells [39].
Interestingly, not all circRNAs are non-coding RNAs. In Drosophila, several circRNAs were found to be associated with ribosomes and translated into proteins [40]. circSHPRH (circ_0001649), which was downregulated in glioblastomas, suppressed tumorigenesis via its translational product, SHPRH-146aa [41]. In hepatocellular carcinoma, circβ-catenin promoted tumor growth and metastasis through the Wnt signaling pathway, which was mediated through translation into a novel isoform, β-catenin-370aa. This "microprotein" worked as a decoy for GSK3β, preventing the degradation of β-catenin [42]. Mechanistically, it was demonstrated that internal ribosome entry sites and 18S rRNA complementary sequences facilitate circRNA translation [43]. Some examples of circRNA translation that have been experimentally validated might be artifacts [44]; nevertheless, microproteins translated from specific circRNAs further expand the functional diversity and application range of circRNAs.

Circulating circRNAs
As previously described, circRNAs with diverse functions are involved in cancer development, progression, and metastasis. Like those of mRNAs, miRNAs, and proteins, the expression levels of circRNAs vary depending on the cell type; thus, they can be applied as diagnostic or prognostic markers in cancer patients. Numerous circRNAs have been reported to be upregulated or downregulated in various types of cancer. circRNAs have been discovered in plasma, saliva, urine, and cerebrospinal fluid, so that circulating circRNAs can be used as cancer biomarkers [45,46]. Like other linear RNAs, circRNAs can be amplified through reverse-transcription PCR (RT-PCR), which makes them more easily detectible than protein markers. Unlike linear RNAs, circRNAs lack free 5 -and 3 -ends, making them highly resistant to degradation by RNases with exonuclease activity.
It has been demonstrated that circRNAs are abundant and stable in exosomes, suggesting their significant translational potential as circulating biomarkers for cancer diagnosis [47]. In hepatocellular carcinoma, exosomal circPTGR1 was shown to promote cancer progression through the regulation of the miR-449a/MET pathway [48]. Similarly, circNRIP1 was also proven to be transmitted via exosomes and promoted tumorigenesis and metastasis of gastric cancer [32]. In laryngeal squamous cell carcinoma, circRASSF2 was secreted by exosomes and promoted tumor growth through the regulation of the miR-302b-3p/IGF-1R pathway [49]. High circCNOT2 expression was associated with poor progression-free survival of patients with breast cancer, and circCNOT is detectable in cell-free RNAs from patient plasma samples [50]. In addition, circRNAs can also be detected in circulating tumor cells [51]. Furthermore, circRNAs have been shown to be highly enriched in blood platelets compared with nucleated cells, which can be used for cancer diagnosis [52]. Since protein carriers such as high-density lipoprotein and Argonaute 2 transport miRNAs [53], the circulation of circRNAs might also be mediated by certain protein carriers or RNA-binding proteins.
Since there is a multitude of evidence indicating that circRNA could be useful in liquid biopsy, this review will highlight how circulating circRNAs can be used for the diagnosis and prognosis of lung cancer.

Cancer Progression
The presence of F-circEA generated from the EML-ALK fusion gene was verified in non-small cell lung cancer (NSCLC) cells and in the plasma of NSCLC patients [54]. F-circEA promoted cancer cell migration and invasion, suggesting that F-circEA could be a novel liquid biopsy marker for NSCLC. Through circRNA profiling of serum or plasma obtained from patients, the clinical relevance of many circulating circRNAs has been explored ( Table 1). Global circRNA expression can be profiled by using RNA sequencing (RNA-seq) followed by bioinformatic approaches [72]. circFARSA was identified as an upregulated circRNA in NSCLC tissues compared with adjacent normal tissue by analyzing backspliced reads on RNA-seq data [73]. circFARSA expression was higher in the plasma from NSCLC patients than in that from healthy volunteers and showed a good diagnostic value for NSCLC (AUC = 0.71). cDNA encoding circFARSA was cloned into the pLCDH-ciR vector, which was specifically designed to overexpress circular transcripts [55]. circFARSA overexpression enhanced the migration and invasion of A549 cells. Through in silico analyses, circFARSA was predicted to sponge miR-330 and miR-326 and regulate fatty acid synthesis. This is one of the earliest studies investigating the possibility of plasma circRNAs as new biomarkers for NSCLC patients; however, it lacks functional evidence supporting the molecular mechanism of circFARSA in NSCLC.
Through a microarray-based screening, circYWHAZ (circ_0005962) was identified as one of the upregulated circRNAs in lung adenocarcinoma (LUAD) [56]. The knockdown of circYWHAZ by siRNAs significantly suppressed the proliferation of LUAD cells, implying that this circRNA can promote cell proliferation [57]. Moreover, circYWHAZ expression was also upregulated in plasma samples, which illustrates a good diagnostic value for LUAD patients (AUC = 0.73). After surgical resection, circYWHAZ expression in the plasma decreased considerably, which suggests that circ_0005962 is potentially a good noninvasive biomarker for LUAD diagnosis [57]. miRNA-target prediction and functional enrichment analysis showed that circYWHAZ might function as a miRNA sponge to regulate LUAD development, which needs further validation.
circACP6 (circ_0013958) was also upregulated in LUAD tumors compared with nontumor tissues, which was validated by microarray and RT-PCR [58]. High expression of circACP6 was associated with the TNM stage (p = 0.009, Cox analysis) and lymphatic metastasis (p = 0.006) in LUAD patients. Moreover, the plasma expression levels of circACP6 distinguished LUAD from the control (AUC = 0.794, 95% CI = 0.703-0.912). Additionally, knockdown of circACP6 inhibited the proliferation, migration, and invasion of LUAD cells. Mechanistically, circACP6 functioned as a sponge against miR-134, which promoted the upregulation of cyclin D1, a target of miR-134. This study suggests that circACP6 might also be a novel biomarker for LUAD [58]. In this study, the authors showed the functional mechanism and diagnostic value of plasma circACP6. However, since this study was conducted with a relatively small number of patients (n = 30), a follow-up study with a larger, more diversified group of patients is needed.
circCXCR4 (circ_0056616) was identified and detected as a CXCR4-related circRNA in LUAD cells and exosomes [59]. Plasma exosome levels of circCXCR4 were lower in LUAD patients with TNM stage III-IV or with lymphatic metastasis than in those with stage I-II or without metastasis, respectively. This suggests that circCXCR4 might suppress the progression and metastasis of LUAD. Indeed, plasma exosomal circCXCR4 represents a good biomarker to diagnose lymphatic metastasis of LUAD (AUC = 0.812, 95% CI = 0.720-0.903), which also needs to be validated in a larger group of patients.
Through the exoRBase database (http://www.exorbase.org; accessed on 3 June 2020), circSATB2 (circ_0008928) was selected as a highly expressed circRNA in cancer exosomes [60]. The expression of circSATB2 was higher in lung cancer cells than in normal bronchial epithelial cells. Furthermore, overexpression and knockdown experiments showed that circSATB2 promoted the proliferation, migration, and invasion of lung cancer cells. Additionally, the packaging and transfer of circSATB2 by exosomes influenced the proliferation and migration of the recipient cells. circSATB2 directly bound to and inhibited miR-326, which in turn upregulated FSCN1, the presence of which has been reported as a poor prognostic marker for NSCLC patients [61]. Therefore, upregulation of FSCN1 by circSATB2 via sponging miR-326 represents a potential mechanism through which circSATB2 promotes NSCLC progression. In addition, serum exosomal circSATB2 expression was higher in NSCLC patients with metastasis than in those without, demonstrating a good diagnostic value for metastatic NSCLC (AUC = 0.797, 95% CI = 0.698-0.896). This study clearly showed that circSATB2 participated in NSCLC progression and could be a potential diagnostic marker for NSCLC.
In contrast, RNA-seq profiling demonstrated that circ_0102537 was one of the downregulated exosomal circRNAs in LUAD, which was also retrieved from a microarray database (GSE101586). Moreover, circ_0102537 was confirmed by quantitative RT-PCR to be downregulated in both plasma exosomes and tissues from LUAD patients. circ_0102537 knockdown by siRNAs promoted the migration and invasion of lung cancer cells and enhanced the expression of EMT markers such as N-cadherin, Snail, and Vimentin. This suggests that circ_0102537 might function as a tumor suppressor; however, the functional mechanism has not been presented [62]. Although many circulating circRNAs have been linked to lung cancer progression so far, further validation with more diverse groups of patients and in-depth mechanistic studies should be performed.

Anticancer Drug Response
Over a long period, numerous studies have been conducted to find predictive markers for sensitivity to EGFR inhibitors [74], and several circulating circRNAs have been proposed as candidate markers. Microarray analysis of plasma RNAs from NSCLC patients sensitive or resistant to gefitinib, an EGFR inhibitor, revealed that 1377 circRNAs were differentially expressed between the two groups [66]. Among them, circZNF91 (circ_0109320) was upregulated in the gefitinib-sensitive group. The plasma levels of circZNF91 could distinguish the gefitinib-sensitive group from the resistant group (AUC = 0.8054) and were associated with better progression-free survival in NSCLC patients treated with this EGFR inhibitor. Overall, circZNF91 could be a predictive biomarker of the sensitivity to gefitinib treatment in NSCLC patients after comparative verification with other parameters in a wider and larger group of patients. circC3 (circ_0002130) increased in NSCLC cells that acquired resistance to the EGFR tyrosine kinase inhibitor, osimertinib [67]. circC3 knockdown inhibited proliferation, glycolysis, and tumor growth in osimertinib-resistant lung cancer cells [68]. circC3 acted as a sponge against miR-498 to upregulate its targets, GLUT1, HK2, and LDHA, which are glycolysis-related proteins. Furthermore, an increase in circC3 was detected in serum exosomes from osimertinib-resistant NSCLC patients with respect to those from osimertinibsensitive patients. circC3 provided a good diagnostic value to predict the efficacy of osimertinib treatment in NSCLC patients (AUC = 0.792, 95% CI = 0.676-0.909), suggesting circulating circC3 as a novel biomarker. A combination of two or more circulating circRNAs with other variables such as EGFR mutations and gene copy number [74] would be a better biomarker for predicting the sensitivity to EGFR inhibitors.
circCNIH4 (circ_0000190) and circSHPRH were identified by RNA-seq to be upregulated in lung cancer cells compared with normal bronchial epithelial cells. They were also detected in conditioned media from lung cancer cells and in blood plasma samples by droplet digital PCR [63]. Furthermore, their plasma levels exhibited a poor response to immunotherapy, which might be due to the upregulation of soluble PD-L1 caused by these circRNAs [64]. Even though the detailed mechanism underlying the interplay between these circRNAs and antitumor immunity is still elusive, along with PD-L1 expression, their plasma levels could predict immunotherapy efficacy in lung cancer patients.
In-depth analysis of two GEO microarray datasets (GSE101684 and GSE101586) identified circRNAs highly expressed in LUAD samples compared with normal tissues [69]. Among them, circ_002178 promoted PD-L1 expression via sponging miR-34a. circ_002178 was also highly detected in plasma exosomes from LUAD patients compared with those from healthy volunteers, and exosomal circ_002178 had a significant diagnostic value for LUAD (AUC = 0.9967). Intriguingly, circ_002178 was transferred from cancer cells to CD8 + T cells via exosomes and then promoted PD-1 expression via sponging miR-28-5p. This indicates that circ_002178 would be a good target for immunotherapy, since it can modulate the expression of PD-1/PD-L1 in LUAD. As shown in this section, circulating circRNAs are highly likely to be exploited as markers for predicting the responses to anticancer drugs once their mechanisms of action are confirmed and their efficacy is validated in more diverse patients.

Cancer Diagnosis and Prognosis
As noted previously, numerous circulating circRNAs (circFARSA, circYWHAZ, cir-cACP6, circSATB2, circZNF91, circC3, and circ_002178) have significant diagnostic value and are associated with prognosis in lung cancer patients. In addition, RNA-seq and subsequent RT-PCR validation confirmed that circCD226 (circ_0047921) and circRALB (circ_0056285) were downregulated, while circATXN7 (circ_0007761) was upregulated in serum exosomes from NSCLC patients [70]. The combination of these three circRNAs provides a noteworthy diagnostic tool, which distinguishes NSCLC from healthy control (AUC = 0.919, 95% CI = 0.877-0.962) or other lung diseases, and their expression levels were associated with NSCLC progression. However, this study had several limitations such as insufficient sample sizes, samples from a single ethnic population, and lack of the mechanism of action of these exosomal circRNAs.
Plasma circCNIH4 demonstrated diagnostic potentials in lung cancer patients at all TNM stages (AUC = 0.95 for stage I-IV, AUC = 0.896 for stage I-II, and AUC = 0.96 for stage III-IV) [63]. Patients with high plasma levels of circCNIH4 exhibited poorer overall survival rates than those with low levels. Mechanistically, circCNIH4 could modulate the EGFR/ERK pathway by sponging miR-142-5p [65]. circPVT1 was also upregulated in tissues and sera from NSCLC patients. Both tissue and serum levels of circPVT1 showed diagnostic potential, distinguishing NSCLC patients from controls (AUC = 0.803 and 0.794, respectively) [71]. The knockdown of circPVT1 by siRNAs suppressed proliferation, migration, and invasion and promoted apoptosis in lung cancer cells. circPVT1 facilitated E2F2 signaling by functioning as a sponge against miR-125b. Even though the authors did not present the effect of circPVT1 on the survival or prognosis of NSCLC patients recruited in this study, they proved that circPVT1 can be used as a diagnostic marker and elucidated its working mechanism in NSCLC.
The expression levels of several circulating circRNAs are associated with major mutations found in lung cancer. For example, F-circEA, but not its host linear mRNA, could be detected in EML4-ALK + lung cancer plasma; thus, circulating F-circEA would be a novel biomarker to detect EML4-ALK fusion and to determine an effective treatment for EML4-ALK + patients [54]. LUAD patients with high plasma expression of circBNC2 (circ_0086414) were revealed to harbor EGFR mutations more frequently than those with low expression (p = 0.001) [57], suggesting that plasma circBNC2 would be a companion diagnostic marker for EGFR tyrosine kinase inhibitors. Considering the examples described above and the stable structure of circRNAs, circulating circRNAs can be novel biomarkers for diagnosis, prognosis, and treatment monitoring in lung cancer patients.

Conclusions and Perspectives
Since circRNAs have only recently started to attract attention, much is still unknown about their biogenesis and mechanisms of action. Nevertheless, the functional diversity and broad application potential of circRNAs offer new opportunities in cancer diagnosis and prognosis. In particular, the circulating circRNAs presented in this paper have great potential versatility as candidates for non-invasive biomarkers. Unlike mRNAs and miR-NAs [75], the heightened stability of circRNAs derived from their intrinsic structure further supports their applicability in liquid biopsy. However, there are several problems and limitations that must be addressed before their practical clinical application.
Basically, the expression levels of circulating circRNAs are quite low, so amplification and detection techniques with high sensitivity and accuracy are required. Techniques currently used for detecting circulating tumor DNA based on next-generation sequencing, digital-PCR, real-time PCR, or mass spectrometry could be applied to circulating circRNA detection [76]. Upregulation of plasma circCNIH4 and circSHPRH in lung cancer patients was detected using droplet digital PCR [63] which can quantitate nucleic acids with high sensitivity and accuracy [77]. The NanoString nCounter ® platform, which is a molecular barcoding system with target-specific oligonucleotide probes, was also used for circRNA quantification [78]. The tethered cationic liposome nanoparticle biochip has shown high sensitivity and selectivity for exosomal miRNA detection in small volumes of patients' serum [79], so it can also be applied for circRNA detection.
Despite the functional versatility and practical applicability of circulating circRNAs, there are still challenges to overcome before they can be applied in clinical settings, as has occurred for other liquid biopsy markers [75]. Since the exact mechanisms of action have not yet been elucidated, the expression levels of circRNAs that will actually lead to meaningful functional aspects are unknown. Depending on the purposes and scopes of application, it is necessary to determine whether qualitative or quantitative, or absolute or relative levels of circulating circRNAs are critical. The association of circRNAs with previously well-defined prognostic and diagnostic variables such as the Lung Cancer Prognostic Index (including stages, histology, mutation status, performance status, etc.) [80] should be considered in more depth. The relationship between circRNAs and key signaling players (e.g., KRAS, EGFR, BRAF, LKB1, MET, etc.) [81], which are closely related to lung cancer progression, should be studied in more detail.
Essentially, circulating circRNAs will face the same limitations and problems that previous candidate markers for liquid biopsy have encountered [82]. Firstly, since these cancer-associated markers are present in body fluids at very low levels, fast, cheap, and ultrasensitive detection techniques are required. The detection of circRNAs in exosomes or body fluids is also more time-consuming than assays based on conventional protein tests, which may limit the widespread use of circRNAs as biomarkers in clinical applications [51]. Proper samples should be accurately collected under tightly controlled environments with well-organized, systematic, efficient, and detailed protocols. Further, it should be clearly distinguished whether cancer cells or normal cells are the sources of the biomarkers. Some circRNAs can be expressed differently in tissues, although their serum levels are not significantly different. As in the case of alpha-fetoprotein in hepatocellular carcinoma, one of the most intensively studied cancer biomarkers [83], even though circRNA candidates are upregulated or downregulated in cancer tissues compared with non-cancer tissues, the difference between their plasma or serum levels in patients and healthy controls may not be apparent. These circRNAs are not suitable for diagnosing or predicting diseases in a non-invasive manner [51]. To be successfully integrated into clinical practice, the entire process from sample acquisition to marker analysis should be standardized and validated.
As mentioned above, various circulating circRNAs control the progression of lung cancer by sponging cancer-related miRNAs or by affecting the expression of cancer-related genes. In order to apply circulating circRNAs to the treatment of lung cancer, additional research on the detailed mechanisms through which circRNAs inhibit cancer progression is necessary. Furthermore, even if treatment targets and the related pathways are identified, biological safety should be guaranteed for their practical application in cancer treatment. For example, nanoparticles are convenient for carrying circRNA plasmids or circRNAtargeting agents in animal models, but it has not been fully revealed how safe they are for clinical use [84]. In addition, synthetic circRNAs can boost the immune system in vivo because foreign circRNAs can be distinguished from endogenous circRNAs due to their lack of N6-methyladenosine modification [85]. There are still many obstacles in developing a treatment for lung cancer using circulating circRNAs; therefore, their application for diagnosis would be prioritized over their use in the treatment of lung cancer. Despite all these challenges, circulating circRNAs have advantages as promising biomarkers and can be actively used for lung cancer diagnosis, prognosis, and treatment monitoring.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare no conflict of interest.