The strains described in this study were isolated from water samples and deposited in the Freshwater Algae Culture Collection at the Institute of Hydrobiology (FACHB collection), Wuhan, Hubei Province, China. The Colemanosphaera charkowiensis (strain FACHB-2326) were collected from a river in Weihe (36°19′7″ N, 119°28′54″ E), Weifang, Shandong Province, China, in September 2017. The Volvulina compacta (strain FACHB-2337) were collected from a river in Tangxihe (31°7′13″ N, 108°49′10″ E), Chongqing, China, in June 2017. The Pandorina colemaniae (strain FACHB-2361) were collected from a river in Wuhe (36°22′25″ N, 119°24′52″ E), Weifang, Shandong Province, China, in September 2017. The Pandorina morum (strain FACHB-2362) were collected from a river in Weihe (36°8′1″ N, 119°25′59″ E), Weifang, Shandong Province, China, in September 2017. The Colemanosphaera angeleri (strain FACHB-2363) were collected from a river in Weihe (36°30′2″ N, 119°24′44″ E), Weifang, Shandong Province, China, in September 2017. The Yamagishiella unicocca (strain FACHB-2364) were collected from a pool in Dichi (47°18′20″ N, 120°28′38″ E), Aershan, Inner Mongolia, China, in August 2017. The strains were grown in a conical flask containing artificial freshwater-6 (AF-6) (Kato, 1982) at 20–25°C under a 14 h light: 10 h dark schedule under cool-white fluorescent lamps at an intensity of 1000–2000 lux. Total genomic DNA was extracted using a Universal DNA Isolation Kit (AxyPrep, Suzhou, China) following the manufacturer’s instructions. Species identification was based on morphological observation and phylogenetic analysis based on five chloroplast genes (rbcL, atpB, psaA, psaB, and psbC), and the sequence data of related species were chosen according to Nozaki et al. (2014). The sequence matrix was aligned by MAFFT v7.394 (Katoh and Standley, 2013), and the ambiguously aligned regions were further manually edited and adjusted by eye using MEGA7 (Kumar et al., 2016). jModelTest v.2.1.7 (Darriba et al., 2012) was used to determine the evolutionary model, which was then analyzed using Bayesian inference (BI) with MrBayes v3.2.6 (Ronquist et al., 2012) and maximum likelihood (ML) with RAxML v8.2.10 (Stamatakis, 2014). Microphotographs were taken by using an Olympus BX53 (Tokyo, Japan) light microscope with an Olympus DP80 digital camera and cellSens standard image analysis software (Tokyo, Japan).

A sequencing library was prepared using an NEBNext Ultra DNA Library Prep Kit for Illumina (New England Biolabs, United States) and sequenced with an Illumina NovaSeq6000 at Novogene (Beijing, China). The data were trimmed using SOAPnuke v1.3.0 (Chen et al., 2017) and then assembled with SPAdes v3.10.1 (Bankevich et al., 2012). The resulting assembly contigs were determined to be form the chloroplast genome based on the following criteria: (1) blast searches of publicly available chloroplast genomes of Chlorophyta algae species with significant e-values (1e-5); (2) the GC content of the contig is less than 45% (the GC content of green algae chloroplast genomes that have been sequenced to date is normally less than 45%); and (3) the sequencing depth is higher than 100×. Then, the trimmed reads were aligned to the resulting assembly contigs by BWA-MEM v0.7.12 (Li, 2013). If reads mapped two contigs at the same time, we determined the order of contigs, and after confirming the orders of the contigs, the sequence we produced was then rechecked by Sanger dideoxy sequencing technology and synteny analysis with related species. The chloroplast genomes were initially annotated using CpGAVAS (Liu et al., 2012). Protein-coding genes were further polished using Blast with genes from the available colonial volvocine chloroplast genes. All chloroplast genome sequences have been submitted to GenBank, the accession number was listed in Table 1.

Our study aims to reveal the evolutionary relationship between colonial and unicellular species in Chlamydomonadales based on the protein-coding genes of chloroplast genomes, we tried to ensure the gene we analyzed are exist in all species, however, we found some species may have limited number of protein-coding genes and some genes of specific species have poor alignment with other species, to keep a balance between the number of genes and the number of species, 12 colonial species and 16 unicellular species were chosen and species information was listed in Table 1. The data set was assembled from the following 55 protein-coding genes: atpA, atpB, atpE, atpF, atpH, atpI, ccsA, cemA, chlB, chlL, chlN, clpP, petA, petB, petD, petG, petL, psaB, psaC, psaJ, psbA, psbB, psbC, psbD, psbE, psbF, psbH, psbI, psbK, psbL, psbM, psbN, psbT, psbZ, rbcL, rpl14, rpl16, rpl2, rpl20, rpl23, rpl36, rpl5, rps11, rps12, rps14, rps18, rps19, rps2, rps3, rps4, rps7, rps8, rps9, tufa, and ycf4. The method of phylogenomic analysis was mainly refer to Lemieux et al. (2015). All genes were aligned using MUSCLE v3.8 (Edgar, 2004), and the alignments of all genes were converted into a codon alignment by TranslatorX (Abascal et al., 2010). The ambiguously aligned regions in alignment were excluded using Gblocks0.91b (Castresana, 2000) with the options −t = c, −b3 = 5, −b4 = 5, and −b5 = half. All alignments were concatenated using Phyutility v2.2.6 (Smith and Dunn, 2008), and then the Degen1.pl 1.2 script (Regier et al., 2010) was applied to the concatenated alignment. jModeltest v.2.1.7 (Darriba et al., 2012) was used to determine the evolutionary model. The data was partitioned by gene, with the model applied to each partition. Phylogenies were inferred using ML and BI methods. ML analyses were carried out using RAxML v8.2.10 (Stamatakis, 2014) and the GTRGAMMA model of sequence evolution. Bootstrap analysis with 1,000 replicates of the dataset for ML was performed to estimate statistical reliability. Bayesian analyses were performed with MrBayes v3.2.6 (Ronquist et al., 2012) with the GTR + I + G model. Markov chain Monte Carlo (MCMC) analyses were run with four Markov chains (three heated, one cold) for 1,000,000 generations, with trees sampled every 500 generations. Each time the diagnostics were calculated, a fixed proportion of samples (burninfrac = 0.25) were discarded from the beginning of the chain. A stationary distribution was assumed when the average standard deviation of the split frequencies was lower than 0.01.

The CODEML program of PAML v4.9 (Yang, 2007) with the ML model (runmode = −2, CodonFreq = 2) was used to measure the values of dS and dN, the analysis was based on 55 chloroplast protein-coding genes. As Chloromonas radiata belongs, with Carteria cerasiformis and Carteria sp., to a clade that is sister to all the other species, so C. radiata was used as reference. Comparisons of the evolutionary rates were conducted using the two-tailed Wilcoxon rank sum test. The multiple testing was corrected by applying the false discovery rate method (FDR) (Benjamini and Hochberg, 1995) as implemented in R.¹ The phylogenetic tree was used as a constraint tree, but branch lengths were inferred by using PAML.

The ML method is a pairwise approach to estimate the dN/dS ratio, a dN/dS ratio may indicate in one or both species, and some specific sites under positive selection may remain undetected (Dussert et al., 2018). So, two precise assessments were used to detected the difference of dN/dS and positive selection.

We use the branch model to test whether unicellular species have a different dN/dS ratio relative to the colonial species. The unicellular species were labeled as the foreground branch. A null model (model = 0), where one dN/dS ratio was fixed across all species, was compared with an alternative model (model = 2), where the unicellular species was allowed to have a different dN/dS. Likelihood ratio tests (LRT) were used to test model fit and the Chi-square test was applied for testing P values. The multiple testing was corrected by FDR. Genes were considered putative fast-evolving genes if they had an FDR-adjusted P < 0.05 and a higher dN/dS ratio in the foreground branch than in the background branches.

We use the branch-site model to find genes that potentially experienced positive selection. The improved branch-site model (model = 2, Nsites = 2) was used to detect signatures of positive selection on individual codons in a specific branch (Zhang et al., 2005). The unicellular species were set as the foreground branch. The null model assumed no positive selection occurred on the foreground branch (fix_omega = 1, omega = 1), and the alternative model assumed that sites on the foreground branch were under positive selection (fix_omega = 0, omega = 2). LRT were used to test model fit and the Chi-square test was applied for testing P values. We performed a correction for multiple testing using an FDR criterion, and Bayes empirical Bayes (BEB) method was used to statistically identify sites under positive selection. Genes were considered putative selected genes if they had an FDR-adjusted P < 0.05.

The overlap gene represent a group of genes that are not only putative fast-evolving genes but also putative positive selected genes, we analyzed the function of positive selected sites of overlap genes, the functional information were derived from the Uniprot.² The three-dimensional (3D) structures were predicted using Phyre2 (Kelley et al., 2015). The 3D structures were visualized by ePlant Web server (Fucile et al., 2011).

Vegetative colonies of strain FACHB-2326 were ellipsoidal in shape, composed by 16 cells of approximately identical sizes embedded by gelatinous matrix forming a hollow colonial structure. Cells spherical, the chloroplast contained more than two pyrenoids of almost identical size (Figure 1A) and prominent longitudinal striations (Figure 1B). Only two contractile vacuoles distributed near the base of the flagella (Figure 1B). Based on the morphology of colony, the number and size of pyrenoids, the number of contractile vacuoles, we identified strain FACHB-2326 as C. charkowiensis. Vegetative colonies of strain FACHB-2337 were ellipsoidal in shape, cells were embedded in a gelatinous matrix forming a hollow sphere. Cells were more or less contiguous and appeared nearly tetragonal by mutual compression. There were no prominent spaces between adjoining cells. Based on the spaces between adjoining cells and the shape of cells, we identified strain FACHB-2337 as V. compacta (Figures 1C,D). Vegetative colonies of strain FACHB-2361 were ellipsoidal in shape and contained eight cells compactly arranged in a gelatinous matrix. The chloroplast had more than two pyrenoids. Based on the shape of cells and the number of pyrenoids, we identified strain FACHB-2361 as P. colemaniae (Figures 1E,F). Vegetative colonies of strain FACHB-2362 resembled strain FACHB-2361, but the chloroplast contained a single, basal pyrenoid, so we can identify strain FACHB-2362 as P. morum (Figures 1G,H). Vegetative colonies of strain FACHB-2363 resembled strain FACHB-2326, but the chloroplast striations were not prominent and contained a large basal pyrenoid and small pyrenoids, so we identified strain FACHB-2363 as C. angeleri (Figures 1I,J). Vegetative colonies of strain FACHB-2364 resembled strain FACHB-2363, but the chloroplast only contained a single, basal pyrenoid, so we identified strain FACHB-2364 as Y. unicocca (Figures 1K,L).

The phylogenetic tree based on five genes (Figure 2) supported our morphological identification with high bootstrap value (both 100) and Bayesian posterior probability (both 1.00), except the phylogenetic position of strain FACHB-2337. We noticed that strain FACHB-2337 clustered with Volvulina pringsheimii form a lineage. The cells of V. pringsheimii are multiangular or circular and not contiguous in surface view, and the colony of V. pringsheimii is a hollow sphere more resembles with Eudorina. But the colony of strain FACHB-2337 was more compact and resembled with Pandorina, and the cells were contiguous in surface view, the morphological observation strongly supported strain FACHB-2337 as V. compacta rather than V. pringsheimii (Starr, 1962; Nozaki and Kuroiwa, 1990). So, we still considered strain FACHB-2337 as V. compacta. In our phylogenetic tree, the phylogenetic position of most species was consisting with the study of Nozaki et al. (2014), except Volvulina steinii. The bootstrap value and posterior probability of V. steinii clade were relatively low in the study of Nozaki et al. (2014), we found the phylogenetic position of V. steinii in our study was consisting with Nakada et al. (2010), and both studies showed high bootstrap value and posterior probability of V. steinii clade. Meanwhile, recent study both show polyphyletic of genera Pandorina and Volvulina (Coleman, 2001; Nakada et al., 2010; Nozaki et al., 2014), so the phylogenetic position of these species may need further study.

We conducted our phylogenomic analysis based on the nucleotide sequence of 55 chloroplast protein-coding genes, ML was carried out using RAxML, Bayesian analyses was performed with MrBayes, the phylogenetic position of most species inferred from both methods are the same except Dunaliella salina. The phylogenetic position of D. salina inferred from both method have high bootstrap value (87) and posterior probability (1.00), but the result of ML method was in accordance with previous study (Lemieux et al., 2015), so the ML tree was used to represent the result (Figure 3). The phylogenetic position of other unicellular species were consistent with previous study with high bootstrap values and posterior probability values (Yumoto et al., 2013; Lemieux et al., 2015). The phylogenetic position of most colonial species were consistent with previous study (Nozaki et al., 2014) except P. morum, we noticed that the P. morum together with V. compacta formed a lineage instead of P. colemaniae, this situation have been reported before (Coleman, 2001), this may mainly due to the polyphyletic of Pandorina (Herron, 2016). In our study, this tree was used as the constraint tree for our evolutionary analysis.

Based on the ML method of 55 chloroplast protein-coding genes, the value of dN and dS were compared between colonial and unicellular species in Chlamydomonadales (Table 2 and Supplementary Table S1). When refer to dN, 27 genes were significantly different between the two group of algae, among these genes, we found 19 genes significantly higher in unicellular species. When refer to dS, 10 genes were significantly different between the two group of algae, among these genes, we found six genes significantly higher in unicellular species. Among genes with statistical significance, both comparisons show more genes have higher substitution rates in the unicellular species.

The branch model was used to compare the dN/dS ratio between colonial and unicellular species based on the chloroplast protein-coding genes, the LRT was used to compare the fit of two models. The null model (H0) assumed that all tree branches evolved at the same rate (the same dN/dS ratio), the alternative model assumed that the foreground branch (the unicellular species) could evolved at a different rate (different dN/dS ratio). We found the FDR-adjusted P value of 16 genes were less than 0.05, this indicates the dN/dS ratio of these genes were significantly different among the unicellular and colonial species (Figure 4 and Supplementary Table S2). Among these genes, the dN/dS ratio of 14 genes (atpA, rpl16, psaB, psbC, atpB, psbE, psaJ, psbA, rps8, rpl2, rps12, rps14, psbN, atpI) were higher in the unicellular species compared with the colonial species, so these genes were considered putative fast-evolving genes.

The positive selection analysis was performed based on the branch-site model, and we also conducted the comparison between the null and alternative models. The null model considered the foreground branch only have dN/dS = 1, and the alternative model considered sites on the foreground branch have dN/dS > 1 (positive selection). We used the Chi-square test to testing P values, after the FDR correction, we found 11 genes (psaB, psbB, psbC, rbcL, tufA, psbA, rps4, rpl5, rpl16, rps12, atpF) have the FDR-adjusted P value lower than 0.05 (Supplementary Table S3), and we considered these genes as the putative positively selected genes. The overlapped genes between the putative fast-evolving genes and positively selected genes were psaB, psbC, psbA, rpl16, and rps12. Based on the BEB method, the positively selected sites for each gene were shown in Table 3. We found the psaB, psbC and psbA have sites may likely under positive selection. For psaB, 134GLN have posterior probability higher than 95%, 253GLN have posterior probability higher than 90%, but there is no related functional sites information of Chlamydomonas reinhardtii in Uniprot, so the positively selected sites of psaB reminds further study. For psbA, site 237ARG (posterior probability higher than 90%) was close to the 215HIS and 272HIS of C. reinhardtii, the 215HIS is the metal binding site of iron and binding site of Quinone (B), the 272HIS is also the metal binding site of iron. For psbC, site 409SER (posterior probability higher than 95%) was close to the 355GLU of C. reinhardtii which was the metal binding site of calcium-manganese-oxide (Figure 5).