Transcription Profiles of Age-at-Maturity-Associated Genes Suggest Cell Fate Commitment Regulation as a Key Factor in the Atlantic Salmon Maturation Process

Despite recent taxonomic diversification in studies linking genotype with phenotype, follow-up studies aimed at understanding the molecular processes of such genotype-phenotype associations remain rare. The age at which an individual reaches sexual maturity is an important fitness trait in many wild species. However, the molecular mechanisms regulating maturation timing processes remain obscure. A recent genome-wide association study in Atlantic salmon (Salmo salar) identified large-effect age-at-maturity-associated chromosomal regions including genes vgll3, akap11 and six6, which have roles in adipogenesis, spermatogenesis and the hypothalamic-pituitary-gonadal (HPG) axis, respectively. Here, we determine expression patterns of these genes during salmon development and their potential molecular partners and pathways. Using Nanostring transcription profiling technology, we show development- and tissue-specific mRNA expression patterns for vgll3, akap11 and six6. Correlated expression levels of vgll3 and akap11, which have adjacent chromosomal location, suggests they may have shared regulation. Further, vgll3 correlating with arhgap6 and yap1, and akap11 with lats1 and yap1 suggests that Vgll3 and Akap11 take part in actin cytoskeleton regulation. Tissue-specific expression results indicate that vgll3 and akap11 paralogs have sex-dependent expression patterns in gonads. Moreover, six6 correlating with slc38a6 and rtn1, and Hippo signaling genes suggests that Six6 could have a broader role in the HPG neuroendrocrine and cell fate commitment regulation, respectively. We conclude that Vgll3, Akap11 and Six6 may influence Atlantic salmon maturation timing via affecting adipogenesis and gametogenesis by regulating cell fate commitment and the HPG axis. These results may help to unravel general molecular mechanisms behind maturation.

(vestigial-like family member 3) explained as much as 39% of the phenotypic variation in maturation age in both males and females (Barson et al. 2015). Other studies of European Atlantic salmon have also observed associations between maturation and the same genome region (Ayllon et al. 2015;Ayllon et al. 2019;Czorlich et al. 2018; M. Sinclair-Waters, J. Ødegård, S. A. Korsvoll, T. Moen, S. Lien, C. R. Primmer and N. J. Barson, unpublished results). However, associations in North American-derived salmon populations and aquaculture stocks have been mixed, possibly due to little or no polymorphism at the locus in North American populations (Boulding et al. 2019;Kusche et al. 2017;Mohamed et al. 2019). In addition to Atlantic salmon, VGLL3 has also been linked with pubertal timing, growth and body condition in humans (Elks et al. 2010;Cousminer et al. 2013;Tu et al. 2015), which indicates that it may have an evolutionarily conserved role in the regulation of vertebrate maturation timing. Age-at-maturity is generally a polygenic trait, controlled by several small-effect loci (Elks et al. 2010;Cousminer et al. 2013;Perry et al. 2014;Day et al. 2017;Zhu et al. 2018a), and, thereby, the identification of a large-effect locus in salmon provides a rare opportunity to investigate the molecular processes behind this association.
Sexual maturation is a biological process stemming from a complex chain of events culminating in the first reproduction. The maturation process commences already in the embryo after fertilization by allocating energy to the growth and differentiation of developing gonads and is completed when gametes are produced (Laird et al. 1978;Okuzawa 2002;Thorpe 2007). Although timing of maturation is known to be mediated by interplay between fat accumulation and activation of the hypothalamic-pituitary-gonadal (HPG) axis (Kaplowitz 2008;Sam and Dhillo 2010;Taranger et al. 2010) the exact molecular mechanisms regulating the process are still obscure.
Maturation requires sufficient fat storage to provide energy for proper gonad development. Therefore, evidence showing that Vgll3 encodes a negative regulator of adipocyte maturation and that its mRNA expression inversely correlates with total body weight and fat content in mice (Halperin et al. 2013) suggests vgll3 is a good candidate for having a role in sexual maturation in salmon. A recent study linking VGLL3 with reduced adiposity indices in the Mongolian human population provides further evidence for general, specieswide role of VGLL3 in adipose regulation (Nakayama et al. 2017). Beyond regulating adipocyte differentiation, VGLL3 has also been shown to have a broader role in mesenchymal-derived cell fate decision. Studies show that Vgll3 overexpression promotes expression of the chondrocyte and osteocyte inducing markers in murine preadipocyte cell line (Halperin et al. 2013) and myogenesis in mouse and human myoblasts (Figeac et al. 2019). Expression pattern of vgll3 during embryonic development (Faucheux et al. 2010;Simon et al. 2016;Simon et al. 2017) and in adult vertebrates (Mielcarek et al. 2009;Faucheux et al. 2010;Kjaerner-Semb et al. 2018;Figeac et al. 2019) in various tissues suggests a broad role for Vgll3 in development. Detection of vgll3 expression in testis (Faucheux et al. 2010;McDowell et al. 2012;Kjaerner-Semb et al. 2018) and ovary (Gambaro et al. 2013;Kjaerner-Semb et al. 2018) further supports the participation of Vgll3 in sexual maturation. The exact molecular mechanisms via which VGLL3 operates on cell fate determination, and also maturation, are unclear, but it is known to be a cofactor for all known TEAD transcription factors (Simon et al. 2017;Figeac et al. 2019). By binding to TEADs, VGLL3 has been shown to influence the Hippo signaling pathway (Figeac et al. 2019) that regulates cell fate commitment and organ growth (Huang et al. 2005;Meng et al. 2016).
In addition to vgll3, two other genes, akap11 (on chromosome 25) and six6 (on chromosome 9), associate with age-at-maturity in Atlantic salmon (Barson et al. 2015). However, association of six6 with maturation timing is only seen before population structure correction. In addition to salmon, SIX6 (SIX homeobox 6) associates with ageat-menarche and adult height in humans (Perry et al. 2014) and puberty in cattle (Cánovas et al. 2014). Six6 encodes a transcription factor whose expression has been widely studied in several vertebrates and detected in the hypothalamus, pituitary gland and testis, organs of the HPG axis (López-Ríos et al. 1999;Jean et al. 1999;Li et al. 2002;Aijaz et al. 2005;Xie et al. 2015). Accordingly, studies in mice show that SIX6 is required for fertility by regulating the maturation of gonadotropinreleasing hormone (GnRH) neurons and expression of GnRH in the hypothalamus (Larder et al. 2011) and repressing transcription of the beta-subunit genes of gonadotropin hormones in pituitary gonadotrope cell line (Xie et al. 2015). In addition to being an important regulator of the HPG axis, SIX6 has an essential role in eye development (López-Ríos et al. 1999;Jean et al. 1999;Li et al. 2002;Aijaz et al. 2005), for example controlling photoreceptor differentiation (Conte et al. 2010), which is crucial for photoperiod sensing in seasonal breeders such as salmon (Melo et al. 2014). Contrary to SIX6 whose role in testis remains so far unknown, AKAP11, although expressed in most studied tissues and cell types, has the strongest signal and suggested functional role in testis and sperm (Reinton et al. 2000). It encodes A-kinase anchor protein 11, which interacts with type I and II regulatory subunits of protein kinase A (PKA) in order to tether PKA to discrete locations within a cell to control many essential functions, such as cell cycle (Grieco et al. 1996) and lipid metabolism (Lee et al. 2016). Expression of AKAP11 mRNA and protein is observed in sperm throughout spermatogenesis and its suggested association with cytoskeletal structure indicates its importance in the sperm function (Reinton et al. 2000), and thus maturation.
Although there is information available about the expression patterns and molecular functions of the age-at-maturity-associated genes vgll3, akap11 and six6 in some vertebrate species, studies covering a range of developmental time points are lacking. Therefore, in order to be able to determine functional molecular mechanisms of these genes during maturation, first, it is crucial to know when and where these genes are expressed, and what are the potential molecular partners and pathways they link with. Since Atlantic salmon has an exceptionally large-effect locus associating with age-at-maturity, it provides an attractive model for studying molecular mechanisms linked with sexual maturation. Thus, our aim was to investigate the expression patterns of vgll3, akap11 and six6, and additional known key genes related to their suggested functions and pathways in a range of Atlantic salmon developmental time points and tissues. Using Nanostring technology, we characterized expression profiles of vgll3, akap11 and six6 paralogs and 205 other genes related to their functions in five Atlantic salmon developmental stages and 15 tissues. Based on our results, we hypothesize a novel role for Vgll3 in participating in cell fate control via actin cytoskeleton regulation, and Akap11 assisting in this process. In addition, we suggest that Six6 may associate broadly with neuroendocrine secretion regulation in the HPG axis, and have a direct role in testis function.

Animals
The Atlantic salmon (Salmo salar) used in this study were derived from hatchery-maintained Neva river strain at the Natural Resources Institute Finland (62°24950$N, 025°57915$E, Laukaa, Finland) and the Inarijoki river, headwater river of the Teno river, near Angeli in northern Finland. Three-year-old wild parr (freshwater stage individuals) from the Inarijoki were caught by electrofishing in early September 2016. The hatchery-maintained Neva strain eggs were fertilized in November 2016 and incubated in circulated water system until hatching. After hatching, alevin (hatched individuals) and fry (individuals able to feed) were grown in standard commercial fish farming conditions. After 14 days from hatching, some alevin were transferred to be grown in hatchery at the Lammi Biological Station (61°0 4945$N, 025°00940$E, Lammi, Finland). Relative ages of the embryos, alevin and fry were calculated in degree days (°d) (number of days after fertilization multiplied by temperature in°C), which are used throughout the text for the embryonic stages. In addition, t s units (time it takes an embryo to form one somite pair in certain temperature in°C) were calculated according to Gorodilov (1996). Three-year-old hatchery-maintained Neva river parr were reared in standard commercial fish farming conditions at the Natural Resources Institute Finland hatchery in Laukaa, Finland.

Sample collections
Pre-hatched embryos from two developmental time points (186°d/177 t s and 258°d/260 t s ) were dissected from eggs, and alevin (436°d/377 t s ) and fry (1119°d, 1320°d and 2212°d) were caught by a plastic Pasteur pipette and net, respectively. From the embryos and alevin, yolk sac was excised with a scalpel. All samples were stored in RNA preservation buffer (25 mM sodium citrate, 10 mM EDTA, 70 g ammonium sulfate/100 ml solution, pH 5.2) at -20°, whereby fry were first killed by anesthetic overdose of Tricaine methanesulfonate. Of these, four 186°d embryos, four 258°d embryos, four alevin, two normal diet fed 1119°d fry, all from Natural Resources Institute Finland, and four 1320°d fry, of which two were fed with normal and two with low-fat diet, from Lammi Biological Station, were chosen for further analysis. In addition, blood samples were extracted from the caudal vein of two normal diet fed and two low-fat diet fed 2212°d fry from Lammi Biological Station and stored in RNAprotect Animal Blood Tubes (Qiagen, Hilden, Germany) at -20°. Samples from adipose, brain, eye, fin (adipose and caudal), gill, gonad, heart, kidney, liver, muscle, pyloric caeca, spleen and skin tissues were dissected from four hatchery-maintained Neva river parr (two males and two females) and four wild Inarijoki parr (two immature males and two mature males) and stored in RNA preservation buffer at -20°. As not all of the above-mentioned tissues were collected from all eight fish, the specific tissue samples assessed for each individual in this study are described in Supplemental Material, Table S1.

RNA extraction
Altogether 96 samples, including whole embryos, alevin and fry, and tissue samples from parr, were disrupted and homogenized in the presence of NucleoZOL (Macherey-Nagel, Düren, Germany) for lysis in 2 ml safe-lock tubes containing one 5 mm stainless steel bead (Qiagen) for 2.5-3 min at 30 Hz using TissueLyzer II (Qiagen). Total RNA was extracted from the lysate using NucleoSpin RNA Set for NucleoZOL (Macherey-Nagel) according to the manufacturer's instructions. After the elution step, RNA was treated with rDNase using the rDNase Set (Macherey-Nagel) to remove any residual DNA and, subsequently, purified with either NucleoSpin RNA clean-up or NucleoSpin RNA clean-up XS kit (Macherey-Nagel) according to the RNA yield. Blood samples were lysed in the RNAprotect Animal Blood Tubes (Qiagen) and total RNA was extracted using the RNeasy Protect Animal Blood System kit (Qiagen) according to the manufacturer's protocol followed by concentration of RNA using NucleoSpin RNA clean-up XS kit (Macherey-Nagel). The quantity and quality of RNA was determined using both NanoDrop Spectrophotometer ND-1000 (Thermo Scientific, Wilmington, DE, USA) and 2100 Bio-Analyzer system (Agilent Technologies, Santa Clara, CA, USA).

Nanostring nCounter mRNA expression panel
A total of 220 genes were chosen for analysis based on information from the literature, the IPA (Ingenuity Pathway Analysis) tool (Qiagen) and other freely available web tools and databases. These included the age-atmaturity-associated genes vgll3a and akap11a on chromosome 25 and six6a on chromosome 9, and their corresponding paralogs vgll3b and akap11b on chromosome 21 and six6b on chromosome 1, as well as 205 genes potentially linked to these age-at-maturity-associated genes based on their functions and pathways. Further, nine commonly used potential housekeeping genes (actb, ef1aa, ef1ab, ef1ac, gapdh, hprt1, rpabc2a, rpabc2b and rps20) were included in the gene panel as candidates for data normalization. Because of the duplicated Atlantic salmon genome, most genes possess one or more paralogs in their counterpart chromosomes. Therefore, all paralogs of each gene of interest were searched using the SalmoBase (http://salmobase.org/) and NCBI RefSeq databases and included in this study. The exception was that for the genes surrounding the age-at-maturity-associated genes in the chromosomes 25 and 9 according to Barson et al. (2015) the possible paralogs were excluded, as those particular genes were only of interest. Paralogs were renamed by adding a letter at the end of their name alphabetically in order to separate them as shown above. Gene accession numbers, symbols, full names and functional categories are listed in Table S2. Genes were grouped into three different functional categories, "Cell fate", "Metabolism" and "HPG axis", based on their known functions and pathways. Many of the genes have several functions and could thus be placed in multiple categories whereby the chosen categories were the most relevant to the current study. In addition, genes surrounding vgll3a and akap11a, and six6a were included in the "Chr 25 candidate region" and "Chr 9 candidate region" categories, respectively. mRNA expression levels of chosen genes were studied using Nanostring nCounter Analysis technology (NanoString Technologies, Seattle, WA, USA). Probes for each gene paralog, targeted at all known transcript variants, were designed using reference sequences in the NCBI RefSeq database. For some paralogs, it was impossible to design specific probes, as sequence similarity between paralogs was too high. Altogether, 96 samples were analyzed using nCounter Custom CodeSet for probes targeting 220 genes and nCounter Master kit (NanoString Technologies). First, 100 ng of the RNA of each sample was denatured, after which the probes were hybridized with the RNA overnight. Post-hybridization purification and image scanning was performed the following day.

Data analysis
Six genes (ef1ac, prkar1aa, rps20, rxrbaa, vdraa and vdrab) (Table S2) were selected from the 220 genes in the gene panel for use in data normalization as they exhibited a low coefficient of variation and average count values across samples. These genes included two (ef1ac and rps20) of the nine genes originally included as potential normalization candidates, as well as four additional genes (prkar1aa, rxrbaa, vdraa and vdrab). Raw count data from the Nanostring nCounter mRNA expression analysis was normalized by RNA content normalization factor for each sample calculated from geometric mean count values of these six genes. After normalization and quality control, two samples (shown in the Table S1) were removed from the data due to too high RNA content normalization factor values (. 3.0). Quality control and normalization of the data were performed using the nSolver Analysis Software (v3.0) (NanoString Technologies). Normalized mean count values of genes of interest were calculated for all four early developmental stages (embryos, alevin and fry) and 15 parr tissues. For fin, mean count value was averaged across adipose and tail fin samples (no difference was apparent, see results). A normalized count value of 20 was set as a background signal threshold. Forty-two genes were on average below background signal across all early developmental stages, which left 178 genes to estimate 507 pairwise correlations (plus one extra, see below) with the age-at-maturity-associated genes. To estimate 507 pairwise correlations between the three age-atmaturity-associated genes and 178 of the 211 studied genes with average expression above the threshold across the early developmental stages, correlation coefficients were determined among residuals of normalized count data that were controlled for different means of genes and developmental stages, using bivariate linear models under residual maximum likelihood. Correlation standard errors were approximated using a Taylor series. To determine significance, models with unstructured covariance vs. diagonal residual covariance structure were compared using likelihood ratio tests (LRT) (Stram and Lee 1994). As this was an exploratory study, P , 0.01 was considered relevant in LRTs. Genes with expression levels correlating with the age-at-maturityassociated genes were included in analyzing temporal transcript variation during early development. Specifically, a linear mixed model was fitted under residual maximum likelihood for normalized count data with a random sample term to account for technical among-sample variation and with fixed effects for 24 retained genes (see results), four stages and their interaction. The genes exhibited different variances (LRT between model with homovs. heteroscedastic residual variance; Χ 2 23 = 830.7, P , 0.001) and residual variance was, therefore, allowed to be heteroscedastic. Model terms were tested using Wald's F-test and gene-specific pairwise comparisons among the four early developmental stages were conducted, whereby p-values were adjusted for the false discovery rate (FDR) according to Benjamini and Hochberg (1995) and comparisons with FDR , 0.05 were regarded as relevant. Linear models were estimated using ASReml-R (Butler et al. 2009) and data analysis was performed using the R statistical software. The custom R script used to analyze the data is shown in File S1.

Data availability
Raw and normalized mRNA expression data for genes reported in the study are deposited in the NCBI Gene Expression Omnibus (GEO) repository (accession number GSE140434).
Supplementary figure, tables and file have been uploaded to figshare.
Supplementary figure 1: Figure S1. mRNA expression profiles of six6a-related genes in fifteen Atlantic salmon tissues.
Supplementary table 1: Table S1. Tissues of eight three-year-old parr chosen for the Atlantic salmon mRNA expression study.
Supplementary table 2: Table S2. Genes included in the Atlantic salmon mRNA expression study.

RESULTS
Vgll3, akap11 and six6 expression during early developmental stages mRNA expression levels of age-at-maturity-associated genes and their paralogs were studied at four early Atlantic salmon developmental stages (whole 186°d and 258°d embryos, alevin and fry). Vgll3a was expressed at low levels during all four developmental stages with the highest expression in alevin ( Figure 1). Six6a expression level was also overall low and declined from the embryonic stages toward fry ( Figure 1). Both akap11 paralogs were expressed at low levels in all stages although akap11b had clearly higher expression throughout development than akap11a (Figure 1). On the contrary, expression of vgll3b and six6b paralogs that do not associate with age-at-maturity was below background count level during all studied developmental stages (Figure 1). This indicates that the paralogs associated with maturation timing are expressed throughout early salmon development.
Genes with expression patterns correlating with vgll3a, akap11a and six6a during early developmental stages In order to identify genes that may be involved in the same molecular networks along with the age-at-maturity-associated genes vgll3a, akap11a and six6a, we determined the expression levels of 205 additional genes (Table S2) previously shown to have related functions or suggested as members of pathways linked with puberty and adiposity (see Methods). Pairwise correlations of these genes with vgll3a, akap11a and six6a were then estimated across the expression levels in samples from the four above-mentioned early developmental stages. Models for seven correlations did not converge due to singularities. Twenty-six correlations yielded results with P , 0.01 (range = 0.0001-0.009). Interestingly, expression of vgll3a correlated positively with one of the other age-at-maturity genes, akap11a, as Figure 1 mRNA expression patterns of the ageat-maturity-associated genes vgll3a, akap11a and six6a, and their paralogs vgll3b, akap11b and six6b during early Atlantic salmon developmental stages. Twenty was set as a count baseline (dashed line) above which counts were considered as detected. The data are represented as individual (dot) and mean (line) counts.
well as with arhgap6e and negatively with rd3l and yap1 (Figure 2). Akap11a correlated positively, in addition to vgll3a, with nr1i2 and negatively with foxp3a, lats1a, sox9d, rd3l and yap1 (Figure 2). Because rd3l and yap1 correlated with both vgll3a and akap11a, the correlation between the two was also estimated, and, indeed, a positive correlation was detected between them ( Figure 2). Six6a correlated positively with twelve genes related to, for example, the HPG axis signaling, eye development and PKA function, and negatively with vdrab and egr1d ( Figure 2). mRNA expression level differences between early developmental stages Age-at-maturity-associated genes vgll3a, akap11a and six6a, and the genes correlating with them were included in the model to estimate transcription expression differences across early developmental stages. Modeling revealed significant gene (F 23, 95.2 = 877.5, P , 0.001), stage (F 3, 13.4 = 31.4, P , 0.001) and their interaction effects (F 69, 140.9 = 25.8, P , 0.001). Controlled for the FDR, pairwise post-hoc comparisons indicated that all of the previously identified genes that correlated with the age-at-maturity-associated genes, but nr1i2 and rtn1, had significant expression changes during early developmental stages. Vgll3a and rd3l were upregulated in alevin compared to 258°d embryo ( Figure 3A). Expression of vgll3a, and also arhgap6e and yap1, was then decreased in fry compared to alevin ( Figure 3A). Akap11a was upregulated in fry compared to alevin ( Figure 3B). In contrast, significant expression changes between developmental stages in genes correlating negatively with akap11a tended to be downregulations, e.g. foxp3a in 258°d embryo, sox9d in alevin, and lats1a, sox9d and yap1 in fry, compared to the previous developmental stage ( Figure 3B). Six6a was downregulated in alevin and fry, compared to the previous stage ( Figure 3C). Similar trend of decreased expression throughout early development was observed with aesb, aesc, otx2b, prkar2ad, six3b, slc38a6, stk3a and tead1b that correlated positively with six6a, but also with egr1d ( Figure 3C). In addition, in alevin, prkar2ad and tead1b expression increased to a higher level than in 186°d embryo ( Figure 3C). In contrast, cyp26b1b was upregulated in 258°d embryo, prkar1aa in alevin and fry, and prkar1bb and vdrab in fry, compared to the previous stage ( Figure 3D).

Vgll3, akap11 and six6 expression profiles in parr and fry tissues
In addition to whole embryos, alevin and fry, we characterized mRNA expression patterns of vgll3, akap11 and six6 paralogs in 14 different tissues of three-year-old Atlantic salmon parr and blood of fry. Vgll3a expression was detected in eight of the 15 tissues: fin (adipose and caudal), gill, heart, liver, muscle, pyloric caeca, spleen and testis ( Figure  4A). Expression of vgll3b was only detected in the gill and ovary, the latter having the highest observed tissue-specific expression of vgll3 paralogs ( Figure 4A). Expression of akap11a and akap11b was detected in all studied tissues, akap11b mostly at a higher level than akap11a ( Figure 4B). Interestingly, the highest akap11a expression level, higher than that of its paralog, was detected in the ovary ( Figure 4B). Six6a expression was observed in four (brain, eye, gill and testis) of the 15 tissues ( Figure 4C), while the six6b paralog was expressed only in eye ( Figure 4C).

DISCUSSION
We used a custom Nanostring nCounter mRNA expression panel designed for Atlantic salmon to examine expression profiles of the vgll3, akap11 and six6 genes that have been recently associated with sexual maturation timing in Atlantic salmon. Our results show that the age-atmaturity-associated paralogs vgll3a, akap11a and six6a, and also akap11b, but not vgll3b and six6b, are all indeed expressed throughout early embryonic and juvenile salmon development, highlighting the potential relevance of the maturation timing-associated paralogs in developmental processes. We were also able to shed light on which other genes are the potential functional partners of these age-at-maturityassociated genes, and thus provide further insights into the molecular pathways and mechanisms behind maturation. Transcript expression correlations of vgll3a, akap11a and six6a with the subset of 205 prechosen genes related to the functions of the age-at-maturity genes in Figure 2 Genes correlating with the age-at-maturityassociated genes during early Atlantic salmon development. The plot shows the 25 gene-pairs for all 24 genes (out of 178 tested) whose expression levels correlated with expression levels of vgll3a, akap11a or six6a with P , 0.01. The data are reported as correlation coefficient (r) 6 95 % confidence interval (CI; 2 Ã se). embryos, alevin and fry suggest that genes associated with salmon maturation timing are linked with cell fate commitment regulation.
Of the studied 205 genes, two (arhgap6e and akap11a) correlated positively with vgll3a, while another two (yap1 and rd3l) correlated negatively. Correlation of vgll3a with arhgap6e (rho GTPase-activating protein 6-like) is noteworthy as ARHGAP6 was shown to be the most downregulated gene in human keratinocytes after VGLL3 knockdown (Liang et al. 2017) emphasizing its relevance in Vgll3-induced signaling also in Atlantic salmon. ARHGAP6 negatively regulates RhoA (the Rho family GTPase), thus causing actin fiber depolymerization (Prakash et al. 2000). This may be important in the context of maturation-related cellular processes as actin cytoskeleton control is crucial in determining cell fate commitments, such as cell proliferation and differentiation (Halder et al. 2012). For example, actin cytoskeleton depolymerization Figure 3 mRNA expression levels of the age-at-maturity-associated genes and genes correlating with them during early Atlantic salmon development. Vgll3a and genes its expression correlates with (A), akap11a and genes its expression correlates with (B) and six6a and genes its expression correlates with (C-D) in early developmental stages (n = 4). The mixed-model estimates were log 2 -transformed for better visualization and log 2 (20) = 4.3 was defined as a detection threshold. The data are represented as mean 6 SE. Expression level differences were tested between all stages but significant comparisons are highlighted by asterisks only for subsequent stages Ã FDR-adjusted P , 0.05, ÃÃ P , 0.01, ÃÃÃ P , 0.001.
is required for cell growth arrest (Halder et al. 2012) and SOX9 transcriptional activity to induce chondrocyte-specific markers and thus chondrogenesis (Kumar and Lassar 2009). Accordingly, Vgll3 overexpression upregulates expression of Sox9 and other genes in chondroand osteogenesis, and downregulates the main genes in adipogenesis and vice versa (Halperin et al. 2013), and, herein, we suggest that Vgll3 could conduct that via actin cytoskeleton regulation. As completion of maturation in Atlantic salmon (Rowe et al. 1991) and another salmon species (Silverstein et al. 1998) is highly dependent on the level of fat storage, the expression status of vgll3 could determine whether to induce adipogenesis or not and, therefore, regulate the timing of maturation. One of the main known regulators of cell fate is the Hippo signaling pathway member YAP (yes associated protein), a transcriptional cofactor that regulates cell proliferation and differentiation based on its cellular location and actin fiber status (Wada et al. 2011;Halder et al. 2012;Seo and Kim 2018) (Figure 5). However, there are contradictory results regarding the role of YAP in cell differentiation decisions (Seo et al. 2013;Deel et al. 2015;Karystinou et al. 2015;Deng et al. 2016). Also, YAP and VGLL3 seem to have overlapping effects on cell fate determination and our suggested outcome of Vgll3 function contradicts with its known role as an inhibitor of adipogenesis (Halder et al. 2012;Halperin et al. 2013;Figeac et al. 2019) ( Figure 5). Nevertheless, the inverse correlation of vgll3a with yap1 in our data implies that these two cofactors could have somewhat opposite roles during development in different stages of cell differentiation processes, such as adipogenesis. And while Yap upregulation has been detected during mammalian puberty (Sen Sharma et al. 2019;Zhang et al. 2019), it is intriguing that vgll3a, not yap1, genetic variation is tightly associated with maturation timing in Atlantic salmon (Ayllon et al. 2015;Barson et al. 2015). Moreover, we observed that the expression levels of vgll3a and akap11a correlate positively, and that expression of akap11a further correlates negatively with yap1. Interestingly, it is known that protein kinase A (PKA), the regulatory subunit to which AKAP11 binds in order to confine the enzyme to discrete locations within the cell (Reinton et al. 2000), can activate the Hippo pathway. It does this by activating LATS kinases, either indirectly by inhibiting Rho GTPase causing actin cytoskeleton depolymerization (Yu et al. 2013) or directly in response to actin disassembly (Kim et al. 2013), which leads into inactivation of YAP and further results into inhibition of cell proliferation and induction of adipogenesis. Our results showing inverse correlation between akap11a and lats1a confirms that PKA-induced Lats regulation is indeed dependent on Akap11. The association of Akap11 with PKA-induced adipogenesis is further supported by our finding that akap11a expression correlated negatively with sox9d, the downregulation of which is required for adipogenesis (Wang and Sul 2009). The relationship between Vgll3 and PKA-regulating Akap11 appears complex, but the observation that expression of these two age-at-maturity-associated genes correlate with each other and the evidence that they take part in the Hippo pathway speaks strongly for cell fate regulation being an important player in maturation. Figure 5 outlines a hypothetical model for the regulation of the Hippo pathway by Vgll3 and Akap11 based on a combination of our results and earlier research.
As akap11a is the adjacent gene of vgll3a in Atlantic salmon with a distance of 57 kb in between, their transcriptional correlation also indicates that they may be under the same transcriptional regulation and hence may be co-evolving, which could be meaningful when genes share the same molecular pathway (Hurst et al. 2002;Hurst et al. 2004). Co-expression of adjacent gene pairs is shown to be conserved across eukaryotes (Arnone et al. 2012), also in zebrafish (Tsai et al. 2009). In addition to the occurrence of missense SNPs in both the vgll3a and akap11a coding regions, the most highly age-at-maturity-associated SNP in the Atlantic salmon vgll3a locus resides in a non-coding region between vgll3a and akap11a (Barson et al. 2015). This led us to hypothesize that these genes could potentially share a regulatory region between them as it has indeed been shown that a bidirectional promoter can control transcription of an adjacent gene pair (Trinklein et al. 2004). Another intriguing link between vgll3a, akap11a and yap1 is rd3l (retinal degeneration 3-like), which correlated negatively with vgll3a and akap11a and positively with yap1. Knowledge of the functional role(s) of rd3l was extremely limited prior to this study. Expression of both vgll3a and rd3l was significantly upregulated at the alevin stage suggesting that their overall expression is induced simultaneously during development, albeit being oppositely regulated at the individual level. High expression of vgll3a in alevin may indicate that Vgll3 starts to induce chondro-and osteogenetic pathways during this stage when lepidotrichia, the bony segmented fin rays, are formed in salmon (Gorodilov 1996). However, it remains to be studied if and how Rd3l possibly participates in cell fate regulation. Altogether, our data suggest that the functions of vgll3a and akap11a are linked and associate with cell fate determination by regulating actin cytoskeleton assembly.
SIX6 and two other homeobox transcription factor proteins, SIX3 and OTX2, regulate the HPG axis signaling in the hypothalamus and pituitary, and eye development (Oliver et al. 1995;Jean et al. 1999;Larder et al. 2011;Larder et al. 2013;Xie et al. 2015). In accordance with these studies, our results show high correlation of six6a with six3b and otx2b, and with aesb and aesc, which encode a corepressor for SIX6 and SIX3 (Zhu et al. 2002;López-Ríos et al. 2003). Further, our results indicate that embryonic stages are crucial time points for the HPG axis and eye development in salmon since six6a, six3b, otx2b, aesb and aesc were all expressed at their highest level during embryonic stages, after which they started to decline. This is supported by a study showing that Six6 expression is detected already early in the development in regions that comprise the eye field and prospective hypothalamus and pituitary gland (López-Ríos et al. 1999). In addition, we detected that expression of slc38a6, which possibly encodes a transporter for glutamine-glutamate cycle in the brain (Bagchi et al. 2014) both correlated and covaried temporally with six6a during early development and was visible in salmon parr brain ( Figure S1). Glutamate is the most prevalent neurotransmitter in the central nervous system, including vertebrate retina (Thoreson and Witkovsky 1999) and the HPG axis where it induces GnRH release (Brann and Mahesh 1995). In addition to slc38a6, rtn1, another gene related to neuroendocrine secretion, correlated with Figure 5 A model for the regulation of the Hippo pathway by Vgll3 and Akap11. By binding to Tead transcription factors, our results suggest that Vgll3 induces expression of Arhgap6 (rho GTPase-activating protein 6) either directly or via other signaling partners. Arhgap6 inhibits RhoA (the Rho family GTPase) activity and thus causes actin fiber depolymerization, which further leads into Lats1/2 activation. Lats kinases phosphorylate the Yap/Taz complex excluding it from the nucleus and causing it to be degraded in the cytoplasm. In addition, PKA (protein kinase A) assisted by Akap11 stimulates Hippo signaling by both inducing and responding to actin depolymerization by activating Lats1/2 kinases either indirectly via RhoA inhibition or directly, respectively. Inactivation of Yap/Taz results in endothelial or epithelial cells to have growth arrest or apoptosis and mesenchymal stem cells (MSCs) to differentiate into chondrocytes or adipocytes. Instead, when the Hippo pathway is inactivated by actin polymerization, Yap/Taz accumulates in the nucleus, binds Teads and induces gene transcription leading to cell growth, proliferation and migration, and osteocyte differentiation. Yap can regulate its own activation via negative feedback loop by upregulating at least Arhgap18 and Arghap29 which suppress RhoA activity. Based on the results of the current study and Halder et al. 2012;Kim et al. 2013;Yu et al. 2013;Karystinou et al. 2015;Porazinski et al. 2015;Qiao et al. 2017;Seo and Kim 2018. Ã Correlation of vgll3a mRNA expression with that of arhgap6e and akap11a (positive), and yap1 (negative).°Correlation of akap11a mRNA expression with that of vgll3a (positive), and lats1a and yap1 (negative).
six6a. Rtn1 is highly expressed in the brain, which our results confirmed in salmon parr ( Figure S1), where it is considered to be a marker for neuronal differentiation (Yan et al. 2006). Interestingly, both slc38a6 and rtn1 are located on the same salmon chromosomal region as six6a. Specifically, slc38a6 locates 123 kb downstream and rtn1 165 kb upstream of six6a, suggesting that these three genes, which all have roles in neuroendocrine secretion, may be co-regulated, possibly by chromatin modification or folding (Hurst et al. 2004). Six6a also correlated and covaried temporally during early development with two Hippo signaling pathway genes, stk3a and tead1b. These genes encode Mst2 kinase needed to inactivate Yap and a transcription factor partner for Vgll3 and Yap, respectively. Accordingly, MST1/2 kinase genes, tead3 and also other major Hippo pathway members have been linked with maturation process both in Atlantic salmon (Christensen et al. 2017;Kjaerner-Semb et al. 2018) and other species (Lyu et al. 2016;Sen Sharma et al. 2019;Zhang et al. 2019). Additionally, six6a correlated with three PKA regulatory subunit genes. It is known that Hippo signaling regulates eye development (Zhu et al. 2018b) and that it is activated in embryonic and postnatal pituitary gland (Lodge et al. 2016). However, to our knowledge, there are no studies so far linking Six6/Six3 and the Hippo pathway. Instead, it has been shown that SIX3, in the embryonic brain (Lagutin et al. 2003) and in eye together with SIX6 (Diacou et al. 2018), and Hippo signaling (Varelas et al. 2010) repress Wnt signaling, another pathway regulating cell fate (Teo and Kahn 2010). Based on our results, it is intriguing that all three Atlantic salmon age-at-maturity-associated genes seem to associate with cell fate determination, especially with the Hippo pathway.
In addition to gaining clues about the early development transcription profiles and molecular mechanisms of vgll3a, akap11a and six6a, we scrutinized tissue-specific expression patterns of these genes and their paralogs in several tissues of three-year-old Atlantic salmon parr. Our results let us to hypothesize that Vgll3 influences maturation timing by regulating both adipogenesis and gametogenesis. Of the studied samples, vgll3a expression was found in fin, gill, muscle, heart, spleen, liver, pyloric caeca -an organ aiding to absorb nutrients specifically in fish -and testis, whereas vgll3b expression was restricted to ovary and gill. This, especially sex-dependent gonad-expression pattern with vgll3a expressed in testis, but vgll3b expressed in ovary, provides evidence for vgll3 paralog sub-functionalization in salmon. Our results are mostly in line with studies performed in other vertebrates showing vgll3 expression in, for example, skeletal muscle, heart, gill, liver, spleen, pyloric caeca and gonads, but some differences occur (Mielcarek et al. 2009;Faucheux et al. 2010;Kjaerner-Semb et al. 2018;Figeac et al. 2019;Simon et al. 2016). Knowledge of vgll3 paralog expression patterns helps to resolve how Vgll3 conducts its function during the maturation process. Based on our results, Vgll3 may participate in inhibiting cell proliferation and organ growth via actin cytoskeleton regulation. This is supported by the study of Kjaerner-Semb et al. (2018) showing that vgll3a is expressed in Sertoli cells and downregulated along with some Hippo pathway genes in mature salmon testis, thereby potentially inhibiting the Hippo signaling pathway and inducing subsequent proliferation of Sertoli cells, which are required to provide support for the developing germ cells (Griswold 1995). However, their finding of vgll3a being upregulated in granulosa cells of early and late puberty ovary (Kjaerner-Semb et al. 2018) does not support the same mechanism in ovary. This, combined with results of the current study, may indicate that gonad development during maturation is differently regulated in the two sexes. In other tissues, Vgll3 may also function as a regulator for cell fate commitment. According to a previous study (Halperin et al. 2013) and our tissue and early development mRNA expression results, Vgll3 could induce chondro-and osteogenic pathways in fin and gill, which would be needed in growing fish. However, in muscle, where we detected very low vgll3a expression, it has been shown that although VGLL3 induces myogenesis it is expressed at low levels in healthy muscle (Figeac et al. 2019). In line with this, our finding that adipose tissue lacks both vgll3 paralog expression is reasonable, as it is known that knockdown of Vgll3 expression promotes pre-adipocytes to differentiate into mature adipocytes (Halperin et al. 2013). In other words, also in salmon, decreased expression of vgll3a may be required to induce adipogenesis, which is critical to ensure enough fat-derived energy for sexual maturation (Rowe et al. 1991;Silverstein et al. 1998).
Unlike the vgll3 genes, both akap11 paralogs were expressed in all studied tissues, which is consistent with Reinton et al. (2000), emphasizing the relevance of Akap11 in basic cellular functions. However, akap11b was expressed at a higher level than akap11a in most tissues, suggesting its higher functional activity. Our results in salmon gonads confirm the ubiquitous nature of akap11 expression but further suggests somewhat specialized functions for each paralog. In contrast to vgll3 whose age-at-maturity-associated paralog was expressed in testis and other paralog in ovary, both akap11 paralogs were expressed in both ovary and testis. However, the age-at-maturity-associated akap11a and several paralogs encoding the regulatory subunits of PKA were more highly expressed in ovary, whereas akap11b expression was higher in testis. This differs somewhat from the study by Reinton et al. (2000) in that AKAP11 mRNA expression in the human ovary was extremely low. So far, it is known that as important as in sustaining interphase during mitosis (Grieco et al. 1996), PKA activation is required to maintain meiotic arrest of oocyte (Kovo et al. 2006). Since the studied salmon ovaries were immature, high akap11a and PKA regulatory subunit mRNA expression could result from the need of PKA to sustain oocytes in meiotic arrest. In testis, however, the AKAP11/PKA complex is suggested to have a dual role in meiosis during spermatogenesis and in motility of mature sperm (Reinton et al. 2000). Overall, our results propose that, in Atlantic salmon, akap11 has evolved paralogspecific functional differences, especially in gonads.
Expression of six6a was detected in brain, eye, testis and gill. To our knowledge, expression in gill has not been reported earlier. Instead, expression of six6b was detected only in eye, which could imply that the expression patterns observed for the age-at-maturity-associated paralog have evolved and expanded to the HPG axis-related tissues and, thus, reproduction. Although mRNA expression of SIX6 has been earlier detected in testis (Aijaz et al. 2005), its specific function remains unknown. Another transcript expressed in brain, but whose expression we detected also in testis, is cgba ( Figure S1). This gene encodes the beta-subunit of Fsh (follicle-stimulating hormone) gonadotropin, which is normally expressed in pituitary and regulated by SIX6 (Xie et al. 2015). These findings hint that Six6 may potentially regulate the expression of the testicular Fsh beta-subunit and thus testis development. This is supported by an earlier observation that FSH beta-subunit mRNA is expressed in mouse testis, where it is suggested to play paracrine or autocrine role in the regulation of testicular function (Markkula et al. 1995), such as Sertoli and germ cell development (Meachem et al. 2005).
We conducted one of the rare follow-up studies of a GWAS association in a wild animal species with a view to better understanding the molecular mechanisms of the previously discovered genotype-phenotype association. Our temporal assessment of mRNA expression of Atlantic salmon age-at-maturity-associated vgll3a, akap11a and six6a during five developmental stages revealed differently regulated expression in a development stage-and tissue-specific manner. Co-expression analysis of the age-at-maturity-associated genes and 205 other pre-selected genes indicated co-regulation of vgll3a and akap11a, and a novel role for Vgll3 in regulating actin cytoskeleton assembly -a process required in cell proliferation and differentiation -and that this regulation could be assisted by Akap11-directed PKA function. In addition, we were able to confirm the same expressional pattern in salmon that Six6 and its partners have in the HPG axis and eye development regulation in other vertebrates. Moreover, we propose that Six6 may associate more broadly with neuroendrocrine regulation in the HPG axis and have a direct role in testis function. Further, our data provide the first evidence that both vgll3 and six6 paralogs have sub-functionalized roles in different tissues. Overall, we conclude that Vgll3, Akap11 and Six6 may contribute to influencing Atlantic salmon maturation timing via affecting on adipogenesis and gametogenesis by regulating cell fate commitment and the HPG axis. Further studies are required to determine the more specific cellular molecular mechanisms of Vgll3, Akap11 and Six6 in sexual maturation processes. This work provides important information for guiding such work, also in other organisms.

ACKNOWLEDGMENTS
We kindly thank Päivi Anttonen and Risto Kannel for help collecting the samples and providing rearing temperature data. We also acknowledge Katja Salminen for help extracting RNA, Victoria Pritchard for suggesting genes for the mRNA expression panel and the staff of the Institute of Biotechnology at University of Helsinki for processing the mRNA samples. Marion Sinclair-Waters is thanked for providing relative chromosomal locations of genes of interest.