Octopus ocellatus samples were obtained from a commercial hatchery. This research was conducted in accordance with the protocols of the Institutional Animal Care and Use Committee of the Ludong University (protocol number LDU-IRB20210308NXY) and the China Government Principles for the Utilization and Care of Invertebrate Animals Used in Testing, Research, and Training (State Science and Technology Commission of the People’s Republic of China for No. 2, October 31, 1988¹).

We caught wild O. ocellatus parents in Rizhao sea area to prepare for this experiment, the latitude and longitude of animals sampling site are 34.91005177008973 and 119.63820016113283. The eggs spawned by parents were collected after temporary feeding, and these were divided into two groups (Pro and Unp). Among them, the eggs in Pro group continued to be brooded by their parents, and the eggs in Unp group were hatched in flowing seawater without parent brood. Then, the primary incubation larvae were temporarily cultured in floating seawater for 24 h. The temperature of seawater was stable at the optimum temperature of eggs hatching (19–20.8°C) in 29 days of hatching. The larvae at 0h, 4h, 12h, and 24h after hatching were collected and stored in liquid nitrogen container until RNA extraction using TRIzol method.

We randomly selected nine primary incubation larvae at each time point in two groups for RNA extraction: brood-care larvae grow for 0h (Pro-C), brood-care larvae grow for 4h (Pro-4h), brood-care larvae grow for 12h (Pro-12h), and brood-care larvae grow for 24h (Pro-24h); brood-careless larvae grow for 0h (Unp-C), brood-careless larvae grow for 4h (Unp-4h), brood-careless larvae grow for 12h (Unp-12h), and brood-careless larvae grow for 24h (Unp-24h). Three of the nine larvae in a group at each time point were randomly selected, and their equal molar mass of RNA was concentrated in a replicate as a template for constructing the transcriptome library; the same method was used to pooled equal molar mass RNA from the three larvae into the second replicate, and the equal molar mass RNA of the remaining larvae was concentrated in the third replicate. We subsequent qRT-PCR verification using remaining RNA that has been preserved.

The method of Li et al. (2017) description was used to construct the library. And Illumina Hiseq 4000 platform was used to sequence samples.

We splice Trinity with RSEM to obtain reference sequences and map the clean reads to these sequences. Then, a great quantity of read counts of each sample mapped to each gene was obtained. To understand the expression and abundance of genes, we transform read counts into FPKM. The result showed that there was a positive correlation between FPKM and the expression level of samples. The unigenes were annotated by searching the sequences against the NR, NT, GO, SwissProt, and KOG databases using BLASTX with a cut-off of Evalue ≤ 1e-5. Meanwhile, we annotated unigenes into Pfam database with Hmmer 3.0 package (Evalue ≤ 0.01).

Differentially expressed genes (DEGs) were screened out using the DESeq2 package for R, and the filter parameters were | log2 Fold Change| ≥ 4 and q-value ≤ 0.01. We identified DEGs distribution and multiple GO terms by analyzing DEGs with GO functional enrichment analysis and conduct the DEGs statistical analyses. We used KEGG enrichment analysis to further study the specific functions of DEGs, and the metabolism-related pathways involved in DEGs were annotated. Finally, the signaling pathways with significant enrichment of DEGs were selected, and we count the species and number of DEGs after selecting.

STRING v11.0 with default parameters was used to construct PPI networks to study relationships between functions of metabolism-related genes (Damian et al., 2018).

We verified the accuracy of our RNA-Seq results using quantitative real-time polymerase chain reaction (qRT-PCR) to validate 10 selected genes. In this experiment, each group contained three biological replicates, and Primer Premier 5.0 software was used to project gene-specific primers. We screened 10 genes for verification, and their names and primer sequences are listed in Table 1. The stability of three genes including 18S, β-actin, and GAPDH were evaluated in different tissues and embryo development stages of O. ocellatus. By comparison, we used O. ocellatusβ-actin as an endogenous control based on its expression level which tended to be stable. This experiment is based on Li et al.’s (2019) description method for qRT-PCR.

Brood-care behavior affects the hatching rate and survival rate of O. ocellatus larvae. The survival rate of brood-careless eggs was significantly lower than that of normal hatching larvae, and the hatching rate was reduced by 30–40%. However, there were no significant differences in morphology and development of larvae between Pro and Unp.

Samples at four time points in two groups were sequenced by RNA-Seq. Table 2 shows the detailed results of sequencing. Raw sequencing reads were submitted to Sequence Read Archive in NCBI; the SRA accession numbers were SRR15204591, SRR15204592, SRR15204593, SRR15204594, SRR15204595, SRR15204596, SRR15204597, SRR15204598, SRR15204599, SRR15204600, SRR15204601, SRR15204602, SRR15927331, SRR15927332, SRR15927333, SRR15927334, SRR15927336, SRR15927337, SRR15927338, SRR15927339, SRR15927340, SRR15927341, SRR15927342, and SRR15927343².

Through differential expression analysis, we found 229, 583, 659, and 2,328 DEGs in Pro group at 0h, 4h, 12h, and 24h after hatching compared with Unp group. Among these, 101 DEGs were up-regulated, and 128 DEGs were down-regulated at 0h after hatching; 298 DEGs were up-regulated, and 285 DEGs were down-regulated at 4h after hatching; 323 DEGs were up-regulated, and 336 DEGs were down-regulated at 12 h after hatching; 1,657 DEGs were up-regulated, and 671 DEGs were down-regulated at 12 h after hatching (Figure 1). Venn diagram helps us find the union of DEGs at three time points (Figure 2). The clustering distribution of DEGs is shown in Figure 3 intuitively.

We identified 269 level-3 terms based on GO enrichment analysis including 151 biological process subclasses, 72 molecular function subclasses, and 46 cellular component subclasses. And the top 10 level-3 terms were shown in Figure 4. Meanwhile, The KEGG signaling pathways enriched with DEGs was analyzed to help us further understand the gene functions. We enriched 179 level-2 KEGG classes pathways containing 1,777 selected DEGs (Figure 5), and 10 significantly enriched metabolism-related pathways (Table 3) were identified in our study.

Proteins exist in all organisms as the main components of cells and tissues. All significant parts of the organism need protein participation, which is the main undertaker of life activities. We can identify key genes in the metabolic gene population by constructing protein–protein interaction networks. In our research, we used 49 gene protein sequences in the signaling pathways in Table 3 to construct the protein–protein interaction networks. The networks were shown in Figure 6. We found that the ITPR1 gene had the highest number of protein-protein interactions, and other genes were interacted with more than 10 proteins. Relevant parameter information of them was listed in Table 4. There were 44 nodes and 197 edges in the network, and each node was interacted with an average of 8.95 nodes. The clustering coefficient was 0.5, and the p-value of this network was ≤ 1.0e-16.

The interaction relationships between key metabolism-related DEGs were studied in our research. We identified 10 genes (Table 5) with high protein interactions or involved in multiple signaling pathways based on KEGG enrichment and PPI networks analyses, and then explored their interaction relationships. We divided these 10 key DEGs into five categories: solute carrier family, purine metabolism signaling pathway, cAMP signaling pathway, protein digestion and absorption signaling pathway, and other significant metabolism-related DEGs. The family and signaling pathways were gone hand in hand to metabolic processes, analyzing them will promote us to further understand the metabolic mechanisms of differences of larval growth between Pro group and Unp group.

We detected the expression of 10 metabolism-related DEGs at each time point in the two groups using qRT-PCR, and verified the consistency of qRT-PCR and RNA-Seq expression. The results of qRT-PCR suggested that all DEGs measured were single products. We compared the gene expression profiles between qRT-PCR results and RNA-Seq results (Figure 7). And the result suggested that their results were significantly correlated, and showed the same trend pattern.

Due to the rich nutrition and excellent taste, O. ocellatus is gradually favored by people (Boletzky, 1975; Wei et al., 2018). As an economically important species employed in mariculture, the artificial breeding of O. ocellatus is an issue of concern. We understand that multiple octopuses have the brood-care behavior based on previous studies, and the behavior directly affects the growth and development of their larvae. Thus, understanding the effects of brood-care behavior is vital to O. ocellatus larval growth and artificial breeding. Our research screened 3,486 DEGs, and these genes are considered to have significant relationships with O. ocellatus larval metabolic processes. A heatmap showed that the genes in two groups have significant differences of the same point in time. It is further indicated that there are quite metabolic differences between the brood-care larvae and the brood-careless larvae. Then, we identified 10 metabolism-related KEGG signaling pathways, they interact with each other and share the metabolic function of organisms. Finally, we construct protein-protein interaction networks using 49 DEGs in these signaling pathways to further explore O. ocellatus larval metabolic mechanisms.

Many metabolism-related terms and pathways were yielded through GO and KEGG enrichment analyses, including cell adhesion term, protein phosphorylation term, cAMP signaling pathway, purine metabolism signaling pathway, and other key pathways and terms. These pathways and terms are directly related to metabolic processes, suggesting that there were several metabolism-related DEGs between larvae in two groups, and thus resulted in complex metabolic mechanisms of larvae. We comprehensively analyzed these terms and signaling pathways. The results will help us to understand the metabolic processes and molecular mechanisms of O. ocellatus under different incubation conditions within 24h of incubation, so as to lay a foundation for the artificial breeding of O. ocellatus.

Protein is an essential component of an organism and is closely related to metabolic processes. Studying their interaction relationships helps us to enrich our understanding of biological metabolic mechanisms We construct protein–protein interaction networks based on 49 key genes contained in metabolism-related signaling pathways. These proteins interacted significantly better than randomly selected groups of proteins of similar size. This conclusion suggests that the above proteins may function as a group. We thus speculated that proteins interact with more other proteins in the network were hub proteins in metabolism. The genes corresponding to hub proteins are further studied.

In our research, we analyze the differences in metabolic mechanisms between larvae under different hatching conditions using transcriptome profiling. The results help us to further understand the effect of O. ocellatus brood-care behavior on metabolic mechanisms. Ten hub genes were speculated and screened involved in multiple protein interactions or metabolism-related signaling pathways, and these hub genes, pathways, and the larval metabolic differences between two groups at the same time point were investigated.

Purine metabolism is one of the key ways for aquatic animals to detoxify and obtain energy, which plays a major part in maintaining the metabolic stability of organisms. In the liver of marine organisms, purines are converted by degrading enzymes into urea and glyoxylate, and then excreted by other metabolic pathways (Noguchi et al., 1982). In recent studies, we found that there is a urease that can decompose urea into NH3 in marine invertebrates. NH3 is very soluble in water, organisms can excrete urea in the form of NH3 and CO2 through purine metabolism to prevent poisoning caused by excessive urea in the body (Yeldandi et al., 1996). In addition to detoxification, purine metabolism will release ATP, ADP, and other derivatives, providing the energy necessary for biological life activities (Furchgott, 1966; Sakuraba et al., 1996). In our research, key metabolism-related genes such as NPR1 and AMPD2 were enriched in this pathway. Compared with the Unp group, NPR1 was continuously up-regulated with the growth of O. ocellatus larvae, AMPD2 was down-regulated slightly in the first 4 h of growth, and up-regulated significantly in the last 20 h. These results indicate that brood-care larvae have higher detox and energy production abilities than brood-careless larvae, and that they have higher metabolic capacity; cAMP is an essential substance in biological cells. It regulates many biological functions, such as proliferation, differentiation, and apoptosis of cells. Meanwhile, cAMP also plays a role in regulating the balance of energy and nutrition. cAMP is widely distributed in many organelles, such as nucleus and cytoplasm. In mitochondria, cAMP can stabilize the function of mitochondria, and promote the growth of O. ocellatus larvae by regulating the stress response of cells and mitochondrial activity (Zhang et al., 2016). cAMP signaling pathway acts a key part in regulating the metabolism of saccharides and fat. It maintains the metabolic homeostasis by regulating the content of various saccharides. At the same time, it can regulate the formation of fat and the distribution of lipid in adipose tissue, control the differentiation of adipocytes, and maintain the fat content in vivo (Ravnskjaer et al., 2015). The maladjustment of cAMP signaling pathway may lead to cell metabolic disorders and the occurrence of various diseases (Tang et al., 2021). In our research, a large number of DEGs were significantly enriched in this pathway. And most DEGs were up-regulated in Pro group compared with Unp group, suggesting that the metabolic mechanisms of the brood-care larvae were more complex, and the ability to regulate this metabolic process was stronger. This phenomenon explains that the O. ocellatus larvae in Pro group had a stronger metabolic capacity and higher survival rate compared with those in Unp group; Protein digestion and absorption signaling pathway regulates the decomposition of proteins by digestive organs to produce favorable carbohydrates and energy (Freeman and Kim, 1978). This pathway was significantly enriched in our research, indicating that the digestive organs of O. ocellatus larvae were active, and metabolized rapidly. Moreover, we found multiple metabolism-related genes enriched in this pathway were significantly up-regulated compared with brood-careless larvae, including COL22A1, COL12A1, and other key genes. However, COL9A1 were down-regulated in our research, the specific reasons need to additional explore. All of these pathways are important pathways for the regulation of biological metabolism, and they are significantly enriched. Compared with the Unp group, most of the genes that are beneficial to biological metabolism enriched in above three important pathways were up-regulated, which suggests that the metabolic efficiency of O. ocellatus larvae in the Pro group was higher and more stable. And it was beneficial to the growth and development of O. ocellatus larvae.