Protein and CDS sequences of E. siliculosus V2016 version were downloaded from the website http://bioinformatics.psb.ugent.be/orcae/overview/Ectsi. Those of S. japonica were downloaded from NCBI under the accession code JXRI00000000. Sequences of C. okamuranus were downloaded from http://marinegenomics.oist.jp/algae/. Other genomes were downloaded from the following sites: green plants, Phytozome V12.0 https://phytozome.jgi.doe.gov/pz/portal.html. Emiliania huxleyi, oomycetes, http://genome.jgi.doe.gov/pages/tree-of-life.jsf. Other genes were downloaded at NCBI https://www.ncbi.nlm.nih.gov/. The LOX domain PF00305 was downloaded from the Pfam website http://pfam.xfam.org/ (Finn et al., 2015). HMMER was used to search for this domain in the predicted proteins of each species. All of the obtained LOXs were further manually examined using the NCBI online BLAST tool https://blast.ncbi.nlm.nih.gov/Blast.cgi to confirm the genes’ annotation. The online InterProScan program http://www.ebi.ac.uk/interpro/search/sequence-search was used to confirm the presence of LOX domain in each sequence.

Sequence alignments were performed on both the full-length proteins and the conserved PF00305 domain using MUSCLE embedded in MEGA 7.0 (Kumar et al., 2016). The maximum likelihood (ML) and Neighbor-joining (NJ) phylogenetic trees were constructed with best substitution models and rate parameter WAG + G + I found by MEGA 7.0. Bootstrap with 1000 bootstrap replicates was performed to obtain the confidence support level. Bayesian inference of phylogeny was carried out using BEAST2 (Bouckaert et al., 2014). Analyses were run for 100 million generations, with sampling every 10000th generation using the same substitution model as above and Yule Model as a prior. The first 10000 samples were discarded as burn-in.

Intron and exon information of the three brown algal LOX genes was obtained from their GFF files and was used for the graphic display using the Gene Structure Display Server of Peking University (Hu et al., 2015). The protein transmembrane helices were predicted by CCTOP¹ (Dobson et al., 2015). Structural motifs were searched for using the MEME program http://meme-suite.org/ (Bailey et al., 2009). The motif number was set to 15. The most conserved motif of LOXs was displayed on the WebLogo website http://weblogo.berkeley.edu/logo.cgi (Crooks et al., 2004). Signal peptides were predicted using the website http://www.cbs.dtu.dk/services/SignalP/ (Petersen et al., 2011). Subcellular protein localization was predicted using ASAFind (Gruber et al., 2015). Mass weight (MW) and isoelectric point (pI) was calculated using the website http://www.expasy.org/.

In total, 18 protein sequences of LOXs in three brown algae Ectocarpus, S. japonica, C. okamuranus were aligned using MUSCLE. A ML tree was constructed with MEGA 7.0 as described above. Brown algae LOXs were divided into two clades C1 and C2. In each clade, the branches were divided again into Ectocarpales clades (C-E) and Laminariales clades (C-S). To evaluate variation in selective pressures between different clades (C1-E vs. C1-S; C2-E vs. C2-S), a branch model in PAML v4.8 (Yang, 2007) was used to estimate the ω (dN/dS) ratio under two assumptions: a one ratio model (M = 0) that assumes the same ω in C1 and C2 clade, respectively; a two-ratio model (M = 2) that assumes that Ectocarpales clade (C-E) and Saccharina clade (C-S) have different ω ratios. To test which of the models best fit the data, likelihood ratio tests (LRTs) were performed by comparing twice the difference in log-likelihood values using a χ2 distribution with one degree of freedom. We also ran site models (M = 0, NSsites = M0, M3, M7, M8) in Codeml to test if there were sites of positive selection in brown algal LOXs. LRTs were performed (M0 vs. M3, M7 vs. M8) to verify which model fit the data using a χ2 distribution with two degrees of freedom. To further explore the evolutionary dynamics of expanded LOXs in S. japonica, pairwise dN, dS and dN/dS values between duplicated pairs of S. japonica LOXs were calculated using the online PAL2NAL program http://www.bork.embl.de/pal2nal/index.cgi (Suyama et al., 2006).

Previous microarray data of the E. siliculosus transcriptome (Dittami et al., 2009; Ritter et al., 2014) was used to explore the expression levels of LOXs in response to abiotic stressors, notably hyposaline stress, hypersaline stress, oxidative stress, and Cu stress. The expression level was determined by averaging the expression values (previously quantile normalized by Roche NimbleGen, Madison, WI, United States) of four replicates for each experimental condition. Furthermore, expression levels of LOXs were examined in different life cycle stages. In E. siliculosus we compared previously published RNA-Seq data of parthenosporphytes, male gametophytes, and female gametophytes (Lipinska et al., 2015). Expression levels of S. japonica LOXs were examined in sporophytes, male gametophytes, and female gametophytes. Their transcriptomes were sequenced by our group using an Illumina HiSeq 2500, with three biological replicates in each life stage. Gene expression levels were estimated by FPKM (Fragments Per Kilobase of transcript per Million fragments mapped). Differential expression analysis among life stages was performed using the DESeq (Love et al., 2014) R package using a model based on the negative binominal distribution. The resulting p-values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate (Benjamini and Hochberg, 1995). Genes with an adjusted p-value < 0.05 found by DESeq were assigned as differentially expressed. The raw sequence data was deposited at NCBI database. The SRA accession numbers are SRR5860560, SRR5860561, SRR5860562, SRR5860563, SRR5860564, SRR5860565, SRR5860566, SRR5860567, and SRR5860568.

To predict the 3D structures of S. japonica LOX proteins, the X-ray structure of a P. aeruginosa LOX protein (PDB accession numbers 4g33.1, 5lc8.1, 4rpe.1) was used as template, as it has the highest coverage and similarity with S. japonica LOXs. Each LOX was aligned with the template and the structure was automatically built using SWISS-MODEL https://www.swissmodel.expasy.org/ (Biasini et al., 2014). Secondary structure was displayed using the Espript web server http://espript.ibcp.fr/ESPript/cgi-bin/ESPript.cgi (Gouet et al., 2003).

A comprehensive genome-wide search identified 12, 4, and 3 LOX genes in genomes of the S. japonica, E. siliculosus, and C. okamuranus, respectively. The number of LOX genes in S. japonica is approximately three-fold that in E. siliculosus and C. okamuranus. The accuracy of these gene models was manually checked by comparison with publicly available RNA-Seq data. We found a large number of reads spanning the SJ21862 and SJ21863, suggesting they correspond to a single transcript encoding a 518 residue protein. Furthermore, SJ16633 and SJ16634 were truncated at the N-terminal end and manually extended increasing their length by 25 and 91 residues, respectively. Additionally, SJ08291 and SJ05859 have the shortest protein length, and no expression data available, suggesting they are probably pseudogenes. Detailed information about the LOX family in brown algae is listed in Table 1. The length of these proteins ranged from 346 to 747 amino acid residues. The molecular masses ranged from 39 to 82 kDa, with SJ08291 and Ec-03_000010 being the smallest and largest, respectively. The pI values of these proteins are similar, ranging from 4.65 to 5.79. In S. japonica, E. siliculosus, and C. okamuranus, three, one, and two LOXs were predicted to be located in the chloroplast, respectively. In order to perform phylogenetic analysis, we also searched for the LOX genes in other organisms, including oomycetes, Phaeodactylum tricornutum, Thalassiosira pseudonana, Fragilariopsis cylindrus, Nannochloropsis gaditana, and Emiliania huxleyi, all of which, together with brown algae belong to the Stramenopiles and Haptophytes. Besides, LOX in prokaryotes, model plants and animal were included. Sequences corresponding to the LOX domain (PF00305) were aligned using MUSCLE imbedded in MEGA 7.0. Then ML and NJ trees were constructed using MEGA 7.0, and a Bayesian tree was constructed using BEAST 2. The three trees exhibited nearly the same topology. As shown in Figure 1, brown algal LOXs formed one robust clade, with 100% bootstrap support. It is interesting that LOXs of brown algae are most closely related to those of prokaryotes with high bootstrap support but not with other Stramenopiles, plants, or animals. From this tree topology, we propose that LOXs in both brown algae and prokaryotes, as well as LOXs in plant/animal or haptophytes/oomycetes originate from several different ancestral LOXs, and that each lineage retained different copies.

To gain further insight into the diversity of brown algal LOXs, we constructed a separate phylogenetic tree using only the brown algal LOX protein sequences. The tree showed that the 18 LOXs were divided into two distinct clades (C1 and C2, Figure 2). C1 and C2 contained 9 genes each. The C1 clade was further divided into two subclades (C1-S and C1-E). Similarly, the C2 clade also contained two subclades (C2-S and C2-E), in which the ‘S’ and ‘E’ denote Saccharina and Ectocarpales, respectively. Subclades C1-E and C2-E contained LOXs from E. siliculosus and C. okamuranus, whereas subclasses C1-S and C2-S contained only LOXs from S. japonica. Striking gene structure conservation was found within 18 LOXs of brown algae in terms of either intron numbers or exon length. A total of 15 out of 18 genes contain from 15 to 20 introns, of which ten genes contain 17 introns with exons of similar length. Fewest exons were seen in SJ05859 and SJ08291, with only ten exons, which further supported that they are possibly pseudogenes. Intron length in S. japonica varies from 39 bp in SJ16634 to 9,722 bp in SJ12567 whereas that in E. siliculosus varies from 159 bp in Ec-03_000010 to 6,520 bp in Ec-20_003610. Pairwise comparison of the 18 LOX full-length protein sequences revealed some features (Figure 3). C1 clade proteins showed 25.4–85.5% pairwise sequence identity. C2 clade showed 14.1–78.1% pairwise sequence identity. However, the protein sequences between the two LOX clades were significantly more different (Mann–Whitney test, p < 0.01) and the pairwise sequence identity was 8.2–58.2% (Supplementary Table S1).

To explore the detailed motif composition of the brown algal LOX family, we predicted putative motifs using the online tool MEME (Bailey et al., 2009). In total, 15 motifs were identified (Figure 4 and Supplementary Table S2). The presence of conserved motifs in LOX proteins provides additional support for the results of the phylogenetic analysis. Compared with the C2 clade, two distinct motifs occurred only in the C1 clade, that is, motifs 14 and 15 in the N-terminal region of the protein. However, whether these motifs confer a specific function to members of the C1 clade remains to be surveyed. Motifs 1 and 6 contain the representative 38 amino acid residues with the structure His-(X)4-His-(X)4-His-(X)17-His-(X)8-His. In the protein sequence alignment of S. japonica, all of the nine sequences have the conserved motif in this His rich region, except two sequences (SJ03316 and SJ16634) where the first His is replaced by Tyr (Figure 5). Motifs 4 and 5 contain an Asn and Ile residue respectively, which are essential for the binding of iron. Notably, some LOXs lack critical motifs. For instance, SJ05859 lacks the six motifs in C-terminal whereas SJ08291 lacks the eight motifs in N-terminal. One gene in C. okamuranus lacks motif 1 and four motifs in C-terminal.

The ratio of non-synonymous (dN) to synonymous (dS) nucleotide substitutions ω (dN/dS) provides information about the evolutionary forces operating on a gene or gene region (Yang et al., 1998). ω < 1 suggests negative selection and the removal of adverse mutation. ω = 1 and ω > 1 suggest neutral selection and positive selection, respectively. The LOX family of brown algae was divided into two clades (C1 and C2), and they were further divided into subclades (C1-S and C1-E, C2-S and C2-E). To infer the selection pressure within each clade, we estimated ω values with branch models using the CODEML program in the PAML software package (Table 2) (see Materials and Methods). The log-likelihood values under the one ratio model were lnL = -8197.2974 and -7447.3493, with estimates of ω = 0.10530 and 0.1262 for the C1 and C2 clades LOX, respectively. Under the two-ratio model, the log-likelihood values were lnL = -8193.6592 and -7441.3471 for the C1 and C2 clade LOX genes, respectively. Likelihood ratio tests (LRTs) of the two assumptions indicated that the two-ratio model was significantly more likely (p < 0.01) for both the C1 and C2 LOX clades, suggesting that selective pressure varied among LOX subclades (C1-S vs. C1-E, C2-S vs. C2-E). LOX genes in S. japonica have experienced more relaxed purifying selection than those in E. siliculosus and C. okamuranus. All of the estimated ω values were <1, suggesting that the two LOX clades have been under purifying selection but that the selection constraint on the C2 clade was more relaxed than that on the C1 clade, which is consistent with the lower sequence identity in the C2 clade. However, the ratio averaged over all sites is unlikely to exceed 1, because positive selection is not likely to affect all codon sites of a gene (Yang, 2007). Site model allows the ω ratio to vary among sites over all branches (Roux et al., 2014). Therefore, we employed site models (M0 vs. M3, M7 vs. M8) to test whether episodic positive selection had affected specific amino acid sites in each gene. M0 assumes one same ratio among all sites whereas M3 assumes discrete ratios among sites. Under the M7 model, ω follows a beta distribution between 0 and 1, while no ω > 1 site allowed. The M8 model allows an additional site class with ω > 1 (Yang, 2007). LRT was performed to test which model fits the data best. The site model analysis between M0 and M3 indicated that selection pressure varied among sites both for C1 and C2 clades (p < 0.001). However, there is no significant difference between the lnL values of M7 and M8, indicating that no site was positively selected, which further supported the results from branch model.

In order to illustrate the expansion of LOXs in S. japonica, we compared the LOX number across different organisms. Proteome size-dependent expansion of metal-binding proteins was observed across all domains of life (Dupont et al., 2010). Thus, the relationship between LOX number and proteome size was taken into consideration (Figure 6). LOXs are rare among most sequenced algae. The green alga Volvox carteri has four members, while other green algae, red algae, and diatoms have maximally one LOX. The two brown algae E. siliculosus and C. okamuranus have four and three LOXs, respectively. Notably, S. japonica, which possesses 11 LOXs, has expanded its LOX family compared with other algae. Furthermore, this expansion goes beyond the mere effect of proteome size, as indicated by the different slopes for brown algae and other algae in the linear regressions in Figure 6A. Moreover, a similar effect was not observed for other metal-binding protein families such as the cytochromes P450 (CYP) family (Figure 6B).

According to the physical locations of the LOXs on chromosomes or scaffolds, gene duplication events were constructed. In E. siliculosus, the four LOX genes were assigned to two chromosomes. Within each chromosome, tandem duplication occurred once. However, LOX genes in S. japonica have undergone more rounds of duplication and formed a more expanded LOX gene family compared to E. siliculosus. Tandem duplication plays a major role in the S. japonica LOX expansion history. Three tandem clusters account for 7 out of 11 LOX genes, including SJ05858/SJ05859, SJ21859/SJ21862, and SJ16632/SJ16633/SJ16634. Other genes likely came from segmental duplication after a series of tandem duplications. The dN, dS and dN/dS ratio of these putative paralogous gene pairs were calculated to investigate the divergence fate after duplication (Table 3). All dN/dS ratios were <0.5, suggesting strong purifying selection, which is consistent with the branch/site model results. The largest ratio was observed in the gene pair SJ05859/SJ12567 with the value of 0.4393, indicating the relatively relaxed selection pressure between them. Duplicate paralogs diverge from the ancestral state and accumulate synonymous mutations at similar rates over time (Holland et al., 2017), so the dS value can be used to measure the time since gene divergence. The lowest dS (0.0520) was observed between Ec-20_003610/ Ec-20_003620, suggesting that they were created most recently, whereas the highest dS (1.6936) was found between SJ21859/SJ21862, suggesting an ancient duplication, likely derived from the first round of ancestral tandem duplications. By reconciliation of the physical position of genes, gene tree, and the dS values, we attempted to reconstruct the most parsimonious history of S. japonica LOX gene duplication (Figure 7). For the C1 clade, it seems that one tandem duplication created two ancestral genes (b and c), after which other events, likely segmental duplications, occurred to form two gene pairs; one pair is SJ21859 and SJ21862, the other pair is ancestral d and e. After that, two duplication events occurred sequentially, to create two more pairs of LOXs. On the other hand, in the C2 clade, ancestral gene A duplicated to form SJ08291 and ancestral B, followed by two rounds of tandem duplications to create SJ16632 and ancestral C, D, which further produced SJ16633, SJ16634, and SJ03316.

Duplicated genes may evolve into different fates, which can be inferred by the divergence of their expression patterns. The expression patterns of LOX gene family members were examined using the publically available microarray data of E. siliculosus and RNA-sequencing data of E. siliculosus and S. japonica. Three of the four LOX genes in E. siliculosus are present in the microarray data, and all the three genes have two contigs/singletons. Hierarchical clustering showed that expression levels under Cu stress were clustered with oxidative stress (Figure 8A), which is consistent with the meta-analysis in Ritter et al. (2014). One gene (Ec-20_003610) was significantly influenced by hyposaline, hypersaline, and oxidative stress (log2-fold change > 1 compared to control). The other two genes were up-regulated by hyposaline, hypersaline, and Cu stress, though the log2-fold change did not reach one. This result indicates that LOX genes in E. siliculosus are responsive to these stress conditions, and may explain the observed accumulation of oxylipins under Cu stress (Ritter et al., 2014). Expression divergence was also observed in different life stages. In E. siliculosus, mature male and female gametophytes have significantly higher LOX mRNA levels than sporophytes or immature gametophytes. Notably, all of the four LOXs were male-biased in immature gametophytes (fold change > 2, adjusted p-value < 0.1; Figure 8B), suggesting their potential roles in sexual dimorphism. For the LOX members in S. japonica, more variation in expression patterns was found in different life stages (Figure 8C and Table 3). Three genes were sex-biased, including one female-biased (SJ16633) and two male-biased (SJ16632 and SJ12567). Two genes in the C1 clade (SJ05858, SJ21859) and two genes in the C2 clade (SJ16632, SJ16633) were ranked as the most highly expressed in sporophytes. However, expression levels of other genes were relatively low, especially in gametophytes. Two genes (SJ21862, SJ03316) were exclusively expressed in sporophytes, whereas no expression was detected in gametophytes, suggesting that they are likely only functional in sporophytes. The other genes also exhibited higher expression levels in sporophytes than in gametophytes, except for one gene (SJ12567), which showed the highest expression in male gametophytes. Another gene (SJ08291) was not expressed in any life stage, and SJ05859 had very a low FPKM of only 0.08. These two genes are likely to be pseudogenes, considering their gene structure and expression level. Notably, the two pairs of LOX genes with lowest expression levels (SJ05859 vs. SJ12567, SJ03316 vs. SJ16634) have higher dN/dS values than other LOX gene pairs (Table 3). Likewise, in E. siliculosus, the lowest overall expression level (Ec-20_003610) and great expression difference was observed in the gene pair with the highest dN/dS value (0.3843 compared to 0.1577 in the other pair of genes). The expression patterns of LOX suggested that functional diversification has occurred between the duplicated LOX gene pairs.

Three-dimensional protein structure modeling across all eleven LOXs in S. japonica revealed that, in general, these proteins are best modeled by a 15-LOX in the bacterium P. aeruginosa (coverage 79–89%, identity 27.83–41.51%), which further supports the close relationships between them. The simulated structures of S. japonica LOX showed remarkable differences to plant or animal LOXs. For most plant and animal LOXs, a two-domain structure has been identified, but LOXs in brown algae fold into a single domain. The helical lid structure in the P. aeruginosa LOX was suggested to be the functional equivalent of the N-terminal β-barrel in eukaryotic LOXs (Garreta et al., 2013). In LOXs of S. japonica, the N-terminal part seems to also contain the helical lid structure (α2A, α2B), though having shorter helices than the P. aeruginosa LOX (Figure 9 and Supplementary Figure S1). The X-ray structure of 4g33.1 showed that five sites, His377, His382, His555, Asn559, and Ile685 are essential for binding of iron (Garreta et al., 2013; Kalms et al., 2017). Sequence comparison showed that the five sites were conserved in all the S. japonica LOXs, except for protein SJ05859 due to the short sequence length. The predicted 3D structure of SJ16632 showed that the five sites could interact with the iron. The simulated structures of six LOXs in the C1 clade and five LOXs in the C2 clade were superimposed to evaluate the goodness of fit of the overall topologies (Supplementary Figure S2). The result shows that the structure of the α-helix is conserved across all members, whereas structures in the loop regions are variable and less conserved, which might be due to relaxed functional constraints.