Comparison of colonial volvocine algae based on phylotranscriptomic analysis of gene family evolution and natural selection

ABSTRACT

This study is the first to determine the transcriptomes of eight colonial volvocine algae: Pandorina morum, P. colemaniae, Volvulina compacta, Eudorina elegans, E. cylindrica, Yamagishiella unicocca, Colemanosphaera charkowiensis and C. angeleri. By analysing these data together with the whole-genome data of Chlamydomonas reinhardtii and other colonial volvocine algae, we determined 26 979 gene families and 105 single-copy orthogroups. Using the single-copy orthogroups, we adopted a supermatrix and coalescent-based approach to perform phylotranscriptomic analysis of colonial volvocine algae. The phylogenetic relationships of most species were consistent with those of a previous study, and the results showed paraphyly of the genus Pandorina. Then, we analysed gene family expansion and contraction in colonial volvocine algae and found that they exhibited a large number of expanded and contracted gene families. In particular, the common ancestor of the Volvocaceae had the largest number of expanded gene families, which may be associated with the increased volume of extracellular matrix. Finally, using the ML method to analyse the substitution rates among colonial volvocine algae based on single-copy orthogroups, we showed a significantly higher non-synonymous substitution (dN) rate in Tetrabaena socialis. By analysing the differences in the dN/dS ratio, we found that more orthogroups exhibited a significantly higher dN/dS in T. socialis; we speculate that the higher dN/dS could indicate relaxed selection pressure and that this lower constraint allows the accumulation of more non-synonymous mutations in T. socialis.

KEYWORDS:

• Coalescent-based approach
• colonial volvocine algae
• gene family evolution
• phylotranscriptomic
• relaxed purifying selection
• substitution rates

Introduction

The transition from unicellular organisms to multicellularity provided advantages such as better utilization of scarce resources, new relationships between cells and the environment, and escape from predators (Umen, 2014). The transition to multicellularity has happened many times in different species (Knoll, 2011; Parfrey & Lahr, 2013), and is associated with considerations of the life cycle, cell biology and ecology of simpler ancestor(s) that gave rise to multicellular descendant lineages. As more whole genomes have been sequenced, comparative genomic analysis provides insights into the evolutionary transition. Regarding the origin of metazoans, choanoflagellates evolved before the origin and diversification of metazoans (Lavrov et al., 2005; Rokas et al., 2005). King et al. (2008) analysed the genome of the unicellular choanoflagellate Monosiga brevicollis, finding several protein domains exclusive to choanoflagellates and metazoans; these domains are important in, for example, cell adhesion, the extracellular matrix and cell signalling, suggesting an important role in the origin of multicellularity. Srivastava et al. (2010) compared the genome of Amphimedon queenslandica with related species and found that many gene families are associated with six hallmarks of metazoan multicellularity. In Protista, algae have acquired complex multicellularity several times independently (De Clerck et al., 2018). Within Streptophyta, the expanded transcription factors families, rare domain combinations, diverse protein domains and others all contribute to multicellularity in Ulva (De Clerck et al., 2018). Within the Chlorophyta, the colonial volvocine algae are an important model system to study the transition to multicellularity (Leliaert et al., 2012). Comparative genomic studies showed that the retinoblastoma cell cycle regulatory pathway and cyclin-dependent kinases are important to multicellularity (Hanschen et al., 2016), and expanded gene families and the ubiquitin proteasomal pathway are also associated with the origin of multicellularity (Featherston et al., 2017).

The colonial volvocine algae belong to the order Chlamydomonadales (sometimes referred to as Volvocales); they are all multicellular species with colonies composed of Chlamydomonas-like cells and surrounded by extracellular matrix. All these algae were initially assigned to the family Volvocaceae. Based on the tripartite boundary of the matrix (which appears as a dark line when stained and examined by electronic microscope; Coleman, 2012), Nozaki & Kuroiwa (1992) assigned the genera Gonium and Astrephomene to the family Goniaceae. Then, Nozaki & Itoh (1994) used cladistic analysis to deduce phylogenetic relationships within the colonial volvocine algae, and the four-celled species Gonium sociale (renamed Tetrabaena socialis) together with Basichlamys was placed in the new family Tetrabaenaceae. Therefore, the colonial volvocine algae are now composed of three families, Tetrabaenaceae, Goniaceae and Volvocaceae. The simplest colonial species are in the family Tetrabaenaceae, with only four cells forming a square colony (Nozaki & Ohtani, 1992; Nozaki & Itoh, 1994; Nozaki et al., 1996). Species in Goniaceae have a colony composed of more than four cells: species of Gonium form planar colonies of eight or 16 cells (Hayama et al., 2010); members of the genus Astrephomene form ovoid colonies of 32 or more cells (Nozaki, 1983). The family Volvocaceae is composed of eight genera (Pandorina, Volvulina, Platydorina, Yamagishiella, Colemanosphaera, Eudorina, Pleodorina, Volvox), with species forming a spheroidal colony of more than eight cells (Nozaki & Kuroiwa, 1992), the evolution of which has been investigated (Hallmann, 2011).

The taxonomy of colonial volvocine algae is mainly based on morphology and phylogenetic analysis. Regarding morphology, the colonial volvocine algae show a series of extant organisms with organizational complexity (Hallmann, 2011). Thus, their morphology makes colonial volvocine algae an ideal model for studying the evolutionary transition to complexity. Phylogenetic analyses have mainly used markers such as the internal transcribed spacer (ITS) (Angeler et al., 1999; Coleman, 2001; Yamada et al., 2008; Hayama et al., 2010; Nozaki et al., 2015), the large subunit of the ribulose-1,5-bisoxyg/oxygenase (rbcL) gene (Nozaki et al., 1997a, 1997b, 1998; Yamada et al., 2008; Hayama et al., 2010), and five chloroplast protein-coding genes (Nakada et al., 2010; Nozaki et al., 2000, 2006, 2014, 2015; Nozaki & Coleman, 2011). Though phylogenetic analyses based on one or a few genes have been widely performed in recent decades, they have some limitations, such as the existence of paralogues, different evolutionary rates, incomplete lineage sorting, horizontal gene transfer and gene duplication (Maddison, 1997; Huynen & Bork, 1998; Doolittle, 1999; Degnan & Rosenberg, 2006; Zou et al., 2016; Gernandt et al., 2018), meaning that the gene tree may not represent the species tree in some cases. Generally, increases in the information in a dataset and the use of a more appropriate analysis method reduce the influence of such errors (Wolf et al., 2002; Mirarab et al., 2014; Davidson et al., 2015; Mirarab, 2019). With the application of high-throughput sequencing, genome archives are now available to perform phylotranscriptomic analyses on different species. High-throughput sequencing would provide sufficient information for our dataset, and the incongruence of phylogenetic results would therefore disappear, resulting in a robust resolution. Though phylotranscriptomics have been employed for different species (Wickett et al., 2014; Johnson et al., 2016; Jackson et al., 2017), little is known about the phylotranscriptomics of colonial volvocine algae.

High-throughput sequencing also makes it possible to analyse gene family evolution. The expansion and contraction of gene families have been widely used to analyse the origin, adaptation and evolution of different species (Hahn et al., 2007; Kim et al., 2016). For colonial volvocine algae, previous studies have proposed that the expansion of gene families is involved in extracellular matrix formation and cell-cycle regulation (Prochnik et al., 2010; Hanschen et al., 2016; Featherston et al., 2017). Thus, gene family evolution analysis may help us understand the evolutionary transitions in this group of algae. There are five colonial volvocine species for which whole genomes are available, but the nucleotide and amino acid sequence annotation of Eudorina and Yamagishiella are currently restricted to sex-determining regions (Prochnik et al., 2010, Hanschen et al., 2016; Featherston et al., 2017; Hamaji et al., 2018). So, gene family evolution analysis has only been studied in three species (Tetrabaena socialis, Gonium pectorale and Volvox carteri) in colonial volvocine algae.

Natural selection plays an important role in shaping organismal diversity (Lenski et al., 2003); generally, most genes are under purifying selection to maintain their function (Zhang, 2003). However, changes in selection may lead to evolutionary innovations; for example, Zhen & Ungerer (2008) found that the freezing tolerance of Arabidopsis thaliana is associated with relaxed purifying selection on a group of genes. The dN/dS ratio is a measurement of natural selection acting on the genes. Normally, values of dN/dS < 1, = 1 and > 1 mean negative purifying selection, neutral evolution and positive selection, respectively (Yang, 2007). dN/dS analysis in colonial volvocine algae would further contribute to revealing the evolutionary relationships in this group of algae.

In this study, the transcriptomes of eight colonial volvocine algae were sequenced. By using these genome archives, we conducted phylotranscriptomic analysis among the colonial volvocine algae to obtain a robust resolution of the phylogenetic relationships between them. We also analysed the expansion and contraction of gene families among these species, which offers useful information about their adaptation. Finally, we analysed dN/dS in these algae, which provides further information on the evolutionary relationships in this group.

Materials and methods

Cultures

Eight colonial volvocine algae were obtained from the Culture Collection of Freshwater Algae at the Institute of Hydrobiology, Chinese Academy of Sciences. The strain numbers were as follows: Pandorina morum (strain FACHB-2362), P. colemaniae (strain FACHB-2361), Volvulina compacta (strain FACHB-2337), Colemanosphaera angeleri (strain FACHB-2363), C. charkowiensis (strain FACHB-2326), Yamagishiella unicocca (strain FACHB-2364), Eudorina elegans (strain FACHB-2321) and E. cylindrica (strain FACHB-2322). All algae were cultured in BG11 medium, and the cultures were grown at 20–25°C under a 14 h light:10 h dark schedule under cool-white fluorescent lamps at an intensity of 10–20 µmol m^–2 s^–1.

Library preparation and sequencing

Total RNA was extracted using TRIzol reagent (Thermo Fisher Scientific, USA), and poly-A+ mRNA was isolated using oligo dT magnetic beads. The mRNA was fragmented with divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5×) and was used as a template for random hexamer-primed first-strand cDNA synthesis, after which the second strand of cDNA was synthesized. The sequencing library was generated via the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, USA) and sequenced on the NovaSeq 6000 platform (Illumina, USA). Finally, 150 bp paired-end reads were generated.

Quality control, de novo assembly and sequence annotation

The quality of the raw reads was first checked by using FastQC v0.11.6 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and the raw reads were subjected to quality control with Trimmomatic v0.38 (Bolger et al., 2014) (LEADING:3, TRAILING:3, SLIDINGWINDOW:4:15, HEADCROP:10, MINLEN:36). Then, de novo assembly of the clean reads was conducted using Trinity v2.8.4 (Grabherr et al., 2011) with default parameters and CD-HIT v4.7 (Fu et al., 2012) was used to remove redundant transcripts with the default parameters. The completeness of the final transcripts was assessed by BUSCO v3.0.2 (Simão et al., 2015). Then, the TransDecoder v5.5.0 was used to predict the open reading frame (ORF) of each transcript. BLASTP (Camacho et al., 2009) searches of the longest ORFs (the longest coding region in each transcript) were conducted against the Uniref90 database with Diamond v0.8.22.84 (Buchfink et al., 2015) and Pfam searches of the longest ORFs in the Pfam database were conducted using HMMER v3.1b2 (Mistry et al., 2013). Finally, TransDecoder was used to integrate the BLASTP and Pfam search results into coding regions, and the nucleotide sequences (CDS) and amino acid sequences (PEP sequence) of these regions were used for the following analyses. Additionally, the CDS and PEP sequence of Chlamydomonas reinhardtii, Tetrabaena socialis, Gonium pectorale and Volvox carteri were downloaded from NCBI; the assembly accession numbers of each species were GCF_000002595.1, GCA_002891735.1, GCA_001584585.1 and GCF_000143455.1, respectively. All raw reads were deposited in the NCBI Sequence Reads Archive (BioProject PRJNA532307). The CDS and PEP sequences of Yamagishiella unicocca and Eudorina sp. are currently restricted to sex-determining regions, so these species were not included in our analysis.

Orthologous group identification and phylotranscriptomic analysis

The PEP sequences of all taxa were used to determine single-copy orthologues with OrthoFinder v2.2.7 (Emms & Kelly, 2015) and the PEP sequences of single-copy orthologues were selected for phylotranscriptomic analysis. Each single-copy orthologue was aligned using MAFFT v7.394 (Katoh & Standley, 2013) with the options -maxiterate 1000 and -globalpair. Regions showing poor alignment were trimmed with trimAl v1.2 (Capella-Gutiérrez et al., 2009) using the parameter -automated1. The trimmed alignment of orthologous groups was used for the following phylotranscriptomic analysis.

We used supermatrix analysis and coalescent-based analysis to construct the phylogenetic tree. For the supermatrix analysis, we first used Phyutility (Smith & Dunn, 2008) to concatenate all orthologous groups and the evolutionary models and partitioning of the PEP dataset were determined with PartitionFinder 2 (Lanfear et al., 2016). Phylogenies were inferred from the PEP dataset using Maximum likelihood (ML) and Bayesian (BI) methods. ML analyses were carried out using RAxML v8.2.10 (Stamatakis, 2014). Rapid bootstrap analysis with 1000 replicates of the dataset for ML was performed to estimate statistical reliability. Bayesian analysis was performed with MrBayes v3.2.6 (Ronquist et al., 2012). Markov chain Monte Carlo (MCMC) analysis was run with four Markov chains (three heated, one cold) for 1 000 000 generations with trees sampled every 1000 generations. Each time the diagnostics were calculated, a fixed number of samples (burnin = 1000) were discarded from the beginning of the chain. For the coalescent-based analyses, ML analysis of each single-copy orthologue was performed in RAxML by conducting 1000 rapid bootstrap replicates with the PROTGAMMAGTR model. Then, the resulting best trees were used to infer the coalescence-based species tree (ST) phylogeny with ASTRAL v5.6.3 (Zhang et al., 2018).

Gene family evolution analysis

The MCMCtree Bayesian program in PAML v4.9 (Yang, 2007) was used to estimate divergence times under relaxed-clock models; these analyses were conducted on the CDSs, and the fossil calibration information mainly referred to Herron et al. (2009) with RootAge specified as ‘>200<300’. The phylogenetic ML tree, BI tree and ST tree all showed the same topology, which was used for divergence time estimation. Based on the estimations, the CAFE v3.1 (De Bie et al., 2006) program was used to calculate the expansion and contraction of a gene family; we used the criterion of p < 0.05 for determining significantly changed gene families. All PEP sequences were imported into InterProScan v5.32 (Philip et al., 2014) for Gene Ontology (GO) term mapping (http://www.geneontology.org/). GO enrichment analysis was performed using clusterProfiler v3.12.0 (Yu et al., 2012), and we controlled for multiple testing using the false discovery rate method (Benjamini & Hochberg, 1995).

dN/dS analysis

The CODEML program of PAML v4.9 (Yang, 2007) was used to estimate substitution rates and dN/dS ratios. The ML model (runmode = −2, CodonFreq = 2) was used to measure non-synonymous (dN) and synonymous (dS) substitution rates; the analysis was based on all single-copy orthologues, Chlamydomonas reinhardtii was used as a reference and orthologues with dS values > 5 were discarded from further analysis. Comparisons of evolutionary rates used the Wilcoxon rank sum test. The phylogenetic tree was used as a constraint tree, but branch lengths were inferred by using PAML. The ML model showed significantly higher dN rates in Tetrabaena socialis than in other colonial volvocine algae. To explore the reason for this phenomenon, we used the branch model (two-ratio; Yang, 1998) to test whether T. socialis presents a different dN/dS relative to the other species. T. socialis was labelled as the foreground branch. A null model (model = 0) in which dN/dS ratio was fixed across all species was compared with an alternative model (model = 2) in which T. socialis was allowed to have a different dN/dS. Likelihood ratio tests (LRT) were used to test the model fit and the chi-square test was applied for testing p values. Multiple testing was corrected by applying the false discovery rate method (FDR) (Benjamini & Hochberg, 1995) as implemented in R (http://www.R-project.org). GO annotation of the T. socialis sequence was based on InterProScan v5.32 (Philip et al., 2014).

Results

De novo transcriptome assembly and orthology detection

The related sequencing and assembly information is listed in Table 1. Illumina paired-end sequencing technology revealed that the number of raw reads yielded by each species ranged from 50362006–71 111684. After filtering for adaptors and low-quality sequences, the number of clean reads ranged from 48 284776–69311428. Then, Trinity was employed for de novo assembly; the number of assembled contigs ranged from 60366–130 254 and the average contig length and N50 values of all species were greater than 1100 bp and 2200 bp, respectively, which is indicative of high contiguity. After removing redundant transcripts via CD-HIT, BUSCO software was applied for further quantitative assessment of assembly and annotation completeness. Searching for the BUSCOs defined for Chlorophyta revealed that the percentages of conserved genes in our eight species transcriptomes recovered via BUSCO analysis were >86%, indicating a high completeness level. Finally, TransDecoder was used to predict coding sequences from the assembled transcripts; the number of coding sequences ranged from 39 532 to 81 820, and the nucleotide sequence (CDSs) and amino acid sequence (PEP sequence) of the coding sequences were used for the following analysis.

Table 1. Summary of sequencing and assembly.

			Sequence assembled by Trinity
Species	Numberof raw reads	Number of clean reads	Number of contigs	Average contig length	N50	Number of contigs after filtered by CD- HIT	Complete BUSCOs	Number of coding sequence predicted by TransDecoder
Pandorina morum	59675614	57558362	130254	1453	3003	109148	94.80%	81820
Pandorina colemaniae	50362006	48284776	64700	1411	2711	56751	88.90%	39532
Volvulina compacta	55086070	53180216	76435	1780	3855	67035	94.00%	49013
Colemanosphaera angeleri	56239620	54588034	105381	1447	3231	92967	90.90%	59379
Colemanosphaera charkowiensis	55935452	53811642	95450	1535	3189	83987	90.50%	59008
Yamagishiella unicocca	61845914	59687248	60366	1682	3211	51406	92.70%	42748
Eudorina elegans	71111684	69311428	103779	1103	2252	94254	87.70%	59656
Eudorina cylindrica	52884220	50672466	93611	1329	2532	81562	86.60%	64418

The CDSs and PEP sequences of C. reinhardtii, T. socialis, Gonium pectorale and Volvox carteri were downloaded from the NCBI database. Together with our eight species, we used the PEP sequences of 12 species in total to detect orthology. The total number of genes that we used was 515 100, and only 69.3% of the genes were assigned to orthogroups. The number of orthogroups was 26 979, where each orthogroup represented a gene family, and all orthogroups were used for the gene family evolution analysis. The number of single-copy orthogroups was 105, and all single-copy orthogroups were used for phylotranscriptomic analysis and dN/dS analysis.

Phylotranscriptomic analysis

The phylogenetic tree (Fig. 1) was constructed for all 12 taxa via the PEP sequences of single-copy orthogroups. The ML tree, BI tree and ST tree all showed the same topology reflecting the relationships among the colonial volvocine algae. When rooted with C. reinhardtii, T. socialis is sister to the clade formed by G. pectorale and other colonial volvocine algae, and G. pectorale is sister to the clade formed by other colonial volvocine algae. T. socialis belongs to family Tetrabaenaceae, G. pectorale to the Goniaceae and all other colonial volvocine algae are in the Volvocaceae. The relationships between these three families in our phylogenetic tree are consistent with those of previous studies mainly based on morphological data, ITS or five chloroplast genes (Nozaki & Itoh, 1994; Coleman, 1999; Nozaki et al., 2000; Herron & Michod, 2008), further indicating their validity.

Fig. 1. Phylogenetic tree of the colonial volvocine algae based on 105 single-copy orthogroups. Numbers on branches represent Bayesian posterior probabilities, bootstrap values and support values of ASTRAL, respectively. Branch lengths are proportional to genetic distances, which are indicated by the scale bar. The grey background indicates the outgroup; the yellow background indicates family Tetrabaenaceae; the blue background indicates family Goniaceae; and the green background indicates family Volvocaceae.

All of the Volvocaceae constitute a larger monophyletic group; in this clade, the phylogenetic position of all species is consistent with the study of Nozaki et al. (2014) (mainly based on five chloroplast genes) except for the Pandorina–Volvulina clade. Our phylotranscriptomic analysis showed that P. morum and V. compacta form a clade and that this clade is sister to P. colemaniae; this phylogenetic relationship has previously been reported, largely based on ITS or five chloroplast genes (Coleman, 2001; Herron & Michod, 2008; Herron, 2016; Hanschen et al., 2018).

Gene family evolution analysis

We used the birth-death model of CAFE to analyse gene family evolution, and this method assumes at least one gene in the root of the species tree. Among all 26 979 gene families, only 9497 gene families met this criterion, so we only used 9497 gene families to analyse expansion and contraction.

As shown in Fig. 2, the gene families of all colonial volvocine algae have experienced a large number of expansions and contractions; Volvox carteri displayed 1144 contracted gene families, which is greater than the number in all other nodes. Referring to the ancestor node, the common ancestor of Volvocaceae displayed 2112 expanded gene families, more than for all other nodes.

Fig. 2. Gene family expansion or contraction in colonial volvocine algae. Branch numbers indicate the number of gene families that have expanded (+) and contracted (−) since the split from the common ancestor. The timelines indicate divergence times among the species.

To investigate whether the observed gene family evolution is associated with morphological diversity, we performed Gene Ontology (GO) enrichment analysis based on the expanded gene families. For the common ancestor of Volvocaceae, the expanded families were enriched in 202 GO terms, among which 21 GO terms belonged to the cellular component category (Fig. 3). These GO terms were mainly associated with different biological complexes; in particular, the expanded gene families of the common ancestor of Volvocaceae were enriched in the term ‘extracellular matrix’ (GO:0031012, adjusted p value = 3.51E-09), which may be related to an increased volume of extracellular matrix.

Fig. 3. Dot plot showing the enrichment of cellular component GO terms for expanded gene families of the common ancestor of Volvocaceae (FDR adjusted p < 0.05). The dot sizes represent the numbers of genes.

dN/dS analysis

The ML method was used to calculate the rate of non-synonymous (dN) and synonymous (dS) substitutions for 105 single-copy orthogroups between the colonial volvocine algae; after discarding the orthogroups with dS >5, there were 20 orthogroups. The boxplot of each species based on the dN and dS rates of these 20 orthogroups is shown in Fig. 4. From the boxplot, we can see that the dN rate of T. socialis was higher than those of other species, and there were no significant differences between dS rates of different species. We compared the dN rate of T. socialis with those of all other colonial volvocine algae using the Wilcoxon rank sum test and it was significantly different (p = 9.844e-05).

Fig. 4. Boxplots of the non-synonymous (dN) and synonymous (dS) substitutions for each species of colonial volvocine algae. For each species, the box represents values between quartiles, outliers are shown as black points, and black lines inside the box show median values.

We used the branch model (two-ratio; Yang, 1998) to test the differences in the dN/dS ratio between T. socialis and other colonial volvocine algae. The results are shown in Table 2. Five orthogroups showed significant differences in their dN/dS ratios (FDR corrected p value <0.05), among which four orthogroups had significantly higher dN/dS ratios in T. socialis. By using InterProScan, we found that these orthogroups mapped to the following GO terms: nucleus (GO:0005634), chromatin assembly or disassembly (GO:0006333), structural constituent of ribosome (GO:0003735), ribosome (GO:0005840) and translation (GO:0006412). All the GO terms of these orthogroups were mainly associated with the process of protein forming.

Table 2. Differences in dN/dS between Tetrabaena socialis and other colonial volvocine algae. Orthogroups with FDR-adjusted p values < 0.05 are shown in bold.

Orthogroups	dN/dS of other colonial volvocine algae	dN/dS of Tetrabaena socialis	df	2LnL	p value	FDR-corrected
OG0010594	0.1142	0.0920	1	0.5212	0.4703	0.5856
OG0010642	0.0626	0.0988	1	2.0005	0.1572	0.3145
OG0010651	0.0185	0.1945	1	25.2248	0.0000	0.0000
OG0010676	0.0302	0.0553	1	1.7255	0.1890	0.3436
OG0010733	0.0539	0.1055	1	3.7270	0.0535	0.1530
OG0010753	0.0930	0.2105	1	8.5655	0.0034	0.0137
OG0010764	0.0174	0.0386	1	3.2686	0.0706	0.1765
OG0010774	0.0942	0.0975	1	0.0114	0.9150	0.9632
OG0010775	0.0255	0.0585	1	4.5023	0.0339	0.1128
OG0010790	0.0234	0.0132	1	15.4574	0.0001	0.0006
OG0010840	0.0188	0.1876	1	24.5617	0.0000	0.0000
OG0010856	0.1070	0.1092	1	0.0021	0.9634	0.9634
OG0010860	0.0375	0.0706	1	2.1768	0.1401	0.3113
OG0010864	0.0365	0.0550	1	0.8699	0.3510	0.4995
OG0010901	0.1504	0.2254	1	1.3298	0.2488	0.4147
OG0010904	0.0570	0.0822	1	1.1466	0.2843	0.4373
OG0010926	0.0101	0.0848	1	14.0079	0.0002	0.0009
OG0010939	0.1409	0.1881	1	0.4598	0.4977	0.5856
OG0010975	0.0355	0.0416	1	0.2229	0.6368	0.7076
OG0010995	0.0844	0.1181	1	0.7882	0.3746	0.4995

Discussion

In a recent study, Proeschold et al. (2018) showed that C. reinhardtii was nested within the colonial volvocine clade, casting doubts on choosing it as an outgroup for the colonial volvocine algae clade. Their analysis was based on three datasets: (1) rbcL, (2) ITS and (3) SSU and ITS. However, Zhang et al.’s (2019) phylotranscriptomic analysis showed that C. reinhardtii was sister to the colonial volvocine algae clade. Their analysis was based on supermatrix and coalescent approaches using several unlinked orthologous genes, an approach similar to ours. Two datasets in Proeschold et al. (2018) were based on single gene like ITS or rbcL, but the evolutionary relationships of the species being investigated may not be reflected by a single gene tree (Maddison, 1997; Pamilo & Nei, 1988; Wolf et al., 2002). Previous studies proposed three main limitations for single-gene phylogenetic analysis: insufficient number of informative sites; different lineages having different evolutionary rates; and the impact of lateral gene transfer. These problems could be overcome using the supermatrix approach (Wolf et al., 2002). The phylotranscriptomic analysis in Zhang et al. (2019) and our analysis were based on the concatenated dataset of single-copy orthologues, genes related by speciation (Koonin, 2005). By using single-copy orthologues we can eliminate the impact of gene duplication and loss (Mallo & Posada, 2016). In our analysis, the concatenated data have been partitioned to apply different evolutionary models to different partitions, which would further overcome the limitations of the supermatrix approach (Wolf et al., 2002). Proeschold et al. (2018) also proposed that gene sequences should not be combined in a concatenated dataset if the tree topologies of single genes differ significantly. One of the main reasons for gene tree discordance is the issue of incomplete lineage sorting (Degnan & Rosenberg, 2009). To overcome this limitation, Zhang et al. (2019) and our study both used coalescent-based analyses using the ASTRAL program. ASTRAL enables highly accurate phylogenomic estimation, even in the presence of high levels of gene tree conflict because of incomplete lineage sorting (Mirarab et al., 2014) or horizontal gene transfer (Davidson et al., 2015; Mirarab, 2019). Meanwhile, the phylogenetic analysis of Proeschold et al. (2018) based on the concatenated dataset of SSU and ITS also has limitations: the SSU and ITS are genes or regions in physical proximity, thus the concatenated dataset is a linkage group. Instead, phylotranscriptomic analysis based on single-copy orthologues uses a large number of unlinked loci, which is an appropriate way to reduce the incongruence between gene trees and species trees (Hillis, 1995; Ruvolo, 1997; Sang, 2002). We consider that all these are a priori reasons to follow Zhang et al. (2019) rather than the Pröschold et al. (2018) analysis. Hence, we still consider C. reinhardtii as an outgroup for colonial volvocine algae.

The genera Pandorina and Volvulina are reported to be non-monophyletic (Coleman, 2001; Nakada et al., 2010; Nozaki et al., 2014), but the phylogenetic analyses in these studies have mainly been based on only one or a few genes. In this study, we used 105 single-copy orthogroups to construct the phylogenetic tree. The total amino acid length of single-copy orthogroups was 22 936, and an increased sequence length should improve the support, resolution and accuracy of the phylogenetic analysis (Wortley et al., 2005). Additionally, we used supermatrix and coalescent-based analyses to perform the phylotranscriptomic analysis, two methods each presenting certain advantages as mentioned above; both methods showed paraphyly of Pandorina. Our results showed that we could not resolve the paraphyly of Pandorina by increasing the sequence length or changing the analysis method. In this study, we used only two Pandorina species and one Volvulina species, and this lack of resolution may have mainly been due to insufficient analysed species.

The number of expanded or contracted gene families may be influenced by different factors, such as gene duplication, de novo gene creation, gene loss, the function of gene families and changes in environmental conditions (Prachumwat & Li, 2008; Guo, 2013; Albalat & Cañestro, 2016). Therefore, the different numbers of gene family expansions and contractions observed in colonial volvocine algae may be due to a variety of causes. Thus far, whole genomes have been made available for seven colonial volvocine algae, with genome sizes and gene family evolution (number in this study, noted in parentheses) as follows: Tetrabaena socialis, 135.78 Mb (+62/−73); Gonium pectorale, 148.81 Mb (+69/−143); Yamagishiella unicocca, 134.23-Plus, 140.84-Minus (+726/−57); Eudorina sp. Male-168.62 MB, Female-184.03 MB (+1427/−90 or +781/−99) and Volvox carteri, 137.68 Mb (+46/−1144). A previous study proposed that gene family contraction is one of the reasons for genome reduction (Qiu et al., 2015). In our case, although V. carteri shows the largest number of gene family contractions, it does not have the smallest genome size – its genome is the third smallest. We speculate that the large number of gene family contractions could be only one of the reasons for the reduced genome size of V. carteri.

Prochnik et al. (2010) compared the genomes of V. carteri and C. reinhardtii, finding that two gene families involved in extracellular matrix functions are substantially larger in the V. carteri genome. As V. carteri has a much greater volume of extracellular matrix than C. reinhardtii, these expanded gene families may be an important factor in the size of the extracellular matrix in colonial volvocine algae (Leliaert et al., 2012). For the common ancestor of Volvocaceae, we also found that the expanded gene families were enriched in extracellular matrix related genes. It has been proposed that the common ancestor of Volvocaceae experienced an increase in the volume of extracellular matrix (Kirk, 2005; Herron et al., 2009); as gene family expansion may be associated with the size of extracellular matrix as mentioned above, and the expanded gene families of the Volvocaceae common ancestor are enriched in extracellular matrix related genes, we propose that the common ancestor of the Volvocaceae experienced gene family expansion, which may have contributed to increasing the volume of extracellular matrix.

We also compared the substitution rate among the colonial volvocine algae and found that dN was significantly higher in Tetrabaena socialis. The dN/dS ratios of most orthogroups showed no difference, but T. socialis exhibited four more orthogroups with significantly higher dN/dS ratios. Although the dN/dS ratios of all orthogroups were relatively low (dN/dS <1), which indicated strong purifying selection (Yang, 2007), the relatively higher dN/dS ratio could be interpreted as indicating relaxed purifying selection (Guisinger et al., 2008). The greater number of orthogroups with higher dN/dS ratios in T. socialis could be considered indicative of relaxed purifying selection in this species; this relaxed selection pressure may allow the accumulation of additional non-synonymous substitutions in T. socialis (Björnerfeldt et al., 2006), leading to higher dN substitution rates. The higher dN substitution rates may result in a novel function, for example, Kato-Minoura et al. (2015) studied two actin genes among species in Volvocales (including colonial volvocine algae), a conventional actin (actin) gene and an unconventional action (NAP) gene. They showed a similar dS between the two genes but higher dN in NAP, and proposed that NAP is unique to volvocine species and performs a cellular function with fewer constraints; more importantly, NAP genes might have a function different from that of the actin gene. As the orthogroups under relaxed purifying selection in T. socialis are mainly associated with the process of protein formation, such a phenomenon may foster evolutionary innovation (Wertheim et al., 2014).

Only a limited number of genome data for colonial volvocine algae are available currently, so several aspects of further study are needed. (1) The transcriptome data of more strains are needed to fully resolve the phylogenetic relationships of the Pandorina–Volvulina clade; (2) More well-annotated whole genomes are needed to find out the main reason for the reduced genomic size of V. carteri; (3) Among the colonial volvocine algae, why would T. socialis have experienced relaxed purifying selection and what would be the result of this relaxed constraint?

Our study is the first to perform phylotranscriptomic analysis based on transcriptome data from colonial volvocine algae. By using 105 single-copy orthogroups and two phylogenetic analysis approaches, we found paraphyly of the genus Pandorina yet again, and we propose that more species are needed for phylogenetic analysis to fully resolve this phenomenon. Based on gene family evolution analysis, we found that the expanded gene families of the common ancestor of the Volvocaceae may be associated with an increased volume of extracellular matrix. When compared with other colonial volvocine algae, we found that some single-copy orthogroups of T. socialis have experienced relaxed purifying selection, which allows for significantly higher non-synonymy in T. socialis.

Author contributions

Y. Hu: original concept, culture experiments, analysis of molecular data, drafting manuscript; All authors: manuscript editing.

Acknowledgements

The research was supported by the Wuhan Branch, Supercomputing Center, Chinese Academy of Sciences, China.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the Special Foundation of Science and Technology Basic Work of China [grant no. 2014FY210700], the National Natural Science Foundation of China [grant no. 31670202] and China Agriculture Research System [CARS-50].